Введение

Очень кратко рассмотрев основы механики в предыдущей статье, перейдем к практике, ибо даже той краткой теории что была рассмотрена хватит с головой.



Итак, задача:

Камень бросают вертикально, без начальной скорости с высоты h = 100 м. Пренебрегая сопротивлением воздуха, определить закон движения камня, как функцию высоты камня над поверхностью Земли от времени. Ускорение свободного падения принять равным 10 м/с²

Простая задачка? Да элементарная, имеющее аналитическое решение, которое легко напишет мало-мальски грамотный школьник. Но эта простая задача послужит нам весьма показательным примером

1. Формализация задачи

То, с чего начинается любое моделирование. Под формализацией понимают получение математических выражений, описывающих происходящий процесс. При этом формулируются допущения: перечень упрощений модели за счет факторов, влиянием которых можно пренебречь.

Для этой задачи применимо допущение, согласно которому камень можно считать материальной точкой. К этой точке приложена одна единственная сила — сила тяжести, так мы используем допущение об отсутствии сопротивления воздуха.

Следующим допущением будет полагать что Земля плоская, так как высота, с которой мы бросаем камень пренебрежимо мала в сравнении с размерами планеты, а значит кривизной её поверхности можно пренебречь. Тогда сила тяжести может считаться постоянной, направленной перпендикулярно поверхности Земли и равной по модулю

$$

где g — ускорение свободного падения у поверхности Земли. Теперь пришло время составить уравнения движения камня. Помните эти уравнения

$$

Левая их часть нас пока не интересует, а вот в правой стоят суммы проекций сил, приложенных к точке на оси координат. Пусть оси x и y располагаются на поверхности, а ось z направлена вверх перпендикулярно ей. Сила одна единственная, её проекции на оси x и y равны нулю, а на ось z проекция отрицательна, так как сила направлена против направления оси, то есть

$$

Масса камня, очевидно не равна нулю, значит можно спокойно поделить обе части получившихся уравнений на неё

$$

Не буду занудствовать, доказывая что движение камня будет происходить строго вертикально, хотя это нужно сделать с формальной точки зрения. Ноль в правой части первых двух уравнений совершенно не означает невозможность движения вдоль этих осей — вспоминаем первый закон Ньютона. На этом я остановлюсь в следующей статье более подробно, а пока справедливо положим одномерность движения, выписав окончательное дифференциальное уравнение

$$

То что у нас получилось не много не мало — математическая модель процесса происходящего в задаче. Пафосно, да?

Нет. Анализируя это уравнение мы делаем вывод, например, что масса камня не оказывает влияния на закон его движения, ведь массы в этом уравнении нет. Видите, даже не решив уравнения, мы уже формально доказали справедливость опыта с пером и кусочком свинца в вакууме, который любят показывать в школе (а некоторые повторили его на Луне).

Аналитическое решение получить просто, даже не буду заморачиваться, оно такое

$$

А вот как решить это численно? И что это вообще такое — «численно»?

2. Численное интегрирование дифференциального уравнения первого порядка

Какой такой первый порядок? Я же говорил в прошлый раз, что уравнения движения имеют второй порядок. Всё правильно, но большинство методов решения диффур на компьютере умеют решать только уравнения первого порядка. Есть методы прямого интегрирования уравнений второго порядка (например метод Верле), но о них не сейчас.

Во-первых, это уравнение относится к такому типу, что допускает понижение порядка. Правая часть не зависит от неизвестной функции (там нет z), поэтому вспоминаем, что

$$

проекция ускорения на ось z равна первой производной проекции скорости на ту же оcь z. Ну классно, тогда

$$

вот вам и уравнение первого порядка. Не всегда этот номер проходит (не буду я сейчас про форму Коши!), но в данном случае всё в порядке. Будем искать не координату а скорость точки. Что дальше-то? А дальше

$$

ведь производная, мы же знаем, это отношение бесконечно малого приращения функции (скорости) в к вызвавшему его бесконечно малому приращению аргумента (время). Возьмем очень маленький промежуток времени $$, настолько небольшой что можно считать

$$

Что получается? А вот что

$$

Мы получили приращение скорости. Отрицательное приращение. Как это так, камень, падая вниз будет разгонятся же! Да, будет. Его скорость, вектор его скорости, будет направлен вниз. А значит проекция этого вектора на ось z будет отрицательной. Всё правильно, мы получаем растущую по абсолютному значению проекцию вектора, направленного вниз. Мы знаем, начальное значение скорости — ноль, а значит

$$

Пользуясь тем что мы можем вычислить приращение скорости, посчитаем, какова будет скорость скажем через 0.1 секунды

$$

а ещё через 0.1 секунды

$$

и ещё через 0.1 секунды

$$

Хм, так мы можем продолжать довольно долго, но ограничимся промежутком времени в одну секунду

Время, с	Скорость, м/с
0.0	0.0
0.1	-1.0
0.2	-2.0
0.3	-3.0
0.4	-4.0
0.5	-5.0
0.6	-6.0
0.7	-7.0
0.8	-8.0
0.9	-9.0
1.0	-10.0

То есть, воспользовавшись формулой

$$

