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

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

4.2. Вычислительные методы решения дифференциальных уравнений
Численные методы интегрирования обыкновенных дифференциальных уравнений и их систем вида y ‘ = f(x,y) можно построить только для случая известных начальных значений всех интегрируемых переменных (для задачи Коши). Общее решение дифференциальных уравнений содержит произвольные постоянные, которые недопустимы в математических моделях, поэтому в качестве искомой функции используется определенное начальными условиями частное решение.

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

А) Методы Эйлера. Простейший метод Эйлера основан на аппроксимации производной простейшей разностной схемой вида:

.

Отсюда в силу решаемого дифференциального уравнения y ‘ = f(x,y) выводится разностное уравнение метода:

Рис. 28 дает геометрическое представление об этом методе: в силу вида исходного дифференциального уравнения функция f(x_k,y_k) представляет собой значение производной в левом конце интервала x, называемого шагом интегрирования. Тогда разностное уравнение метода Эйлера просто описывает выходящую из левого конца шага интегрирования касательную к неизвестной искомой интегральной кривой (изображенной пунктиром).

Р
ис. 28.

Очевидно, что небольшую неизбежную погрешность при такой аппроксимации можно обеспечить только малым шагом интегрирования x. Поэтому численное решение задачи Коши на достаточно большом промежутке изменения аргумента – очень кропотливая процедура, немыслимая без вычислительной техники.

Простейший метод Эйлера относится к методам I порядка, поскольку использует в разностной формуле значение функции в одной точке.

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

Простейший метод Эйлера на практике почти не используется. Наибольшее распространение получили модифицированные методы Эйлера II порядка. Идея первой модификации заключается в выполнении шага интегрирования за два полушага и дает уравнение:

.

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

.

Общими недостатками методов Эйлера I порядка являются невысокая точность и слабая устойчивость (погрешность одного шага интегрирования не только не компенсируется на последующих шагах, а растет – см. § 4.4). Например, для интегрирования уравнений динамики полета транспортных самолетов в условиях, близких к установившимся, с практической точки зрения допустимо пользоваться простейшим методом Эйлера. Но при исследовании неустановившихся режимов полета этого недостаточно – следует применять модифицированные методы Эйлера, а для моделирования движения самонаводящихся ракет использование методов Эйлера практически недопустимо.

Б) Методы Адамса используют значения функции в нескольких предыдущих точках (учитывают предысторию поведения функции: y_k–1. ) для исправления направления касательной. Формула метода Адамса I порядка совпадает с формулой простейшего метода Эйлера:

а формулы более высокого порядка строятся наращиванием формул меньшего порядка:

– II порядка,

– III порядка.

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

В) Методы «прогноз-коррекция» осуществляют расчет в два шага: предварительный расчет – «прогноз» («предсказание») и последующее уточнение – «коррекцию» . Для построения формул метода «прогноз-коррекция» определенного порядка используются формулы метода Адамса того же порядка, например, для простейшего метода I порядка:

– «предсказание»,

– «коррекция».

Геометрическая интерпретация этого метода I порядка показана на рис. 29.

Рис. 29.

При предсказании, как и при методе Эйлера (рис. 28), решение разностного уравнения на шаге отклоняется от точного решения в сторону выпуклости функции, так как строится с помощью касательной в начале шага. В свою очередь коррекция приводит к отклонению в сторону вогнутости, так как строится с помощью прямой, проведенной из той же исходной точки шага, но с наклоном, соответствующим наклону касательной в конце шага. Таким образом, разность между и может служить мерой погрешности численного интегрирования на одном шаге. Т.е. методы «прогноз-коррекция» выгодно отличаются от ранее описанных методов тем, что допускают контроль величины погрешности на каждом шаге интегрирования. Это можно использовать для повышения точности расчетов с помощью уменьшения шага или для экономии времени расчетов с помощью увеличения шага x.

Для всех разностных методов справедливо утверждение: чем меньше x, тем меньше погрешность на шаге, тем выше точность интегрирования дифференциальных уравнений. Однако нельзя заранее сказать, какова должна быть величина x для обеспечения заданной точности. Поэтому расчеты с неприемлемой погрешностью просто идут «в корзину». В отношении этого выгодно отличаются разностные методы, которые позволяют не только контролировать погрешность, но и изменять шаг в процессе интегрирования. Этим последним удобством обладают все разностные методы I порядка, но из них только метод «прогноз-коррекция» дает возможность проконтролировать погрешность и подсказать, когда возникает необходимость изменения шага. Из методов более высокого порядка предоставляют возможность изменения шага интегрирования методы Эйлера и Рунге-Кутта.

