[ Буслов, Яковлев ] Часть 2 - Численные методы

= Ay(n+1) (èëè y(n+1)

в) Метод скалярных произведений

	Этот метод является обобщением метода итераций. Пусть					x	è y произвольные начальные векторы.
	Определим итерации
	ym = Aym¡1 = Amy = Xs						¸mPsy ;
	xm = Axm¡1 = Xs				¸mPsx :
	Аналогично методу итераций убеждаемся, что
		hym; xmi	!	¸		max	:
	hym; xm¡1i

2.2.3 Обратные итерации

Поиск минимального по модулю собственного значения

Пусть y некоторый стартовый вектор. Определим обратные итерации как y(n) = A¡1y(n)) , то есть это прямая задача для нахождения максимального собственного значения ¹max матрицы B = A¡1, обратной к исходной матрице. Очевидно, что минимальное по модулю собственное значение матрицы A равно максимальному по модулю собственному числу обратной матрицы.

B = A¡1		; (¹i =	1	) :
¹i	¸i		¸i

Метод обратных итераций со сдвигом

Пусть A невырожденная эрмитова матрица и ¸¤ некоторое пробное число. Рассмотрим матрицу (A¡¸¤I) , ее собственными значениями являются числа (¸i ¡¸¤) , ãäå ¸i собственные значения исходной

матрицы A . У обратной матрицы (A ¡ ¸¤I)¡1				собственные значения это величины 1
					¸i¡¸¤ . Процедура
метода обратных итераций cо сдвигом
			y(n) = (A ¡ ¸¤I)y(n+1) ;
приводит к нахождению maxi j	1		j . Иными словами, мы находим то собственное значение ¸j , которое
	¸i¡	¸¤
является ближайшим к пробному числу ¸¤ . Варьируя пробное ¸¤					и вновь применяя метод обратных
итераций со сдвигом можно, найти все собственные значения матрицы					A .

2.3 Неэрмитовы матрицы

2.3.1 Дополнительные сведения

В случае если алгебраическая и геометрическая кратности собственных чисел оператора A совпадают, то унитарным преобразованием (то есть преобразованием, сохраняющим скалярное произведение: hUx; Uyi = hx; yi ) (â RN ортогональным преобразованием) оператор приводится к диагональному виду, и на диагонали стоят собственные числа A с учетом кратности. Однако нередка ситуация, когда алгебраическая кратность собственного значения превышает геометрическую (обратное, кстати, невозможно).

21

при определенном выборе базиса (называемым жордановым или каноническим базисом оператора

 CN

A ) матрица оператора становится блочно-диагональной. В каждом из блоков (жордановых клеток) матрица

оператора является верхнетреугольной и имеет вид				:: :: ::		0	0 1
0	0 ¸		1
	¸	1	0			0	0
B	0	0	¸	: : :		0	0 C
B				.	..	.	.	C
B . . .								C
B								C
B								C	:	(1)
B	0	0	0	: : :		¸		C
B							1 C
B								C
B								C
B								C
@	0	0	0	: : :		0	¸	A

Размеры жордановых клеток, их количество, также как и числа ¸ (корни характеристического уравнения) являются инвариантами оператора A (то есть не зависят от выбора жорданова базиса).

 RN жорданов базис приводит к клеткам вида (1) если ¸ вещественный корень характеристическо-

го уравнения матрицы оператора A в каком либо базисе. Поскольку коэффициенты характеристического полинома матрицы оператора в RN вещественны, то вместе с каждым комплексным корнем ¸ = ¹ + iº

он обладает и комплексно сопряженным			¯
			¸ = ¹ ¡ iº . Жорданова клетка в этом случае имеет вид
0		º	1		0	0		0	: : :		0		0	1
	¹º ¹		0		1	0		0	: : :		0		0
B	¡	0	¹		º	1		0	: : :		0		0 C
	0
