Численное интегрирование системы обыкновенных дифференциальных уравнений

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

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

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

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. В следующий раз будет именно он.

Численное интегрирование системы обыкновенных дифференциальных уравнений

Системой дифференциальных уравнений называется система вида

где x — независимый аргумент,

y_i — зависимая функция, ,

Функции y_i(x), при подстановке которой система уравнений обращается в тождество, называется решением системой дифференциальных уравнений.

Численные методы решения систем дифференциальных уравнений.

Модифицированный метод Эйлера.

Метод Рунге-Кутта четвертого порядка.

Дифференциальным уравнением второго порядка называется уравнение вида

F(x,y,у’,y»)=0

(1)

y»=f(x,y,y’).

(2)

Функция y(x), при подстановке которой уравнение обращается в тождество, называется решением дифференциального уравнения.

Численно ищется частное решение уравнения (2), которое удовлетворяет заданным начальным условиям, то есть решается задача Коши.

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

.

(3)

Функция f₂(x, y₁, y) в систему (3) введена формально для того, чтобы методы, которые будут показаны ниже, могли быть использованы для решения произвольной системы дифференциальных уравнений первого порядка. Рассмотрим несколько численных методов решения системы (3). Расчетные зависимости для i+1 шага интегрирования имеют следующий вид. Для решения системы из n уравнений расчетные формулы приведены выше. Для решения системы из двух уравнений расчетные формулы удобно записать без двойных индексов в следующем виде:

Метод Рунге-Кутта четвертого порядка.

где h — шаг интегрирования. Начальные условия при численном интегрировании учитываются на нулевом шаге: i=0, x=x₀, y₁=y₁₀, y=y₀.

Контрольное задание по зачетной работе.

Колебания с одной степенью свободы

Цель. Изучение численных методов решения дифференциальных уравнений второго порядка и систем дифференциальных уравнений первого порядка.

Задание. Численно и аналитически найти:

1. закон движения материальной точки на пружинке х(t),
2. закон изменения силы тока I(t) в колебательном контуре (RLC — цепи) для заданных в табл.1,2 режимов. Построить графики искомых функций.

Свободные незатухающие колебания

Затухающее колебательное движение

Предельное апериодическое движение

Вынужденное колебание без сопротивления

Вынужденное колебание без сопротивления, явление резонанса

Вынужденное колебание с линейным сопротивлением

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

Численное интегрирование системы обыкновенных дифференциальных уравнений

Особенности численного моделирования решений обыкновенных дифференциальных уравнений

Основные источники ошибок методов численного интегрирования

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

Таким образом, Вы видите, что за появление d _n «ответственны» два процесса. Первый — процесс возникновения ошибки на шаге. Второй — процесс преобразования уже возникшей ошибки на последующих шагах (накопление ошибок). Первый процесс — более простой, второй — более сложный.

В этом разделе будут проанализированы основные источники ошибок на примере метода Эйлера.

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

Проанализируем процесс численного моделирования в методе Эйлера. Будем пока предполагать, что все вычисления выполняются точно. Попытаемся выразить d _n+1 через d _n.

Учитывая, что y(x_n) = y_n + d _n, и используя теорему о среднем, получим

Первый член определяет преобразование на шаге уже имевшейся ошибки d _n, а второй член e_n – дополнительную ошибку данного шага (ошибка аппроксимации).

Какие выводы можно сделать из (10) и (11) ? Во-первых, если шаг h слишком велик, то основной «вклад» в ошибку в (10) внесет ошибка аппроксимации e_n. Но даже, если h достаточно мало, ошибка аппроксимации (11) все равно может стать очень большой, если вычислительный процесс (x_n,y_n) приближается к области, в которой функция f или одна из её производных ¶ f/ ¶ x, ¶ f/ ¶ y принимают неограниченные значения. Следует обратить внимание на то, что с уменьшением h область, в которой ошибка аппроксимации e_n будет велика, будет уменьшаться, хотя полностью никогда не исчезнет.