мы получили зависимость скорости точки от времени. А всего-то нужно взять значение скорости в текущий момент времени, и добавить к нему приращение, которое скорость получит в новый момент времени, отстоящий от текущего на $$ секунд. Приращение времени называется здесь шагом интегрирования. А приращение мы вычисляем как значение производной искомой функции в текущий момент времени умноженное на шаг. Просто? Да просто конечно. И та формула, которую я написал, имеет название название — явный метод Эйлера для численного решения дифференциальных уравнений. Это так называемая рекуррентная формула, когда новое значение вычисляемой величины зависит от её предыдущего значения.

А что же с высотой точки над Землей? Да всё аналогично, смотрите

$$

ведь проекция скорости есть производная от соответствующей координаты. Применим формулу метода Эйлера для этого уравнения, ведь скорость мы уже знаем

$$

и по этой формуле добавим в нашу таблицу ещё одну колонку

$$

и так далее

Время, с	Скорость, м/с	Высота, м
0.0	0.0	100.0
0.1	-1.0	100.0
0.2	-2.0	99.9
0.3	-3.0	99.7
0.4	-4.0	99.4
0.5	-5.0	99.0
0.6	-6.0	98.5
0.7	-7.0	97.9
0.8	-8.0	97.2
0.9	-9.0	96.4
1.0	-10.0	95.5

Хм, ну, во-первых, заметно, что высота меняется у нас уже неравномерно, так как скорость со временем меняется. Теперь наша производная сама зависит от времени. Но уже на первом шаге, мы замечаем неладное — скорость уже есть, а вот высота по прежнему 100 метров. Как так?

Это вышло потому, что на каждом шаге мы полагаем производную (скорость) постоянной. Метод не дает информации о том, что происходит с решением между шагами. Соответственно накапливается погрешность, сравним полученное решение с точным

$$

Время, с	Скорость, м/с	Высота, м	Точное решение, м
0.0	0.0	100.0	100.00
0.1	-1.0	100.0	99.95
0.2	-2.0	99.9	99.80
0.3	-3.0	99.7	99.55
0.4	-4.0	99.4	99.20
0.5	-5.0	99.0	98.75
0.6	-6.0	98.5	98.20
0.7	-7.0	97.9	97.55
0.8	-8.0	97.2	96.80
0.9	-9.0	96.4	95.95
1.0	-10.0	95.5	95.00

Да, наш камень как будто зависает в воздухе. Численное решение отстает от аналитического, и чем дальше мы считаем, тем выше погрешность счета. Погрешность накапливается, так как на каждом шаге мы берем всё более и более грубое приближение. Что делать?

Во первых, можно уменьшить шаг. Скажем в 10 раз, пусть $$ секунды

Время, с	Скорость, м/с	Высота, м	Точное решение, м
0.0	0.0	100.0	100.00
0.1	-1.0	99.96	99.95
0.2	-2.0	99.81	99.80
0.3	-3.0	99.57	99.55
0.4	-4.0	99.22	99.20
0.5	-5.0	98.78	98.75
0.6	-6.0	98.23	98.20
0.7	-7.0	97.59	97.55
0.8	-8.0	96.84	96.80
0.9	-9.0	96.00	95.95
1.0	-10.0	95.05	95.00

Уже лучше, погрешность в конце счета не превышает 0,05 метров, и это в 10 раз меньше предыдущего значения. Можно предположить, что уменьшив шаг ещё в 10 раз мы получим ещё более точное решение. Я схитрил, выводя значения только для 10 точек с шагом 0.1, на самом деле, чтобы получить такую таблицу нужны уже 100 итерации а не 10. При шаге 0.001 потребуется уже тысяча итераций, а результат будет таким

Время, с	Скорость, м/с	Высота, м	Точное решение, м
0.0	0.0	100.0	100.00
0.1	-1.0	99.9505	99.95
0.2	-2.0	99.8010	99.80
0.3	-3.0	99.5515	99.55
0.4	-4.0	99.2020	99.20
0.5	-5.0	98.7525	98.75
0.6	-6.0	98.2030	98.20
0.7	-7.0	97.5535	97.55
0.8	-8.0	96.8040	96.80
0.9	-9.0	95.9545	95.95
1.0	-10.0	95.0050	95.00

Если вы попробовали выполнить эти расчеты в ручную, то понимаете теперь насколько они однообразны и трудоемки, если нужна высокая точность. Именно поэтому расцвет численного моделирования совпал с появлением компьютеров. Они как раз и нужны для того, чтобы быстро выполнять множество однообразных операций над числами.

Метод Эйлера самый простой из известных методов интегрирования дифференциальных уравнений. Из нашего простого примера видно, что погрешность метода прямо пропорциональна шагу интегрирования, и это действительно так. Такие методы называются методами 1-го порядка точности.

Точность расчетов даже на шаге 0.1 можно улучшить, если мы применим, скажем метод Рунге-Кутты 4-го порядка точности. Но это отдельная история.

Заключение

Задумайтесь… Мы рассмотрели очень простой пример. Мы даже не применяли компьютер, но уже понимаем принцип, по которому работают те самые мощные в мире суперкомпьютеры, что моделируют ранние этапы жизни Вселенной. Конечно, там всё устроено гораздо сложнее, но принцип лежит этот же самый.

Представьте себе, какой мощный инструмент вы получаете в свои руки. Эта последняя статья, где мы не будем применять компьютер. Я обещал Octave. В следующий раз будет именно он.