B														C
B	0	0		º ¹		0		1	: : :		0			C
B			¡										0 C
B														C
B														C
B	0	0	0		0	¹		º : : :			0		0 C		:
B														C
B	0	0	0		0		º ¹ : : :				0			C
B													0 C
B						¡								C
B									.	..	.		.	C
B . . . . . .														C
B														C
B	0	0	0		0	0		0	: : :		0			C
B													1 C
B														C
B	0	0	0		0	0		0	: : :		¹		º	C
B														C
B														C
B	0	0	0		0	0		0	: : :					C
B												º ¹ C
@											¡			A

2.3.2 Метод итераций для максимального по модулю собственного числа кратности 2 в случае жоpдановой аномалии

Остановимся подробно на случае, когда максимальному по модулю собственному значению ¸ оператора A

соответствует жорданова клетка размера		2 £2 . В каноническом базисе u1; u2; : : : ; uN матрица оператора
A имеет вид	0		1	j	0 : : : 0			1
		¸
		0	¸		0	: : :	0
	B	¡	¡	j	¡	¡		C	:
	B			¡			¡ C
	B	0	0					C
	B							C
	B	. .		j		B		C
	B							C
	B			j				C
	B							C
	B	0	0					C
	B							C
	@			j				A

Здесь B матрица, отвечающая оставшимся собственным значениям, конкретный вид которой нас не

интересует. Обозначим дуальный базис через	v1; v2; : : : ; vN : Тогда
Au1 = ¸u1	A¤v1 = ¸v1
	:

Au2 = ¸u2 + u1 A¤v2 = ¸v2 + v1

22

Вектор u1 является собственным для оператора A , соответствующим собственному значению ¸ . Вектор