Г) Методы Рунге-Кутта m-го порядка используют m внутренних точек шага интегрирования x: , которые задаются характерным для определенной модификации этого метода способом и в которых последовательно вычисляются m значений функции:

а затем производится непосредственно сам шаг интегрирования:

.

Простейший метод Рунге-Кутта I порядка (m = 1) – это метод Эйлера. Наиболее распространенный в программном обеспечении алгоритмических языков – «стандартный» метод Рунге-Кутта IV порядка использует 4 значения функции, вычисленные для двух промежуточных точек на шаге (в середине) и обеих крайних, и соответствующий набор коэффициентов _i:

Наиболее экономичным из методов Рунге-Кутта является метод II порядка следующего вида:

который по форме совпадает со вторым из приведенных выше модифицированных методов Эйлера. (Этот метод разработан как улучшение метода «прогноз-коррекция» I порядка, когда в качестве окончательного значения функции на шаге принимается среднее арифметическое между прогнозом и коррекцией – см. рис. 29)

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

ПРИМЕР. Сравнение методов численного интегрирования дифференциальных уравнений проведем на примере решения следующей задачи Коши:

требуется определить y(2). Результаты численного интегрирования рассмотренными выше методами с шагом x = 0,2 сведены в табл. 2. В ней сравниваются: простейший метод Эйлера, простейший метод «прогноз-коррекция» I порядка, метод Адамса II порядка с началом (первый шаг) по методу Эйлера и метод Рунге-Кутта II порядка. Для краткости в табл. 2 обозначено:

В правом крайнем столбце для сравнения приведено точное решение этой задачи Коши:

Из сравнения результатов численного интегрирования видно, что метод «прогноз-коррекция» действительно дает систематическую погрешность в сторону вогнутости графика функции, в то время как метод Эйлера – в сторону выпуклости (см. рис. 30). Хорошо виден процесс накопления погрешности. Методы I порядка, очевидно, проигрывают перед методами II порядка в точности. При подробном анализе этому можно найти объяснение в накоплении погрешности у первых и в ее явной компенсации у последних (в чем и проявляется устойчивость рассматриваемых методов II порядка). Наименее трудоемкими оказались методы Эйлера и Адамса. Метод Адамса проигрывает в точности методу Рунге-Кутта того же порядка, в основном из-за «нестандартного» и более грубого начала.

Сравнительная таблица методов численного интегрирования

обыкновенных дифференциальных уравнений

k		Эйлер		«прогноз-коррекция»				Адамс		Рунге-Кутта				точное реше-ние
											k2		k1
0	1	0	-1	0	-1	0	-1	0	-1	0
1	1,2	-0,2	-0,8(3)	-0,2	-0,8(3)	-0,1(6)	-0,8611	-0,2[Э]	-0,8(3)	-0,2	-0,8(3)	-0,18(3)	-0,8472	-0,18(3)
2	1,4	-0,3(6)	-0,7381	-0,3389	-0,7579	-0,3183	-0,7727	-0,35	-0,75	-0,3528	-0,7480	-0,3429	-0,7551	-0,34286
3	1,6	-0,5143	-0,6786	-0,4728	-0,7045	-0,4592	-0,7130	-0,4917	-0,6927	-0,4939	-0,6913	-0,4876	-0,6953	-0,4875
4	1,8	-0,6500	-0,6389	-0,6018	-0,6657	-0,5923	-0,6709	-0,6245	-0,6531	-0,6267	-0,6519	-0,6223	-0,6543	-0,6(2)
5	2,0	-0,7778	-0,7265	-0,6368	-0,7197	-0,7512	-0,7532	-0,6234	-0,7501	-0,75

Рис. 30.

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

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

Забытая мысль всегда кажется значительной. Эмиль Кроткий
ещё >>

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

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

Рассмотрим несколько численных методов решения дифференциальных уравнений первого порядка. Описание численных методов приводится для уравнения в виде у’=f(x,y).

Рассмотрим два варианта вывода расчетных формул

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Заключение

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

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