Простейшими численными методами для решения задачи Коши для ОДУ являются следующие.
Метод Эйлера. Этот метод основан на разложении искомой функции Y(x) в ряд Тейлора в окрестностях узлов системы x = xi (i = 0, 1, 2, …, n), в котором отбрасываются все члены, содержащие производные второго и более высоких порядков. Как правило, используется равномерная сетка Dx = xi+1 – xi= h = const (i = ). Разложение запишем в виде
. (8.13)
Заменив значение функции Y(x) в узлах сетки xiзначениями сеточной функции и использовав уравнение (8.4), получим
.
Тогда из (8.13) получим
. (8.14)
При i = 0 для узла x = x1: .
Далее по алгоритму (8.14)
;
. . .
.
Геометрическая интерпретация метода представлена на рис. 8.2.
Рис. 8.2
На рис. 8.2 линия 0 – точное решение, линии 1 и 2 – приближенные решения.
Искомая интегральная кривая y(x), проходящая через точку (x0, y0), заменяется ломаной с вершинами в точках (xi, yi). Каждое звено ломаной имеет направление, совпадающее с направлением интегральной кривой (8.4), которая проходит через точку (xi, yi).
Блок-схема алгоритма представлена на рис. 8.3.
Рис. 8.3
Вывод полученных результатов выполняется на каждом шаге, но если необходимо сохранить результаты, то следует ввести массив значений y0, y1, ..., yn.
Локальная погрешность метода Эйлера, как видно из (8.13), оценивается как О(h2). Весь интервал [a, b] разбивается на n частей, тогда общая погрешность
n×О(h2) = О(h2) = О(h) – 1-й порядок.
Для оценки погрешности при машинном расчете пользуются двойным просчетом, т. е. на отрезке [xi, xi+1] расчет повторяют с шагом h/2, и погрешность более точного решения (при шаге hi/2) оценивается как разность .
Метод Эйлера с пересчетом. При данном подходе рекуррентное соотношение (8.14) видоизменяется, а именно, вместо f(xi, yi) берут среднее арифметическое между f(xi, yi) и f(xi+1, yi+1).
Тогда
(8.15)
Это неявная схема. Она реализуется в две итерации: сначала находится первое приближение по (8.14), считая yi начальной:
, (8.16)
затем (8.16) подставляется в правую часть (8.15) вместо yi+1:
(8.17)
Геометрическая интерпретация метода представлена на рис. 8.4.
Рис. 8.4
С помощью метода Эйлера с пересчетом можно производить контроль точности, сравнивая yi+1 и ỹi+1.
Если величина | ỹi+1 – yi+1 | сравнима с заданной точностью e, то шаг можно увеличивать, если больше, то уменьшать, т. е. имеет место схема двойного просчета с оценкой погрешности по величине:
» ,
где yi,– приближенные значения, полученные с шагом h и h/2 соотвественно; y(xi) – точное решение в точке х = xi.
Метод Эйлера с последующей итерационной обработкой. Метод Эйлера можно еще более уточнить, применяя итерационную обработку каждого полученного значения yi. А именно, сначала исходя из первого грубого приближения по (8.16)
,
строят итерационный процесс согласно (8.15) по следующей схеме:
, k = 1, 2, … . (8.18)
Итерации продолжают до тех пор, пока в двух последовательных приближениях , не совпадут соответствующие десятичные знаки, и полагают yi+1 » . Как правило, при достаточно малом шаге h итерации сходятся быстро. Если после трех-четырех итераций не произошло совпадение нужного числа десятичных знаков, то шаг расчетов h уменьшается. После такой обработки значения yi переходят к следующему узлу xi+1.
Пример 8.1. По методу Эйлера составить таблицу решения на отрезке [0; 1] для уравнения с начальным условием y(0) = 1, выбрав шаг h = 0,2.
Результаты вычислений поместим в таблицу, которая заполняется следующим образом:
i
xi
yi
Dyi
Точное у =
1,0000
0,2000
1,0000
0,2
1,2000
0,1733
1,1832
0,4
1,3733
0,1561
1,3416
0,6
1,5294
0,1492
1,4832
0,8
1,6786
0,1451
1,6124
1,0
1,8237
1,7320
В первой строке при i = 0 записывается x0 = 0, y0 = 1,000 и по ним вычисляется f(x0, y0) = 1, а затем Dy0 = hf(x0, y0) = 0,2. Тогда по формуле (8.14) получаем y1 = 1 + 0,2 = 1,2.
Значения x1 = 0,2 и y1 = 1,2000 записываются во второй строке при i = 1. Используя их, можно вычислить
При i = 2, 3, 4, 5 вычисления ведутся аналогично. В последнем столбце таблицы для сравнения помещены значения точного решения.
Из таблицы видно, что абсолютная погрешность для y5, что составляет 5 %.
Замечание. Метод Эйлера легко распространяется на системы дифференциальных уравнений и на ДУ высших порядков при их предварительном приведении к системам ДУ первого порядка.
Рассмотрим систему двух уравнений первого порядка:
(8.19)
с начальными условиями y(x0) = y0 и z(x0) = z0.
Тогда приближенные значения y(xi) » yiи z(xi) » zi вычисляются по формулам
(8.20)
Пример 8.2. Применив метод Эйлера, составить на отрезке [1; 1,5] таблицу значений решения уравнения
(8.21)
при начальных условиях y(1) = 0,77 и y'(1) = –0,44, выбрав шаг h = 0,1.
Решение. Заменим уравнение (8.21) системой уравнений первого порядка, посредством подстановки y' = z, y" = z':
при начальных условиях y(1) = 0,77 и z(1) = –0,44. Таким образом, получим
Результаты вычисления по формулам (8.20) записаны в таблице:
i
xi
yi
DYi
f1i= zi
Dzi
1,0
0,77
–0,044
–0,44
–0,033
–0,33
1,1
0,726
–0,047
–0,473
–0,030
–0,296
1,2
0,679
–0,050
–0,503
–0,026
–0,260
1,3
0,629
–0,053
–0,529
–0,022
–0,222
1,4
0,576
–0,055
–0,551
1,5
0,521
Таблица заполняется следующим образом. Записываем в первой строке i = 0, x0 = 1,0; y0 = 0,77; z0 = –0,44.
Далее находим
Используя формулы (8.20), получаем
Таким образом, во второй строке таблицы можно записать i = 1; x1 = 1,1; y1 = 0,726; z1 = –0,473. По этим значениям находим
И, следовательно,
Заполнение таблицы при i = 2, 3, 4, 5 производятся аналогично.
Метод Рунге - Кутта. На его основе могут быть построены разностные схемы разного порядка точности. Идея реализации метода состоит в подгонке ряда Тейлора при разложении искомой функции y = y(x) в окрестностях узлов сетки в плане повышения точности этого разложения, а именно, увеличение числа производных высшего порядка без их непосредственного определения из-за сложности аналитических выражений полных производных по x от функции f(x, y).
Рассмотрим наиболее широко применяемую на практике разностную схему четвертого порядка.
Ее алгоритм состоит в следующем:
(8.22)
где ; ;
; .
В данной расчетной схеме Рунге - Кутта на каждом шаге вычисления yi нужно четыре раза обратиться к правой части уравнения f(x, y), т. е. метод Рунге - Кутта (8.22) требует бóльшего объема вычислений, однако это окупается повышенной точностью, что позволяет проводить расчет с большим шагом.
Можно показать, что метод Эйлера и его модифицированный вариант является аналогом метода Рунге - Кутта первого и второго порядков, однако для достижения у них одинаковой точности шаг расчета будет значительно меньше.
Для данного метода шаг расчета можно изменять при переходе от одной точке к другой. Для контроля правильности выбора шага h рекомендуется вычислять дробь
.
Величина Q не должна превышать нескольких сотых. В противном случае h следует уменьшать.
Оценка погрешности метода затруднительна. Чаще всего используется грубая оценка погрешности по формуле , где y(xn) – значение точного решения уравнения (8.4) в точке xn, а y, yn – приближенное решение, полученное с шагом h/2 и h.
При реализации (на ПК) метода Рунге - Кутта с автоматическим выбором шага обычно в каждой точке xi делают двойной просчет сначала с шагом h, потом с h/2. Если полученное yi при этом различается в пределах допустимой точности, то шаг h для следующей точки xi+1 удваивают, в противном случае берут половинный шаг.
В заключение следует отметить, что одношаговые методы Рунге - Кутта успешно могут быть применены к решению систем ДУ первого порядка.
). Разложение запишем в виде
. (8.13)
.
. (8.14)
.
;
.
О(h2) = О(h) – 1-й порядок.
(при шаге hi/2) оценивается как разность 
.
(8.15)
, (8.16)
(8.17)
» 
,
– приближенные значения, полученные с шагом h и h/2 соотвественно; y(xi) – точное решение в точке х = xi.
,
, k = 1, 2, … . (8.18)
, 
не совпадут соответствующие десятичные знаки, и полагают yi+1 » 
с начальным условием y(0) = 1, выбрав шаг h = 0,2.
, что составляет 5 %.
(8.19)
(8.20)
(8.21)
(8.22)
; 
;
; 
.
.
, где y(xn) – значение точного решения уравнения (8.4) в точке xn, а y
, yn – приближенное решение, полученное с шагом h/2 и h.