u2 называется присоединенным. Для сопряженного оператора			A¤ собственным и присоединенным век-
торами, соответствующими собственному значению ¯	N	просто ¸ ) являются векторы дуального
¸ (â R
базиса v1 è v2 , соответственно. Заметим, что u1 = v2 , è			u2 = v1 , то есть собственный вектор для

оператора является присоединенным для сопряженного и наоборот.

Непригодность обычного метода итераций

Будем считать, что собственные значения пронумерованы в порядке убывания модуля и что ¸1 = ¸ . Пусть x произвольный вектор. Разложим его по векторам жорданова базиса и дуального к нему

P

iP

N

N

x = i=1hx; viiui ;

x = =1hui; xivi :

Подействуем на x оператором

A и сопряженным:

iP

P

N

N

Ax = i=1hx; viiAui ; A¤x = =1hui; xiA¤vi ;

èëè

N

Xi

Ax = (¸ < x; v1 > + < x; v2 >)u1 + ¸ < x; v2 > u2

+ < x; vi > Aui ;

=3

N

A¤x = (¸¯ < u1; x > + < u2; x >)v1 + ¸¯ < u2; x > v2

Xi

+ < ui; x > A¤vi :

=3

Аналогично, поскольку

0 0

¸n¡

j

0 : : :

0 1

¸n

n¸n 1

0

: : : 0

An

B

¡

¡

j

¡

¡

C

;

= B

¡

¡ C

B

0

0

C

B

C

B

.

j

Bn

C

B .

C

B

j

C

B

C

B

0

0

C

B

C

òî

@

j

A

Anx = (¸n < x; v1 > +n¸n¡1 < x; v2 >)u1 + ¸n < x; v2 > u2 + : : : :

(2)

Квадратичная форма n-ой степени оператора

A с использованием (2) может быть записана как

N

X

hui; xivii =

hAnx; xi = hAnx;

1

= ¸n n³a + b¸´ + O ([¸0=¸]n)o ;

n

ãäå

a =< x; v1 >< u1; x > + < x; v2 >< u2; x > (= 2 < x; v1 >< u1; x > â Rn),

b =< x; v2 >< u1; x > è

¸0

следующее по модулю за

¸ собственное значение.

В нашей ситуации ¸ вещественно. По аналогии с методом скалярных произведений, применяемым для

эрмитовых матриц, рассмотрим отношение Релея

½n =

h

i

=

h

i

= ¸ (

+ O [¸0=¸]

)

=

Anx; A¤nx

A2nx; x

a +

2nb

h

i

h

i

¸

³

´

An+1x; A¤nx

A2n+1x; x

a +

(2n+1)b

2n

¸

23

= ¸ (1 +

+ O [¸0

=¸]

) = ¸ 1 + O

:

a +

2

nb

n

³

´

µ

¶¾

¸

½

1

a + ¸b

2n

Èòàê, ½n = ¸f1 + O(1=n)g , то есть сходимость при

n ! 1 настолько неудовлетворительная, что теряет

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

Модифицированный метод итераций

Составим квадратное уравнение, для которого ¸ является корнем кратности 2

(t ¡ ¸)2 = t2 + pt + q = 0 p = ¡2¸ ; q = ¸2 :

Коэффициенты

p

è q заранее неизвестны, поскольку неизвестно само

¸ . Попытаемся их определить.

Обозначим xn = Anx и рассмотрим выражение

xn+1 + pxn + qxn¡1 =< x; v1 > f¸1n+1 + p¸1n

+ q¸1n¡1g u1+

2

n

n 1

n 2 1|

{z

}

n

n 1

2

2

n+1

=0

+ < x; v

> f(n + 1)¸1 + pn¸1

¡ + q(n ¡ 1)¸1¡ gu + < x; v > f¸1

+ p¸1 + q¸1¡ g u + : : : =

=< x; v2 > fn¸n¡2 (¸2 + p¸ + q) +¸n ¡ q¸n¡2gu1 + : : : =< x;|

2{z(

1}

+ : : : ;

v1 > ¸n¡

¸2 ¡ q) u

=0

|

{z

}

|

{z

}

.

=0

=0

2

p2

n+1

n

n

1

n+1

поскольку

(¸

¡ q) =

¡ q = 0

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

x

+ px

+ qx

¡

= o(x

)

. При этом, координаты

4

n-ой итерации

xk

ведут себя как соответствующая степень

¸ : xin = (Anx)i » ¸nxi , поэтому естественно

ввести три вектора

yn+1;n;n¡1

=

xn+1;n;n¡1

¸n+1 . Для координат этих векторов, как следует из предыдущего,

выполнено

ykn+1 + pykn + qykn¡1 = O([¸0=¸]n+1) :

Выпишем соответствующие равенства для пары координат, скажем k è l.

yn

yn+1

+ pyn

+ qyn¡1

= O [¸0=¸]n+1

0

yn¡1

l

k

k

k

³

´

»

ln

:

n

n+1

n

n

1

n+1

1

yk

yl

+ pyl

+ qyl ¡

= O ³[¸0

=¸]

´

» 0

yk¡

Домножая первое равенство на yn¡1

, а второе на

yn¡1

и вычитая из первого равенства второе, получаем

l

k

p =

¡

ykn+1yln¡1 ¡ yln+1ykn¡1

+ O [¸0=¸]n+1

=

yknyln¡1 ¡ ylnykn¡1

³

´

=

¡

xkn+1xln¡1 ¡ xln+1xkn¡1

+ O [¸0

=¸]n+1

´

:

xknxln¡1xlnxkn¡1

³

Аналогично, домножая первое равенство на yn

yn

и вычитая их первого равенства второе,

l , а второе на

k

получаем

xkn+1xln ¡ xln+1xkn

q =

+ O [¸0=¸]n+1

:

¡xkn¡1xln ¡ xln¡1xkn

³

´

Заметим, что необходимое количество итераций в предложенном методе, можно контролировать исходя из того, что должно выполняться равенство p2=4 = q .

24

Глава 3

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

3.1 Общие сведения

Уравнение

F (x; u; u0; : : : ; u(n)) = 0

называется обыкновенным дифференциальным уравнением n-го порядка, если F определена и непрерывна в некоторой области G 2 Rn+2 (n ¸ 1) и, во всяком случае, зависит от u(n). Его решением является любая функция u(x), которая этому уравнению удовлетворяет при всех x в определенном конечном или

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

				u(n) = f(x; u; : : : ; u(n¡1)) :	(1)
Решением этого уравнения на интервале I = [a; b] называется функция					u(x) , такая что
1) u(x) 2 Cn[a; b] ;
2) (x; u(x); : : : ; u(n¡1)(x)) 2 D(f) 8x 2 I ,
3) u(n)(x) = f(x; u(x); : : : ; un¡1(x))	8	x	2	I :

