метод последовательных приближений

М етод последовательных приближений (или метод Пикара) является аналитическим, т. е. позволяет получить приближённое решение задачи Коши, определяемой дифференциальным уравнением (1) с начальным условием (2), в виде формулы. Возник метод в связи с доказательством теоремысуществования и единственности решения задачи Коши (гл. 1).

Пусть в условиях данной теоремы требуется найти решение уравнения (1) с начальным условием (2). Проинтегрируем обе части уравнения (1) от х ₀до x:

, откуда

у (х) = у ₀+ . (7)

Очевидно, что решение интегрального уравнения (7) будет удовлетворять уравнению (1) и начальному условию (2). Действительно, при х = х ₀ получим

у (х ₀) = у ₀+ = у ₀.

Применим к интегральному уравнению (7) метод последовательных приближений. Заменим в равенстве (7) неизвестную функцию у данным значением у ₀, получим первое приближение

у ₁(х) = у ₀+ .

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

у ₂(х) = у ₀+

и т. д. В общем случае итерационная формула имеет вид

у _n(х) = у ₀+ (n =1, 2,...). (8)

Применив неоднократно формулу (8), получим последовательность функций

у ₁(х), у ₂(х),..., у _n(х),.... (9)

Можно доказать [1, 2, 3], что эта последовательность сходится и

= у (х),

т.е. предел последовательности является решением интегрального уравнения (7), а следовательно, и дифференциального уравнения (1) с начальным условием (2). Это означает, что k -й член последовательности (9) является приближением к точному решению уравнения (1)
с определённой степенью точности.

Погрешность k -го приближения можно оценить формулой

, (10)

где L — постоянная Липшица; М — верхняя грань модуля функции f, т.е. ;

величина d для определения окрестности вычисляется по формуле , числа а и b — из неравенства Липшица (гл. 1).

Пример 1. Найти три последовательных приближения решения дифференциального уравнения у' = x + y ²,удовлетворяющего начальному условию у (0) = 1.

Решение. В качестве начального приближения возьмём

у ₀= у (0) = 1,

тогда:

первое приближение у ₁(х) = у ₀+ = 1+ ,

второе приближение у ₂(х) = у ₀+ = 1+ ,

третье приближение у ₃(х) = у ₀+ = 1+ .

Вычисления интегралов и построение графиков полученных функций у ₁(х), у ₂(х), у ₃(х) проведём в системе MathCAD. Результаты решения представлены на рис. 14.

Оценим погрешность третьего приближения.

Для определения области G, заданной неравенствами (6), примем а = 1, b = 2. Получим

G: – 1 х 1,–1 y 3.

В прямоугольнике G функция

f (x, y) = x + y ²

определена и непрерывна, причём: ,

= .

По формуле (10) получим

Рис. 14

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

1) записать интеграл и выделить его в рамку;

2) выбрать команду Evaluate (Вычислить) из меню опции Simbolic (Символика) главного меню.

Существует и другой способ вычисления несобственных интегралов в программе MathCAD, по которому следует:

1) записать интеграл и выделить его в рамку;

2) выбрать команду Simplify (Упростить) из меню опции Simbolic (Символика) главного меню.

Пример 2. Найти пять последовательных приближений решения дифференциального уравнения

у' = х – 2 у,

удовлетворяющего начальному условию у (0) = 0.

Сравнить полученные приближения с точным решением.

Решение. В качестве начального приближения возьмём

у ₀= у (0)= 0.

Решение данного уравнения, проведённое в системе MathCAD, показано на рис. 15.

Рис. 15

МетодЭйлера

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

Суть метода заключается в том, что искомую интегральную кривую y = y (x) заменяют ломаной M₀M₁M₂..., звенья которой являются касательными к интегральным кривым (рис. 16).

Рис. 16

Пусть требуется решить задачу Коши, т.е. найти решение дифференциального уравнения (1) с начальным условием (2) в виде функции y = y (x). Выбрав шаг h, построим, начиная
с точки х ₀, систему равноотстоящих точек:

х _i = х ₀ + ih (i = 0, 1, 2, ...).