Что же представляют собой области больших e_n с точки зрения дифференциального уравнения? Возможно, Вы уже заметили, что существует непосредственная связь между областями больших ошибок e_n и областями нарушения достаточных условий существования и единственности решения в теореме Коши. Действительно, как следует из теоремы Коши-Пикара [2] неограниченность f(x,y) нарушает достаточное условие существования решения вида y=y(x), а неограниченность ¶ f/ ¶ y — нарушает достаточные условия его единственности. Таким образом, в окрестностях особых решений дифференциальных уравнений, а также в окрестностях изоклин вертикальных наклонов можно ожидать ошибок аппроксимации при численном интегрировании. Возможна также геометрическая интерпретация особенностей таких областей. Действительно, в области больших значений f в методе Эйлера (4) происходит шаг по x на величину h, а y изменяется вдоль прямой большого наклона, равного f (рис.1). В областях с неограниченными значениями частных производных ¶ f/ ¶ x, ¶ f/ ¶ y происходит «бесконечно» быстрая «переориентация» поля направления. При этом ломанная Эйлера как бы не успевает отслеживать изгиб интегральной кривой.

Рассмотрим теперь первое слагаемое в (10). Оно определяет преобразование ранее возникшей ошибки на текущем шаге.
Если ç 1+h( ¶ f/ ¶ y)(x_n,y_n) ç ¶ f/ ¶ y, то ошибка начинает неограниченно усиливаться.

Процесс усиления, неограниченного нарастания (в процессе численного интегрирования) однажды возникшей ошибки называют вычислительной неустойчивостью. Рассмотрим этот процесс отдельно, а именно, как процесс преобразования однажды возникшей ошибки. Пусть на шаге n ошибка равняется d _n. Тогда новое значение y_n+1 будет подсчитано с ошибкой

Считая d _n , h малыми и отбрасывая нелинейные по d _n ,h члены в разложении Тейлора, в линейном приближении получим

Вычислительная неустойчивость будет наблюдаться при таких h, при которых
выражение ç 1+ ( ¶ f/ ¶ y)(x_n,y_n) ç принимает значения, больше единицы. Для систем уравнений аналогичное условие заключается в том, что определитель матрицы (E+(Df/Dy)(x_n,y_n)) (где E — единичная матрица, а (Df/Dy) — матрица Якоби по переменным y) должен стать больше единицы. Например, для системы уравнений

(штрих означает дифференцирование по t), матрица Якоби будет иметь вид:

Её определитель D = 1,01 × h 2 — 0,2 × h + 1 > 0 . Условие ê D ê > 1 будет выполнено при h > 0,2/1,01 » 0,2 . Следовательно, начиная с этого значения, при численном интегрировании по методу Эйлера для данной системы возникнет вычислительная неустойчивость (полезно это проверить в вычислительном эксперименте).

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

Проанализировать в полном объеме этот процесс сложно, поэтому рассмотрим только влияние ошибок конечно-разрядной арифметики при выполнении операции сложения в (4). Предположим, что дифференциальное уравнение dy/dx = f(x,y) имеет частное решение y = a , где величина a ¹ 0 и имеет, например, порядок 1. В этом случае y = a есть изоклина горизонтального наклона. Предположим также что решение задачи Коши для x_o,y(x_o)=y_{o ¹} a стремится с ростом x к значению y=a (рис.3). Тогда в процессе численного интегрирования f(x_n, _n) ® 0 при n ® ¥ .

Предположим также, для определенности, что интегрирование производится с шагом h » 10 -4 , а разрядность вещественной переменной на компьютере составляет семь разрядов (в десятичной системе). Начиная с некоторого шага f(x_n,y_n), как очевидно, станет меньше 10 -3 . При этом произведение h × f(x_n,y_n) в (4) станет меньше 10 -7 . При сложении этого числа
с y_{n »} a » 1 произойдет выравнивание порядков слагаемых и старший значащий разряд произведения h × f(x_n,y_n) выйдет за пределы разрядной сетки и будет утерян. Следовательно, к y_n будет прибавлен чистый нуль! При всех последующих шагах будет сохраняться постоянное (!) значение y_{n ¹} a.

Заметим, что с уменьшением h эта ошибка станет только больше. Например, при h=10 -7 изменения y_n прекратятся начиная уже со значений f(x_n,y_n) » 1 (рис.4).

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