3.1.1 Задача Коши

Задачей Коши (начальной задачей) для уравнения (1) называется задача нахождения такого решения уравнения (1), которое удовлетворяет начальным условиям

u(x0) = u0 ; u0(x0) = u00 ; : : : ; u(n¡1)(x0) = u(0n¡1) ;

ãäå u(i)

0 некоторые заданные числа. Справедлива

Теорема Пеано. Åñëè f непрерывна в D, тогда для любой точки x0; u0; : : : ; u(n¡1)

0 , принадлежащей области D, существует решение уравнения (1), определенное в некоторой окрестности точки x0 2 I.

Замечание. Теорема Пеано не гарантирует единственности.

25

существует единственное решение уравнения

Теорема Коши-Пикара. Åñëè f непрерывна в D и удовлетворяет условию Липшица по переменным u, u0, : : :, u(n¡1), òî åñòü

Xn

jf(x; ¹1; ¹2; : : : ; ¹n) ¡ f(x; º1; º2; : : : ; ºn)j < L j¹k ¡ ºkj ;

k=1

то для любой точки (x0; u0; : : : ; u(0n¡1)) 2 D существует единственное решение (1), определенное в некоторой окрестности точки x0 2 I.

Любое уравнение типа (1) можно свести к равносильной ему системе

dudxi = fi(x; u0; u1 : : : ; un¡1) ; i = 0; 1; ; : : : ; n ¡ 1 ;

дифференциальных уравнений первого порядка путем замены высших производных неизвестными функциями (ui(x) = u(i)(x)).

Теорему Коши-Пикара несложно доказать воспользовавшись теоремой о неподвижной точке сжимающего отображения [15]. Действительно, уравнение первого порядка

(u0 = f(x; u) u(x0) = u0

эквивалентно интегральному уравнению

u(x) = u0	+xZ0	f(t; u(t))dt :
	x

По условию f непрерывна и, следовательно, jf(x; u)j · M в некоторой области D0 ½ D, содержащей точку (x0; u0). Выберем ± > 0 так, чтобы:

1) (x; u) 2 D0, åñëè jx ¡ x0j · ± è ju ¡ u0j · ±M;

2) ±L < 1, ãäå L константа, фигурирующая в условии Липшица.

Пусть C0 пространство всех непрерывных функций u, определенных при jx ¡ x0j · ± и таких, что

j	u(x)	¡	u	0j ·	±M с естественной для непрерывных функций метрикой	½(u1; u2) =	max u (x)	¡	u (x)		. Êàê
							x j 1		2	j

замкнутое подпространство полного пространства C[x0¡±;x0+±], пространство C0 является полным. Убедим- ся, что отображение y = Au, определяемое формулой

y(x) = u0	+xZ0	f(t; u(t))dt ;
	x