Вместо искомой интегральной кривой y = y (x) на отрезке [ х ₀, х ₁]рассмотрим отрезок касательной L ₁к ней в точке М₀(х ₀, y ₀). Уравнение касательной L ₁, в силу (1), имеет вид

y = y ₀ + f (х ₀, y ₀) (х – х ₀).

При х = х ₁из уравнения касательной L ₁получим

y ₁ = y ₀ + hf (х ₀, y ₀),

откуда видим, что приращение функции на первом шаге имеет вид

у ₀ = hf (х ₀, y ₀).

Аналогично, проводя касательную L ₂к некоторой интегральной кривой семейства в точке М₁(х ₁, y ₁), получим

y = y ₁ + f (х ₁, y ₁) (x – х ₁).

При х = х ₂имеем y ₂ = y ₁ + hf (х ₁, y ₁), т.е. y ₂получается из y ₁добавлением приращения

у ₁= hf (х ₁, y ₁).

Таким образом, значения искомой функции y (x) могут быть определены по формулам:

y_i ₊₁ =y_i + у _i, у _i = hf (х _i, y_i), (11)

где i = 0,1,2,..., которые называются вычислительными формулами метода Эйлера.

При этом искомую интегральную кривую y = y (x), проходящую через точку М₀(х ₀, y ₀), приближённо заменяем так называемой ломаной Эйлера M₀M₁M₂..., звенья которой M _i M _i ₊₁прямолинейны между прямыми x = x_i, x = x_i ₊₁ и имеют подъём

Метод Эйлера является простейшим численным методом, удобным в применении, однако он имеет ряд существенных недостатков. Основной из них — малая точность. Она равна порядку h ², причём с каждым шагом погрешность возрастает, т.е. происходит систематическое накопление ошибок. Поэтому на практике часто используют способ двойного счёта — с шагом h и с шагом h / 2. Совпадение десятичных знаков в полученных двумя способами результатах даёт естественные основания считать их верными.

Пример.

1. Найти методом Эйлера численное решение дифференциального уравнения
у' = x ³ + y,удовлетворяющее начальному условию у (0) = 1, на отрезке [0, 1] с шагом h = 0,1.

2. Найти точное решение уравнения у' = x ³ + y и сравнить его с приближённым на отрезке [0, 1].

Решение.

1. Для данного уравнения вычислительные формулы (11) имеют вид:

y_i ₊₁ =y_i + у _i, у _i = 0,1(х _i ³ + y_i),

где i = 0,1,2,..., 10.

Учитывая, что погрешность метода имеет порядок h ²= 0,01, достаточно в промежуточных результатах брать три цифры после запятой, а во всех y_i сохранять только две цифры.

Результаты вычислений оформим в виде таблицы.

Таблица

i	х _i	y_i	y_i = hf (х _i, y_i) = 0,1(х _i ³ + y_i)
0	0	1	0,1
1	0,1	1,1	0,110
2	0,2	1,21	0,122
3	0,3	1,33	0,136
4	0,4	1,47	1,634
5	0,5	1,62	0,175
6	0,6	1,79	0,201
7	0,7	1,99	0,233
8	0,8	2,22	0,273
9	0,9	2,49	0,322
10	1	2,82	—

2. Данное уравнение у' = x ³ + y является линейным дифференциальным уравнением первого порядка. Решим его методом Бернулли.

Полагая y = uv, имеем

= 0.

Сгруппируем члены, содержащие u в первой степени, получим

= 0.

Полагаем = 0, откуда . Интегрируя, находим , или (постоянную интегрирования не вводим, так как достаточно найти какое-либо частное решение этого вспомогательного уравнения).

Для нахождения u имеем уравнение

или

Разделим переменные, получим , откуда

Интегрируем по частям три раза:

Таким образом, общее решение данного уравнения

y = uv = ,

или y = .

Используя начальное условие у (0) = 1, получим 1 = ‑ 6 + С, откуда С = 7. Следовательно, искомое частное (точное) решение имеет вид

у = .

Вычислим значения полученного точного решения на отрезке [0, 1] с шагом h = 0,1. Результаты округлим до 0,01 и запишем в таблицу.

