8. Лабораторная работа №8. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

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

8.1.1. Постановка задачи

Рассмотрим первого порядка.

Задача Коши:

ìy'(x) = f (x, y),

í

îy(x0 ) = y0

обыкновенные

дифференциальные

уравнения

x Î[x0 , b]

.

(8.1)

Пусть требуется найти решение y(x) на отрезке [a,b], где x0 = a .

Применим

к

отрезку[a,b]

равномерное

разбиение, т. е. получим

h =

b - a

и xk = x0 + kh , где xn

= b , xk

– узлы сетки, h – шаг сетки.

n

Обозначим через y(xk )

точное

значение функции y(x) в

точке

xk , через

yk

приближенное

вычисленное значение

функцииy(x) в

точке xk .

8.1.2. Метод Эйлера

Разложим

в

ряд

Тейлора

в

x

значение

функции

точкеk

y(xk + h) = y(xk +1 ) .

h2

y(x

+ ) = y(x

k

)

+ hy' (x

) +

y"(x)

, где

x

k

< x < x

k +1 .

k 1

k

2

Согласно

задаче

Коши(8.1) y'(xk ) = f (xk , y(xk )) ,

тогда

h2

разложение Тейлора y(xk +1 ) = y(xk ) + hf (xk

, y(xk )) +

y" (x ) .

2

Полагая, что значение функции в следующем узле получается,

таким образом,

99

yk +1

yk +1 = yk + hf (xk , yk ) .

(8.2)

Эта формула и определяет метод Эйлера.

Замечание. Метод Эйлера имеет первый порядок точности О(h).

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

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

Проинтегрируем		дифференциальное уравнение(8.1) на отрезке
[xk ; xk +1 ], получим
xk +1	xk +1	xk +1
ò y' (x)dx = ò f (x, y)dx , y(xk +1 ) - y(xk ) = ò f (x, y)dx .
xk	xk	xk
Воспользуемся формулой трапеций, тогда получим

y(xk +1 ) = y(xk ) + h (f (xk , y(xk )) + f (xk +1 , y(xk +1 ))).
2
Эта формула дает приближенное значение
yk+1 = yk + h (f (xk , yk ) + f (xk +1, yk +1 )).	(8.3)
2

Формула (8.3) – это неявная формуламетода Рунге-Кутта II

порядка.

Воспользуемся методом предиктор– корректор для избавления от «неявности». Заменим в правой части равенства(8.3) по формуле Эйлера на

yk*+1 = yk + hf (xk , yk ) .				(8.4)
Затем подставим yk*+1 в формулу (8.3)			вместо yk +1	в правой части
yk +1 = yk +	h	(f (xk , yk ) + f (xk +1, yk*+1 )).		(8.5)
2

Формула (8.5) – явная формула Рунге-Кутта II порядка. Формула (8.4) – предиктор, формула (8.5) – корректор. Точность метода О(h2).

100

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

y

k +1

= y

k

+

h

(F + 2F + 2F + F )

.

(8.6)

6

1

2

3 4

F1 = f (xk , yk ) ,

æ

h

h

ö

F

=

f ç x

k

+

, y

k

+

F

÷

,

2

2

2

è

1

ø

æ

h

h

ö

F

=

f ç x

k

+

, y

k

+

F

÷

,

2

2

3

è

2

ø

F4 = f (xk +1 , yk + hF3 ).

Замечание. Точность метода О(h4).

8.1.4. Выбор шага интегрирования

Точность расчетов существенным образом зависит от величины

шага интегрирования h,	поэтому	важно правильно выбрать		его
начальное значение h0.
Выбор начального	шагаh0	проведем на	примере	метода
Рунге-Кутта IV порядка.	Итак, пусть ε – заданная		точность счета.

Поскольку метод Рунге-Кутта имеет точность четвертого порядка относительно шага h, должно выполняться условие h4 = ε. Кроме того,

