Решение дифференциальных уравнений методами Эйлера и Милна

 

 

 

 

 

Достигли нужной точки, процесс  вычислений закончен.

 

 

Уточненный  метод Эйлера (4.6)–(4.7)

 

 

Начинаем процесс вычислений

 

, 

 

 

Осуществляем «разгон», вычисляя значения и по формулам (4.7):

 

 

Далее вычисляем  по формуле (4.6):

 

 

 

Достигли нужной точки, процесс  вычислений закончен.

Занесем результаты вычислений в таблицу:

Напомним, что 

Методы								Фактическая ошибка
М.Э.(3.3)		0	1
Неявный М.Э.(4.1)		0	1
М. трапеций (4.2)		0	1
М. Хойна (4.3)		0	1
Уточнен. М.Э.(4.6)-(4.7)		0	1

 

Последний столбец в этой таблице со всей очевидностью показывает бо́льшую точность методов второго порядка (метод трапеций, метод Хойна, уточненный метод Эйлера).

5. Исправленный метод Эйлера

Пусть найдено приближенное значение решения задачи (2.1)–(2.2) и требуется вычислить , где . Запишем разложение решения по формуле Тейлора p-го порядка, принимая за базовую точку (т.е. по степеням ) и положим в этом разложении . Имеем

    (5.1)

Если ограничиться двумя  слагаемыми в правой части разложения (5.1), то, как уже было показано выше, получим обычный метод Эйлера (3.3). Посмотрим что будет, если учесть третье слагаемое.

 

При из (5.1) следует равенство

(5.2)

Значение первой производной  в точке , в силу связи (2.1), приближенно известно:

(5.3)

Дифференцируя (2.1), по формуле полной производной

 

находим приближенное значение второй производной:

(5.4)

Подставляя приближенные выражения , и в равенство (5.2), получаем следующую формулу для вычисления при

:

(5.5)

Определяемый ею метод  называется исправленным методом Эйлера.

 

Так как при формулы (5.3) и (5.4) точны, а , согласно начальному условию (2.2), то на первом шаге вычислений по формуле (5.5) будет совершаться ошибка, связанная только с усечением ряда Тейлора. Следовательно, локальная ошибка или, иначе, шаговая погрешность метода (5.5) составляет величину , а это означает, что исправленный метод Эйлера относится к методам второго порядка [2].

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6. Пошаговый контроль точности

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

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

, другой раз  совершается за два шага длины  («транзитом» через точку со значением

 ) и дает значение

 

Поправка Ричардсона в таком случае будет составлять величину

(6.1)

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

При условии  стоит попытаться двигаться дальше с более крупным шагом (например, удвоить ).

 

Пример 2. (продолжение примера 1., см. главу 4)

Посмотрим, что дает применение принципа Рунге к нескольким простым  методам численного решения того же уравнения  с начальным условием . Из точки перейдем в точку за один шаг используя явный и неявный методы Эйлера, метод трапеций и метод Хойна. С помощью полученных значений и найденных ранее теми же методами в примере 1 значений

 подсчитаем  поправки Ричардсона

 

при для явного и неявного методов Эйлера и для методов трапеций и Хойна. Найдем уточненные приближенные значения решения путем прибавления к значениям поправок Ричардсона . Найдем их истинные погрешности

 

Решение:

Метод Эйлера (3.3)

 

 

 

 

 

Начинаем процесс вычислений

 

 

 

 

 

Достигли нужной точки, процесс  вычислений закончен.

 

 

 

 

Неявный метод Эйлера (4.1)

 

 

 

 

 

Начинаем процесс вычислений

 

 

 

 

 

 

Достигли нужной точки, процесс  вычислений закончен.

 

 

 

 

Метод трапеций (4.2)

 

 

 

 

 

Начинаем процесс вычислений

 

 

 

 

 

 

Достигли нужной точки, процесс  вычислений закончен.

 

 

 

 

Метод Хойна (4.3)

 

 

 

 

 

 

 

Начинаем процесс вычислений

 

 

 

 

Достигли нужной точки, процесс  вычислений закончен.

 

 

 

 

