Моделирование динамических систем: численные методы решения ОДУ

  • Tutorial

Введение


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



Итак, задача:
Камень бросают вертикально, без начальной скорости с высоты h = 100 м. Пренебрегая сопротивлением воздуха, определить закон движения камня, как функцию высоты камня над поверхностью Земли от времени. Ускорение свободного падения принять равным 10 м/с2
Простая задачка? Да элементарная, имеющее аналитическое решение, которое легко напишет мало-мальски грамотный школьник. Но эта простая задача послужит нам весьма показательным примером

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


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

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

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

$P = m \, g$



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

$ \begin{cases} &m \, \ddot x = \sum F_{kx} \\ &m \, \ddot y = \sum F_{ky} \\ &m \, \ddot z = \sum F_{kz} \end{cases} $



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

$ \begin{cases} &m \, \ddot x = 0 \\ &m \, \ddot y = 0 \\ &m \, \ddot z = -m \, g \end{cases} $



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

$ \begin{cases} &\ddot x = 0 \\ &\ddot y = 0 \\ &\ddot z = -g \end{cases} $



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

$\ddot z = -g$



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

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

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

$z(t) = h - \frac{g \, t^2}{2}$



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

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



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

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

$\ddot z = \dot v_z$



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

$\dot v_z = -g$



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

$\frac{dv_z}{dt} = -g$



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

$\frac{dv_z}{dt} \approx \frac{\Delta v_z}{\Delta t}$



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

$\Delta v_z = -g \, \Delta t$



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

$v_z^{(0)} = 0$


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

$v_z^{(1)} = v_z^{(0)} - g \, \Delta t = 0 - 10 \cdot 0.1 = -1.0$


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

$v_z^{(2)} = v_z^{(1)} - g \, \Delta t = -1.0 - 10 \cdot 0.1 = -2.0$


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

$v_z^{(3)} = v_z^{(2)} - g \, \Delta t = -2.0 - 10 \cdot 0.1 = -3.0$


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

Время, с Скорость, м/с
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

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

$v_z^{(i+1)} = v_z^{(i)} - g \, \Delta t$



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

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

$\frac{dz}{dt} = v_z$



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

$z^{(i+1)} = z^{(i)} + v_z^{(i)} \, \Delta t$



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

$\begin{align} &z^{(1)} = z^{(0)} + v_z^{(0)} \, \Delta t = 100 + 0.0 \cdot 0.1 = 100 \\ &z^{(2)} = z^{(1)} + v_z^{(1)} \, \Delta t = 100 - 1.0 \cdot 0.1 = 99.9 \\ &z^{(3)} = z^{(2)} + v_z^{(2)} \, \Delta t = 99.9 - 2.0 \cdot 0.1 = 99.7 \\ \end{align}$



и так далее

Время, с Скорость, м/с Высота, м
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 метров. Как так?

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

$z(t) = 100 - 5 \, t^2$



Время, с Скорость, м/с Высота, м Точное решение, м
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 раз, пусть $\Delta t = 0.01$ секунды

Время, с Скорость, м/с Высота, м Точное решение, м
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. В следующий раз будет именно он.
Поделиться публикацией

Комментарии 11

    0
    Отличный пример «на пальцах» для начинающих!
      0
      Я как-то писал решение простейшего ОДУ, но при этом в уравнение еще входило показание одного датчика. Естественно, датчик шумит.
      Я сравнивал метод Эйлера, методы Рунге-Кутты и Адамса — Башфорта (из Википедии) разных порядков.
      Так вот метод Эйлера оказался лучше всех. Точнее, и, самое главное, использованные методы высокого порядка иногда расходились. Как раз из-за того, что в уравнении присутствует шумный сигнал датчика. Метод Эйлера отклонялся от решения, но всегда к нему возвращался.

      Есть ли какие-то простые численные методы, учитывающие неидеальность реальных сигналов?
        0
        Здесь лучше смотреть в сторону фильтров. Классика — фильтр Калмана.
        Но если датчик инерциальный, то по высоте (направлению g) ничего толкового не получится практически никогда.
          0
          Может это неявный Эйлер был?
            0
            что такое неявный Эйлер?
              0
              Когда уже сделают формулы в комментах…

              Явный метод Эйлера выглядит так

              y(i+1) = y(i) + f[y(i), t(i)] * dt

              то есть производную мы вычисляем, опираясь на текущее значение фазовой координаты. В неявном методе делают так

              y(i+1) = y(i) + f[y(i+1), t(i+1)] * dt

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

              y(i+1) — y(i) — f[y(i+1), t(i+1)] * dt = 0

              относительно y(i+1). Для систем оду это приводит к необходимости вычислять матрицу Якоби. Но устойчивость метода для жестких систем коллосальна
                0
                Аааа, вспомнил.
                Я неявные методы сразу отбросил, из-за того, что в решение входит метод Ньютона с неопределенным количеством итераций, а решение надо выдавать в реальном времени.
                Может, я в чем-то не прав.
          0
          Очень познавательно
          хех
          для студентов 1 курса
          Надеюсь, в последующих выпусках будут фигурировать стохастические методы, а не хеллоуворлд от мира численного моделирования:D
            0
            Это и есть хелло ворлд, так что вы пока что не по адресу
            0
            Про методы повышения точности без уменьшения шага будет статья?
              0
              Еще не читал, но после прочтения первой части добавил вторую в закладки. Спасибо.

              Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

              Самое читаемое