Особенно интересные эффекты, связанные с численным интегрированием, могут наблюдаться в существенно-нелинейных системах, свойства которых значительно изменяются при переходе фазовых траекторий из одних областей фазового пространства в другие. Поясним это, сравнив два примера. В качестве первого рассмотрим линейную систему (12). Её состояние равновесия находится в точке ноль, а квадрат расстояния до него r 2 = x 2 + y 2 .

Вычисляя в силу системы (12) производную от времени, получим

d r 2 / dt= -2 × 0,1 r 2 r 2 на шаге вычислительного процесса, то

r 2 _n+1 — r 2 _n /h = (-2 × 0,1 + (0,1 2 + 1) h) r 2 n (14)

Чтобы величина (14), также как и (13) была отрицательна (только в этом случае при численном счете будет наблюдаться устойчивость состояния равновесия) необходимо, чтобы шаг h был достаточно мал 0,2 h/ 1,01> h (Заметьте, что мы получили несколько другим способом ранее выведенное условие устойчивости процесса численного интегрирования для (12)). Выбирая такой шаг, Вы будете получать качественно правильные результаты при счете. Проведем аналогичное рассмотрение существенно-нелинейной системы

x = y — (x 2 + y 2 )x
y = -x — (x 2 + y 2 )y,

(15) которая отличается от (12) тем, что коэффициент 0,1 заменен на r 2 .

Нетрудно показать, что для неё

d r 2 / d t = -2 r 4 r 2 r 4 × 0,1 r 2 . Это означает, что скорость приближения фазовой траектории к состоянию равновесия в системе (15) убывает существенно быстрее, чем в (12). Это должно приводить к тому, что шаг h , достаточно малый для устойчивости вычислительного процесса при больших r , может стать недостаточно малым при приближении к состоянию равновесия. Действительно, для (15)

Пренебрегая величиной r _n 4 по сравнению с единицей, получим приближенное условие убывания расстояния до состояния равновесия при численном интегрировании: r _n 2 > h/2.Таким образом, при сколь угодно малом h найдется такая окрестность состояния равновесия системы (15), в которой точка (x_n,y_n) будет удаляться от него в процессе интегрирования. Следовательно, при численном интегрировании методом Эйлера у системы (15) будет обнаруживаться предельный цикл, которого в действительности нет. Его радиус будет примерно равен (рис.6).

Важно отметить, что природа описанного эффекта скрыта не в численном методе, а в негрубости самого уравнения (15).

Исследования в программной лаборатории.

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

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

При выполнении вычислительных экспериментов необходимо понять характер проявления различных источников ошибок численного интегрирования для метода Эйлера и уметь объяснить причины возникновения ошибок в конкретных примерах. Рекомендуется также выполнить экспериментальное исследование метода Штермера, Рунге-Кутта, а также метода Эйлера с регулировкой шага.

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

1. dy /dx = (1 — y 2 ) 1/2 dx;
2. dy / dx = (y + (y 2 — 4x 2 ) 1/2 )/x;
3. dy / dx = (4 — y 2 ) 1/2 /y;
4. dy / dx = 0,5(y 2 — 1) ;
5. dy / dx = xy(1 + 2y)/ (1 — y 2 );
6. dy / dx = — x + (x 2 + 2y) 1/2 ;
7. dy / dx = x (4 — y 2 );
8. dy / dx = 1 — y 4 ;
9. x / = 1 — y 2 ; y / = xy(1 + 2y);
10. x / = y — 0,1x ; y / = — x — 0,1y
11. x / = y — x(x 2 + y 2 ); y / = — x — y(x 2 + y 2 )
12. x / = — x + 2y ; y / = -2x + y
13. dy/ dx = — xy(1 — x 2 )

Литература

1. Компьютерная программа “Элементы математического моделирования”.//Рег.ном. ОФАП ВШ России 1992г. 025.7000.290. ЦИФ Гос ФАП 5092000058. Городецкий С.Ю.,
   Котва Т.А., Скорнякова Б.Л., Хентов А.А., Шпилькерман Б.М.
2. Степанов В.В. Курс дифференциальных уравнений. М.:Государственное издательство технико-теоретической литературы.1950.