Занесем результаты вычислений в таблицу:

Методы
Явный М.Э.
Неявный М.Э.
М. трапеций
М. Хойна

 

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

 

 

 

 

 

7. Методы прогноза и коррекции. Метод Милна

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

 

Для случаев, когда требуется  достичь вполне определенной точности , применяют так называемые методы прогноза и коррекции (схемы предиктор-корректор, методы предсказания и уточнения), заключающиеся в совместном применении явных и неявных методов одинакового или смежных порядков. По явной формуле значение решения задачи (2.1) – (2.2) в текущей (расчетной) точке прогнозируется, т.е. находится его, быть может, достаточно грубое приближение, а с помощью неявной формулы, в правую часть которой подставляется спрогнозированное значение, оно уточняется (корректируется) [2]. Примерами методов прогноза и коррекции могут служить ранее рассмотренные метод Хойна, метод Эйлера-Коши с итерационной обработкой, уточненный метод Эйлера.

 

Далее рассмотрим еще один широко известный метод прогноза и коррекции – метод Милна. Сначала дадим описание этого метода, затем рассмотрим вывод необходимых формул, применяемых в методе.

 

Пусть дано дифференциальное уравнение (2.1) с начальным условием (2.2). Выбрав шаг положим, как обычно, ,  ,

 

Обозначим для краткости

Первые четыре значения искомого решения  (т.н. «начальный отрезок») находим, используя начальное условие (2.2) и применяя какой-либо численный метод решения ОДУ (например, метод Эйлера). Тем самым будут известны , .

 

Дальнейшие значения последовательно определяются по следующей схеме: предполагая, что известны,

1. вычисляем первое приближение для ближайшего следующего значения по формуле

 

2. значение  подставляем в дифференциальное уравнение (2.1) и определяем соответствующее значение
3. находим второе приближение по формуле

 

Милн показал, что абсолютная погрешность значения приближенно равна

 

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

 

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

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

 

Далее приступим к выводу формул метода Милна (7.1) – (7.3). Для этого запишем общий вид первой интерполяционной формулы Ньютона:

 

 

где , , а выражения вида – конечные разности

 

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

 

или

, (7.4)

где 

.   (7.5)

Полагая в формуле (7.4) и почленно интегрируя полученную формулу по в пределах от до , будем иметь

 

Отсюда, учитывая, что  и , находим

 

Так как

,

,

,

 

то, подставляя эти значения в формулу (7.6), после обычных упрощений получим первую формулу Милна:

 

 

Для вывода второй формулы  Милна положим  в формуле (7.4) и проинтегрируем обе части получившегося выражения по в пределах от до , будем иметь

 

Отсюда, выполняя квадратуры, получим

 

Подставляя в формулу  (7.8) известные выражения

,

,

придем ко второй формуле Милна:

 

 

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

 

Отсюда, считая, что четвертая  разность постоянна на интервале длины , получим .

Так как, очевидно, и , то

 

Следовательно, имеем контрольную формулу Милна:

 

Заметим, что если шаг  достаточно мал, то приближенно можно положить

 

Пусть – интересующий нас отрезок, на котором строится решение дифференциального уравнения (2.1), и Тогда из формул (7.10) следует, что предельная абсолютная погрешность на отрезке приближенного решения выражается следующим образом:

 

где при .

Таким образом, глобальная ошибка метода Милна есть величина порядка [4].

 

Пример 3. (продолжение примера 1., см. главу 4)

Дано уравнение  с начальным условием . Найдем методом Милна приближенное значение решения в точке с точностью до

 

Решение:

Метод Милна имеет глобальную ошибку , это значит, что взяв , получим погрешность результата порядка , таким образом, заданная точность практически достигается.

(из начального  условия)

Значения  найдем явным методом Эйлера.

 

 

 

Найдем значения , и

 

 

 

Далее используем метод Милна.

 

 

 

 

 

Проверка:

 

 

 

 

Будем заносить результаты расчетов в таблицу

0		-	-
1		-	-
2		-	-
3		-	-
4