является сжатием в C0. Действительно, пусть u 2 C0								è jx ¡ x0j · ±, тогда
		y(x)		u0		=	¯ x f(t; u(t)dt¯			±M
	j		¡		j		¯Z	¯	·
							¯x0	¯
							¯	¯
и, следовательно, A переводит C0	в себя. Далее,						¯	¯
							¯	¯
		x
jy1(x) ¡ y2(x)j ·xZ0			jf(t; u1(t) ¡ f(t; u2(t)jdt · L±jju1 ¡ u2jjC0 ;

и поскольку ±L < 1, òî A сжатие и, следовательно, в C0

u = Au. Аналогично доказывается однозначная разрешимость задачи Коши для системы уравнений первого порядка, а, следовательно, и для задачи Коши произвольного порядка.

26

®1u(a) +

3.1.2 Краевая задача

Сформулируем краевую задачу только для уравнений второго порядка, являющуюся одной из самых суще-

ственных. Такая задача имеет вид:		8	u00		= f(x; u; u0); x [a; b];
			®1u(a) + ¯1u0(a) =2°1;
		>
		<	® u(b) + ¯						u0(b) = °2;
где в краевых условиях считается, что		>		2		i		,2		. В отличие от задачи Коши здесь значительно
	j	:i	j	+ ¯			= 0 i = 1; 2
		®
					j		j 6

сложнее исследуется вопрос о существовании решения. Очень важный и наиболее часто встречающийся случай линейное дифференциальное уравнение второго порядка

u00 + p(x)u0 + q(x)u = f(x) ;

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

3.1.3 Задача Штурма-Лиувилля

Задача Штурма-Лиувилля или задача на собственные функции и собственные значения является одновременно и краевой задачей (с однородными краевыми условиями) и обычно записывается в, так называемом,

самосопряженном виде:	¡dx			·k(x)dx¸			+ [q(x) ¡ ¸r(x)] u(x) = 0 ;
		d		du

¯1u0(a) = 0 ; ®2u(b) + ¯2u0(b) = 0 :

Здесь требуется найти те ¸, при которых задача разрешима (собственные значения), и соответствующие им

решения u¸(x) собственные функции, определяемые с точностью до постоянного множителя.

3.1.4 Что понимается под численным решением

Точные (аналитические) методы решения такие методы, когда решение дифференциального уравнения можно получить в виде элементарных функций или квадратур от них, что, естественно, возможно не всегда. Численные методы методы нахождения решений не на всем промежутке изменения независимой переменной, а лишь в дискретном наборе точек x0; x1; : : : ; xN 2 [a; b]. Здесь, правда, следует отметить, что

можно искать решение в виде разложения в ряд по некоторой полной системе функций (скажем, в ряд Фурье) и обрезать его на некотором члене. Однако, вопрос о том, какую систему функций использовать и какое количество членов разложения использовать, является одновременно и численным и аналитическим. Численные методы применимы к очень широкому классу дифференциальных уравнений. В соответствии с двумя типами задач для дифференциальных уравнений, численные методы тоже делятся на два класса: Численные методы решения задачи Коши и численные методы решения краевой задачи и задачи ШтурмаЛиувилля.

3.2 Задача Коши

Рассмотрим задачу Коши для уравнения первого порядка на отрезке	[a; b] :
u0 = f(x; u) ; u(a) = u0 ;	(2)

27

Разобь¼м промежуток [a; b] íà N частей				a = x0 < x1; < : : : ; < xN . Обозначим u(xi) = ui , ãäå u(x)
точное решение задачи Коши, и через yi				значения приближенного решения в точках xi . Существует два
типа численных схем :
1.	явные : yi = F (yi¡k; yi¡k+1; : : : ; yi¡1)			(à);
2.	неявные : yi = F (yi¡k; yi¡k+1; : : : ; yi)			(á).
Здесь		F	некоторая функция, связывающая приближения. В явных схемах приближенное значение yi â
точке		xi	определяется через некоторое число k уже определенных приближенных значений. В неявных
схемах		yi	определяется не рекуррентным образом как в явных схемах, а для его определения возникает

уравнение, поскольку равенство (б) представляет из себя именно уравнение на yi. Явные схемы проще, однако зачастую неявные схемы предпочтительнее.

3.2.1 Получение явных схем

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

	u(x + h) = u(x) + hu0(x) +	h2	u00	(x) + : : : +	hn	u(n)(x) + : : : :
		2			n!

Åñëè u(x) решение задачи (2), то

u0(x

) = f(x ; u )

и, следовательно

u00(x

) =

d

f(x; u)

= f0

(x ; u

) +

i

i

i

i

dx

jxi

x

i i

f(xi; ui)fu0 (xi; ui) . Поступая далее таким же образом, можно выразить все производные

u(k)

через произ-

водные известной функции f(x; u) :

u

= u

+ hf(x

; u

) +

h2

[f0

(x ; u ) + f(x ; u )f0 (x ; u )] + : : : :

(3)

i+1

i

i

i

2

x

i

i

i i

u

i i

Обрывая (3) на том или ином члене, получаем различные явные схемы для вычисления приближенного решения с определенной степенью точности по h.

3.2.2 Схема Эйлера (метод ломаных)

Оставляя в (3) только члены первого порядка по h, получаем приближенное равенство: ui+1 ¼ ui+hf(xi; ui) :

Заменяя в нем точные значения ui = u(xi) на приближения yi, получаем приближенную схему:

8

< y0 = u0

; i = 0; 1; : : : ; N :

: yi+1 = yi + hf(xi; yi)

Указанная процедура является методом Эйлера и имеет пеpвый порядок сходимости по h , åñëè f(x; u) ограничена и ограничены ее первые производные по обоим аргументам. Убедимся в этом. Действительно,

пусть c = max

f

; f0

;

j

f0

. Обозначим разность между истинным решением u

j

в точке x

и найденным

x;u fj

j

j xj

ujg

j

по методу Эйлера приближением yj

через

vj

, тогда

h2

vj+1 = vj + h [f(xj; uj)

¡

f(xj; yj)] +

u00(xj) + O(h3) ;

2

|

fy0

;y˜ )v

j

}

(x{zj j

ãäå y˜j0

некоторая точка между

uj è

yj

. Заметим, что поскольку y0 = u0 , òî

v0

= 0 . Тогда

v1 = 1 h2u00 + O(h3) , и далее

2

0

1

v

= v

(1 + hf0

(x ; y˜0 )) +

h2u00 + (h3) =

2

1

u

1

1

2

1

28

= 12h2 (u001 + u000 [1 + hfu0 (x1; y˜1)]) + O(h3) ;

.

					j	j
			1		X	Y
v	j+1	=		h2	u00	[1 + hf0	(x	; y˜ )] + O(h3) =
		2			k=0	i=k+1

Поскольку u00 = fx0 + ffu0 , òî

jvj+1

= 2h2 k=0 uk00

·1 + i=k+1 hfu0 (xi; y˜i)¸+O(h3) :

j

j

1

X

X

|

{z

}

j

k+1

·c(xj¡xk+1)

u00

c + cc

|

c

·, èf

{z

¡

g

}

exp

c(x

x

)

j

j ·

´

1

j

xj

x0

1

X

1

j ·

2

hc1 k=0 hec(xj¡xk)

=

2

hc1

Z ec(xj¡t)dt + o(h) =

= h2c1c[ec(xj¡x0) ¡ 1] + o(h) = O(h) :

Таким образом, метод Эйлера имеет первый порядок точности по h и при достаточно малом шаге приближенное решение близко к точному.

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

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

Выпишем ряд Тейлора для решения дифференциального уравнения u(x) с точностью до квадратичных

членов

h2

u

= u

+ hf(x ; u

) +

[f0

(x ; u

) + f(x ; u

)f0

(x

; u

)] + : : : :

(4)

j+1

j

j j

2

|

x

j j

j j

u

j

j

}

u00{z(xj)

Сама по себе такая схема уже годится для приближенного решения дифференциального уравнения, однако

ее неудобство состоит в том, что приходится дифференцировать функцию f(x; u)	по обоим аргументам.
Если заменить эти производные разностными, то формально можно записать
uj+1 = uj + h[®f(xj; uj) + ¯f(xj + °h; uj + ±h)] + : : : ;	(5)

где константы ®; ¯; °; ± необходимо определить исходя из того, что эти два представления должны совпадать с точностью до O(h3) . Для этого разложим в (5) f(xj + °h; uj + ±h) в ряд Тейлора:

uj+1 = uj + h(® + ¯)f(xj; uj) + ¯h2[°fx0 (xj; uj) + ±fu0 (xj; uj)] + O(h3) :

Сравнивая с (4), получаем 3 уравнения на 4 неизвестных коэффициента: ® + ¯ = 1 , ¯° =					1	,	¯± =
					2
1 f(xj; uj) . Выразив их через ¯ и заменив истинные значения uj = u(xj) на приближенные yj					и отбросив
2
кубические члены, получаем набор разностных схем Рунге-Кутта 2-го порядка
	h		h
yj+1 = yj + h[(1 ¡ ¯)f(xj; yj) + ¯f(xj +		; yj +		f(xj; yj))] ; 0 < ¯ · 1 :
	2¯		2¯

Обычно полагают ¯ равным 1/2 или 1.

29

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

Изложенным выше способом можно строить схемы типа Рунге-Кутта различного порядка точности по h . Â

частности, метод Эйлера является схемой Рунге-Кутта 1-го порядка. Наиболее удобной и употребительной является схема 4-го порядка. Она имеет следующий вид

	h
yj+1 = yj +		(k1 + 2k2 + 2k3	+ k4) ;
	6
k1 = f(xj; yj) ; k2 = f (xj + h=2; yj + hk1=2) ;
k3 = f (xj + h=2; yj + hk2=2) ; k4 = f(xj + h; yj + hk3) :
На каждом шаге величины km рассчитываются заново.
Интересно отметить, что если f есть функция только от			x , то решение уравнения есть u(x) =
x

u0 + R f(t)dt , и формулы Рунге-Кутта превращаются в формулы приближенного интегрирования. Методу

x0

Эйлера соответствует формула левых прямоугольников, методу Рунге-Кутта 2-го порядка с ¯ = 1 ñîîò-

ветствует формула средних, а с ¯ = 1=2 формула трапеций. Наконец, методу Рунге-Кутта 4-го порядка соответствует формула Симпсона с шагом h=2 . Это косвенно свидетельствует о порядке точности той или

иной схемы.

Естественным образом схемы Рунге-Кутта обобщаются на случай систем уравнений 1-го порядка при помощи формальной замены функций y(x) è f(x; y) на вектор-функции y(x) è f(x; y) . Ïðè ýòîì,

поскольку уравнение n-го порядка эквивалентно системе из n уравнений 1-го порядка, то методы Рунге-

Кутта можно применять к задаче Коши для уравнений порядка выше 1-го. В частности, рассмотрим задачу

Коши для уравнения 2-го порядка

u00 = f(x; u; u0)

>

<

8u(x0) = u0

:

>u0

(x0) = u00

u:

Обозначим u0 = v и введем вектор

u = Ã v ! , тогда система принимает вид

u = f(x; u)

8 dx

à v

! = Ãf(x; u; v) !

>

d

u

v

:

>

<

=

!

>

à v(x0) ! Ãu00

(u(x0) = u0

>

u(x

)

u

>

0

0

>

>

>

yj

:

km

!

Если ввести вектор yj = Ã zj ! приближений к истинному решению uj в точке xj и вектора km = Ã qm

расчетных коэффициентов, то метод Рунге-Кутта 4-го порядка принимает вид

y(j+1) = Ã zj+1

!

= Ã zj

+ h(q1

+ 2q2

+ 2q3

+ q4)=6

!

yj+1

yj

+ h(k1

+ 2k2

+ 2k3

+ k4)=6

h

h

k1 = zj ; k2 = zj +

q1

; k3 = zj +

q2

; k4 = zj + hq3 ;

2

2

h

h

h

q1 = f(xj; yj; zj) ; q2 = f(xj +

; yj +

k1; zj +

q1) ;

2

2

2

h

h

h

q3 = f(xj +

; yj

+

k2; zj +

q2); q4 = f(xj + h; yj + hk3; zj + hq3) :

2

2

2

30