отрезок [a,b] должен быть разбит на четное число частей. Поэтому начальный шаг h0 должен быть определен из двух условий:

h0 = 4		,	b - a
	e			– четно.				(8.7)
			h0
Наибольшее h0, удовлетворяющее						условиям (8.7),		является
грубым	приближением				начального	.шагаДля	его	уточнения

поступаем следующим образом. Находим решение задачи Коши в точке х0 + 2h0 по формулам Рунге-Кутта с шагамиh0 и 2 h0, получаем

два значения и ~ . Путем увеличения или уменьшения шага в два y2 y2

раза (не обязательно однократного) подберем наибольшее значение

101

f (xk , yk ) .

h0, при котором будет выполнено неравенство	1		~
	15		y2 - y2	< e	. Это и

будет величина шагаh, с которым решается задача Коши методом Рунге-Кутта.

	8.1.5. Многошаговые методы Адамса
Рассмотренные методы			Рунге-Кутта описываются				формулой
yk +1 = Ф(xk , yk ) , т. е. только в одной точке используется только один
шаг.
Идея:	получить	значение,		на	следующем	шаге используя
значение не только в одной точке, но и в точках, стоящих перед ней
(N + 1 шаг)
yk +1 = Ф(xk ,..., xk -N ; yk ,..., yk - N ) .
Рассмотрим интегральное				представление		дифференциального
уравнения (8.1)
	xk +1
y(xk +1 ) = y(xk ) + ò f (x, y)dx .							(8.8)
	xk
Идея:	построить		интерполяционный			полином	Лагранжа
L(xk -N , xk -N +1 ,..., xk -1, xk )		для функции			f (x, y) .
I. Явные методы Адамса
Пусть известны значение yk , yk -1 ,..., yk -m , т. е. известны значения
f (xk , yk ),...,	f (xk -m , yk -m ) . На				данных	значениях	построим

интерполяционный полином Лагранжа степениm. Рассмотрим два случая.

а) Пусть m = 0, т. е. известно

L0 (x) = f (xk , yk ) .

Тогда из (8.8) получим метод Эйлера yk +1 = yk + hf (xk , yk ) .

б) Пусть m = 1, т. е. известны две точки. Получим

102

L ( x) =

x - xk -1

f

k

-

x - xk

f

k -1 .

1

h

h

Тогда из (8.8) следует, что

y

= y

é(x - xk -1 )2

f

(x - xk )2

f

ùxk +1

= y

4h

f

h

f

h

f

+

+ ê

-

- ú

+

-

-

-

,

2h

2h

2

2

2

k 1

k

ë

k

k 1 ûx

k

k

k

k 1

k

yk+1 = yk

+

h

(3 f (xk , yk ) - f (xk -1, yk-1 )).

(8.9)

2

Формула (8.9) определяет явный метод Адамса. Метод Адамса

имеет второй порядок точности. II. Неявные методы Адамса

Построим

полином

Лагранжа

на:

f (xk +1, yk +1 ),..., f (xk -m , yk -m ) .

Пусть m = 0, тогда

L	(x) = -	x - xk +1	f	k	-	x - xk	f	k +1
1		h				h
и из (8.8) получаем неявный метод Адамса

yk +1 = yk

é

(x - xk +1 )2

fk

(x - xk )2

ùxk+1

h

h

,

+ ê-

2h

+

2h

fk +1 ú

= yk

+

2

fk +

2

fk +1

ë

ûxk

yk+1	= yk	+	h	(f (xk , yk ) + f (xk +1, yk +1 )).		(8.10)
		2
Так	как (8.10) неявная				формула (требуется	разрешение
относительно yk +1 ), применяют схемы предиктора – корректора.
То есть предиктор по (8.9)
yk*+1	= yk	+	h	(3 f (xk , yk ) - f (xk -1, yk-1 )),
		2
и корректор по (8.10)
yk+1	= yk	+	h	(f (xk , yk ) + f (xk +1, yk*+1 )).
		2

103