Таблица

i	х _i	Приближённые значения y_i	Точные значения y (х _i)
0	0	1	1
1	0,1	1,1	1,11
2	0,2	1,21	1,22
3	0,3	1,33	1,35
4	0,4	1,47	1,5
5	0,5	1,62	1,67
6	0,6	1,79	1,86
7	0,7	1,99	2,08
8	0,8	2,22	2,35
9	0,9	2,49	2,66
10	1	2,82	3,03

Сравнение приближённого (численного) решения данного дифференциального уравнения с точным на промежутке [0, 1] проведём с помощью системы MathCAD.

Результаты сравнения, а также численное решение данного уравнения, проведённое методом Эйлера в системе MathCAD, представлены на рис. 17.

Рис. 17

МодификацииметодаЭйлера

Существуют различные уточнения метода Эйлера, повышающие его точность. Цель модификаций — более точно определить направление перехода из точки (х _i, y_i) в точку (х _i ₊ ₁, y_i ₊ ₁). Так, метод Эйлера-Коши предлагает вычислять значения искомой функции y (x) по формулам:

= y_i + Д у _i, Д у _i = hf (х _i, y_i),

y_i ₊₁ = y_i + h , i = 0,1,2,....

Геометрически это означает, что мы определяем направление интегральной кривой в исходной точке (х _i, y_i) и во вспомогательной точке (х _i ₊ ₁, ). В качестве окончательного берём среднее этих направлений.

Другой модификацией метода Эйлера является усовершенствованный метод ломаных, при котором сначала вычисляют промежуточные значения:

и находят значение направления поля интегральных кривых в средней точке (, ), т.е. = f (, ), а затем полагают

y_i ₊₁ = y_i + h .

Метод Эйлера и его модификации являются простейшими представителями конечно-разностных методов (шаговых методов) для приближённого решения задачи Коши.

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

На рис. 18 и 19 показаны решения дифференциального уравнения у' = x ³ + y,удовлетворяющего начальному условию у (0) = 1, полученные модифицированными методами Эйлера (методом Эйлера-Коши и усовершенствованным методом ломаных) с помощью системы MathCAD.

Рис. 18

Рис. 19

Метод Рунге-Кутта

Рассмотренный выше метод Эйлера относится к семейству методов Рунге-Кутта и является их простейшим частным случаем (методом первого порядка точности). Наиболее известным из методов Рунге-Кутта является классический четырёхэтапный метод четвёртого порядка точности. Его расчётные формулы для решения задачи Коши, определённой уравнениями (1) и (2), имеют вид:

y_i ₊₁ =y_i + у _i; у _i = (k ₁⁽ⁱ⁾ + 2 k ₂⁽ⁱ⁾ + 2 k ₃⁽ⁱ⁾ + k ₄⁽ⁱ⁾), (12)

где k ₁⁽ⁱ⁾= h f (х _i, y_i);

k ₂⁽ⁱ⁾= h f (х _i + , y_i + );

k ₃⁽ⁱ⁾= h f (х _i + , y_i + );

k ₄⁽ⁱ⁾= h f (х _i + h, y_i + k ₃⁽ⁱ⁾), i = 0,1,2,....

Погрешность метода на каждом шаге является величиной порядка h ⁵.

Геометрический смысл использования метода Рунге-Кутта с вычислительными формулами (12) состоит в следующем (рис. 20).

Рис. 20

Из начальной точки М ₀(х ₀, y ₀) сдвигаются в направлении, определяемом углом ₁, для которого tg ₁= f (х ₀, y ₀). Идут в этом направлении на полшага, т.е. до вертикальной прямой
х = х ₀ + . На этом направлении выбирается точка Р ₁с координатами

Затем из точки М ₀(х ₀, y ₀)сдвигаются в направлении, определяемом углом ₂, для которого tg ₂= f (х ₀ + , y ₀ + ), и на этом направлении выбирается точка Р ₂с координатами

Далее из точки М ₀(х ₀, y ₀) сдвигаются в направлении, определяемом углом ₃, для которого tg ₃= f . На этом направлении выбирается точка Р ₃с координатами

(х ₀ + h, y ₀ + k ₃⁽⁰⁾). Этим задаётся ещё одно направление, определяемое углом ₄, для которого tg ₄= = f (х ₀ + h, y ₀ + k ₃⁽⁰⁾). Четыре полученных направления усредняются в соответствии с формулой

= (k ₁⁽⁰⁾ + 2 k ₂⁽⁰⁾ + 2 k ₃⁽⁰⁾ + k ₄⁽⁰⁾).

На этом окончательном направлении и выбирается очередная точка М ₁с координатами (х ₁, y ₁) = (х ₀ + h, y ₀ + ).

Теперь, уже исходя из точки М ₁, все построения с помощью усреднений направлений повторяют сначала. Идут в новом усреднённом направлении до вертикальной прямой х = х ₂, получают точку М ₂(х ₂, y ₂) и т.д.

Эффективная оценка метода Рунге-Кутта затруднительна [2, 4]. Поэтому для определения правильности выбора шага h на практике обычно на каждом этапе из двух шагов применяют двойной пересчёт, а именно: исходя из текущего верного значения y (х _i) вычисляют величину y (х _i + 2 h) двумя способами: один раз с шагом h, другой раз — с двойным шагом 2 h.

Если расхождение полученных значений не превышает допустимой погрешности, то шаг h для данного этапа выбран правильно и полученное с его помощью значение можно принять за y (х _i + 2 h). В противном случае шаг уменьшают в два раза.

На практике при вычислениях по формулам (15) обычно пользуются схемой, приведённой в таблице.

Таблица

i	x	Y	k = hf (х, y)	у
0	х ₀ х ₀ + х ₀ + х ₀ + h	y ₀ y ₀ + y ₀ + y ₀ + k ₃⁽⁰⁾	k ₁⁽⁰⁾ k ₂⁽⁰⁾ k ₃⁽⁰⁾ k ₄⁽⁰⁾	k ₁⁽⁰⁾ 2 k ₂⁽⁰⁾ 2 k ₃⁽⁰⁾ k ₄⁽⁰⁾
—	—	—	—
1	х ₁	y ₁	...	...

Пример. Найти методом Рунге-Кутта решение дифференциального уравнения у' = x ³ + y,удовлетворяющего начальному условию у (0) = 1, на отрезке [0, 1] с шагом h = 0,1.

Решение. Учитывая, что погрешность метода имеет порядок h ⁵= 0,00001, в промежуточных результатах следует брать пять цифр после запятой, а во всех y_i сохранять только четыре цифры. Результаты вычислений оформим в виде таблицы.

Таблица

i	х	y	k = 0,1(х ³ + y)	y
0	0 0,05 0,05 0,1	1 1,05 1,0525 1,1053	0,1 1,10501 1,10526 1,11063	0,1 0,21003 0,21053 0,11063
				0,1052
1	0,1 0,15 0,15 0,2	1,1052 1,1604 1,1634 1,2219	0,11062 0,11637 0,11668 0,11136	0,11062 0,23278 0,21121 0,11136
				0,10556
2	0,2 0,25 0,25 0,3	1,2218 1,2717 1,2752 1,3399	0,12188 0,12874 0,12908 0,13669	0,12188 0,25747 0,25816 0,13669
				0,12903
3	0,3 0,35 0,35 0,4	1,3520 1,4081 1,4124 1,4853	0,13668 0,1451 0,14552 0,15493	0,13668 0,2902 0,29105 0,15493
				0,14548
4	0,4 0,45 0,45 0,5	1,4988 1,5628 1,568 1,6512	0,15493 0,16539 0,16591 0,17762	0,15493 0,33078 0,33182 0,17762
				0,16586
5	0,5 0,55 0,55 0,6	1,6661 1,74 1,7465 1,8425	0,17762 0,19064 0,19132 0,20585	0,17762 0,38128 0,38258 0,20585
				0,19122
6	0,6 0,65 0,65 0,7	1,8588 1,9618 1,9699 2,0826	0,20584 0,22199 0,2228 0,24082	0,20584 0,44399 0,4456 0,24082
				0,22271
7	0,7 0,75 0,75 0,8	2,0833 2,1855 2,1955 2,3268	0,24081 0,26074 0,26173 0,28388	0,24081 0,52148 0,52347 0,28388
				0,26161
8	0,8 0,85 0,85 0,9	2,3468 2,4898 2,5021 2,6585	0,28589 0,3104 0,31162 0,33875	0,28589 0,62079 0,62324 0,33875
				0,31145
9	0,9 0,95 0,95 1	2,6582 2,8545 2,8695 3,0566	0,34129 0,37119 0,37269 0,40566	0,34129 0,74238 0,74537 0,40566
				0,37245
10	1	3,0280

Соответствующее решение данного дифференциального уравнения, полученное методом Рунге-Кутта в системе MathCAD, представлено на рис. 21.

Рис. 21

Лабораторная работа

«Численные методы решения задачи Коши
для обыкновенных дифференциальных уравнений»

Задание 1.

1. Для заданного дифференциального уравнения первого порядка у' = f (x, y) c начальным условием у (a) = c найти приближённое решение в виде многочлена пятой степени.

2. Найти численное решение данного дифференциального уравнения на отрезке [ a, b ] с шагом интегрирования h, округляя результат до 0,001.

3. Найти точное решение заданного дифференциального уравнения у' = f (x, y) и сравнить его с приближённым на отрезке [ a, b ]. Построить графики полученных решений.

Исходные данные для 15-ти вариантов содержатся в таблице.

Вариант	f (x, y)	a	b	с	h
1		0	1	0	0,1
2		0	1	1	0,1
3		0	2	0	0,1
4		p	2p	0	p/10
5		1	2		0,1
6		1	3		0,2
7	4 +	1	2	2	0,1
8		1	2	0	0,2
9		0	2	0	0,1
10		0	2	1	0,2
11		1	2	0	0,2
12	–	0	1	p/4	0,1
13	–			е	0,1
14		0	1	1	0,1
15				0	0,3

Указания к выполнению задания 1

1. Для того, чтобы получить приближённое решение заданного дифференциального уравнения в виде многочлена пятой степени, используйте формулу (3) при k = 0, 1,..., 5.

2. При выборе метода для вычисления точного решения учитывайте то, что дифференциальные уравнения вариантов 1- 4 являются линейными дифференциальными уравнениями, уравнение 5-го варианта — уравнение Бернулли, уравнения 6-8-х вариантов — однородные дифференциальные уравнения, а уравнения 9-15-х го вариантов — дифференциальные уравнения с разделяющимися переменными.

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

Задание 2. Решить задачу Коши для обыкновенного дифференциального уравнения первого порядка у' = f (x, y) на отрезке [ a, b ]при заданном начальном условии у (a) = c и шаге интегрирования h:

1) методом Эйлера с шагом 2 h и с шагом h;

2) модифицированным методом Эйлера (методом Эйлера - Коши или усовершенствованным методом ломаных);

3) методом Рунге-Кутта с шагом 2 h и с шагом h.

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

Вариант	f (x, y)	a	b	с	h
1	2 – sin(x + y)²	2	3	2,3	0,1
2	co s (1,5 x – y ²) – 1,3	– 1	1	0,2	0,2
3	cos(1,5 y + x ²) + 1,4	1	2	0,9	0,1
4	cos(0,6 +y) + 2,5 x	1	3	1,5	0,2
5	1,5 + sin(x + y)	1,5	2,5	0,5	0,1
6		2,6	4,6	1,8	0,2
7	+1	0	2	2,9	0,2
8		0	2	1	0,2
9		0	2	0	0,2
10	arcsin x+x –	0	1	3	0,1
11	4,1 x – y ²+0,6	0,6	2,6	3,4	0,2
12	+ 2 x	0	0,5	0,3	0,05
13	+ 2 y	1,5	2	2,1	0,05
14	+ x + 1	0,1	0,5	1,25	0,05
15	– 0,4	3	5	1,7	0,2