book1989

Приведенные здесь методы являются частными случаями методов Рунге —Кутта третьего и четвертого порядков, рассмотренных по­ дробнее в пп. 3 и 4.

2. Доказательство сходимости. Докажем, что методы Рунге — Кутта сходятся и порядок их точности совпадает с порядком ап­ проксимации.

Выпишем уравнение, которому удовлетворяет погрешность zn = = уп—u(tn). Основное уравнение метода Рунге—Кутта имеет вид

К (y) = f(tn,yn).

Подставим в левую часть уравнения (8) вместо г/, выражения ul+zl при l = n, п+ 1, а в правой части этого уравнения добавим и вычтем сумму

т

t—i
	т,
К (и) = / (fn, Un).	(10)
Тогда уравнение (8) примет вид
■2n" ~ Zn- = № + '№,	(П)
где
	(12)

есть по определению погрешность аппроксимации метода (8), (9) на решении исходной задачи (1) (невязка) и

пг

(13)

i=

Будем рассматривать (11) как уравнение для погрешности ме­ тода. Оно выполняется для п = 0, 1, ... Поскольку начальные зна­ чения г/0 задаются точно, у0= и(0), величина z0 равна нулю. Будем считать, что задача (1) решается на ограниченном отрезке вре­ мени 0<Lt-^iT, и, следовательно, при любых п и т выполняется не­ равенство tn = m ^ .T .

221

Предположим, что в рассматриваемой области изменения пе­ ременных (1, и) функция /(/, и) удовлетворяет условию Липшица по и с константой L, не зависящей от t. При этих предположениях оценим сначала функцию фф2’, а затем и решение 2 „+1 уравнения

( И ) .

Из выражений (9), (10), используя условие Липшица, получим

I h (У) — h (и) | ^ L уп — ип| + ^ т I ьц I I ki (У) — ki (“) l^j .

1=2, 3, . . . . т, | kl{y) - kl(u) | г^Ь \ул- и п\.

Обозначим		г = 1, 2, ....
п = | ki (у) — Mu) I,			т,
Ь = шах |bii |,		g = L\yn— un\.		(14)
Тогда согласно предыдущему неравенству будем иметь
t-1	i — 2, 3, . . ., m,		rt s^g,
n < L b '2 } Tri -\-g,
f=l
или, что то же самое,
i
r-u rx^L b ^xrj + g,	1=	0, 1, .. ,,m —1,	r„ = 0.	(15)
;=o

Л е м м а 1. Из неравенств (15) при Lbx> 0 следуют оценки

	П < p‘-1g, 1=		1, 2,	. ..,	т,		(16)
где р= l + Lbr.		Оценка	(16)	при	1=1 совпадает	с оцен­
Д о к а з а т е л ь с т в о .
кой (15) для	1 = 0. Пусть неравенство			(16) выполнено		для	1 =
= 1, 2, ..., k.	Покажем,	что оно выполнено			и для != /г+1. Из		(15)
при i =k получим		k

rk+ L^Lb^ Tr,- + g.

/—L

Согласно предположению индукции имеем

/=1, 2,

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

run < {bin 2 р'-2 + 1j g = [ibx

+ l ) g = pkg,

что и требовалось.

222

Оценим

теперь функцию

ф())\

определенную

согласно

(13).

Из (14), (16) следует неравенство

П1

Ш

I С I -< 2 1

и 1^

og 2 iJl“1 ^ аётРт~\

i--=1

i^i

где а = max

|сь|, р=1+1Ьт,

g —L\zn\.

Итак, окончательно имеем следующую оценку для ф')Ф

|ф«’ \^ a L m ( l

+ Lbz)m-l \zn\.

(17)

Таким

образом, при

возрастании погрешности

|z„|

величина

|фа2) I растет не быстрее первой степени погрешности.

Теперь

уже

несложно

оценить

погрешность

£„= //„—u(tn). Из

уравнения

(11)

имеем

= г п + тф® + тф»),

откуда, учитывая (17), получаем неравенство

\гПп | <(1 + ат) |г„ | + т | фф |,

п = 0, 1, . . . ,

(18)

где

а = а(т) =аЬт(\+ЬЬт)п~1.

(19)

Заметим,

что a (x)-^oLm

при

т->0.

Если

т ^ т а,

то а ( т ) ^

^oL m eLhi,n-l)x°,

т. е. «(т)

ограничена равномерно по т. В качест­

ве т0с большим загрублепием можно взять Т.

Из неравенства (18) следует оценка

I гп+11^(1 + ост)'1+11гч | +

У т (1 + а%)п~’ | ф'/> |,

(20)

/=О

которую легко доказать по индукции.

что 20= 0,

получим

Загрубляя

оценку (20) и учитывая,

| гл+1 К

(я + 1) т (1 +

ил)п шах

| фф | < /Г1+1еа<'1 шах | ф9> |,

где /,,=/гт^У .

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

часть уравнения

(1) f(t, и)

удо­

Т е о р е м а

1. Пусть

правая

влетворяет условию Липшица по второму аргументу с константой L. Пусть ф'Ц —невязка метода Рунге — Кутта (2), определенная

согласно (12). Тогда для погрешности метода при пт^.Т справед­

лива оценка	шах	\ф<?>|,	(21)
|уп— и (tn) | ^ ТеаТ
где
a = oLm (1-fLb%)m_1,
ст — шах |о;|, Ь =	шах	\Ьц\.
	2«0'<0п
			223

С л е д с т в и е . Если метод Рунге — Кутта аппроксимирует ис­ ходное уравнение, то он сходится при т-»-0, причем порядок точно­

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

следует

из оценки

Доказательство

этого утверждения

сразу

(21) и замечания о равномерной ограниченности а(т).

3. М етоды тр ет ь его

п о р я д к а

точ н ости .

П ри

реш ен и и обы к н овен н ы х

д и ф ф е ­

рен ц и ал ьн ы х

уравн ен и й

ч асто

и сп ол ь зую т ся

м етоды тр ет ь его

и ч ет в ер т ого

п о р я д ­

ка точ н ости . П р и в ед ем

вы в од

так и х

м ето д о в .

С н ач ал а

р ассм отр и м

т р ехэтап н ы й

м е т о д

ki =

Уп),

62 =

f(tn +

апт , уп+

blxxki),

/

(/„ .

63 =

/

( tn

4" a 'JT > Уп4"

4

632т 62) ,

(22)

I/n + l = ( / n + T ( ( J i £ i +

0262 + 0363).

В ы ясн и м ,

каким усл о в и я м

д о л ж н ы у д о в л ет в о р я т ь

п ар ам етр ы at,

Ьц,

at

д л я

того, чтобы

дан н ы й

м е т о д

и м ел

третий п о р я д о к ап п р ок си м ац и и .

П огр еш н ость

ап пр ок си м ац и и м е т о д а (2 2 ) д а е т с я в ы р аж ен и ем

1 — —

т

4

4

0262 4

0363г

(23)

k\ =

/

(tn, Un), kn =

f(tn 4

02-г,

ип+

b,jXk 1) ,

f (tn 4

03т , un+ b3ixkx+ - by/ik..)

(24)

6 3

=

и iin=

u(tn) — р еш ен и е

и сх о д н о го

ур авн ен и я

( 1 ) .

П р и м ен яя р а зл о ж е н и е п о ф о р м у л е Т ей л ор а д л я ф ун к ц и и д в у х п ер ем ен н ы х

и уч и ты вая, что k\ = f(tn, ип), п олучи м

62 =

/ 4

02- Д

4

&21Д

/ ц 4

у

-

l°\fu 4

2a,b.Jfta +

b l f f au] 4

0 ( 4 ) .

(25)

З д е с ь

зн ач ен и я ф ун к ц и и

f(t,

и) и

е е

частн ы х

п р о и зв од н ы х

б ер у т ся

при

/ = / „ ,

« = « п . Т оч н о

так

ж е

63 = f 4 аз Д

4 (631614 63262) Д

4

Н—

7^ Ialftt 4

2а 3 (63161 + -

Ьззб,,) ftu+ -

(b3Xkx-+

63262)2 fua1 4

О(т3) .

П о д ст а в л я я сю д а

6l=/,

б2= /+а2Т((+б21Ди+ 0(хг),

п ол уч и м

63 = f 4 х Ш 1 4

(6314

632) ffа\ 4

lalfи 4

2аз (63 14

632) ffta 4

4

( 6 3 1 4

6 3 2 ) 2 }%u 4 2632а Д

/ и 4

2632&2i / Д ) 2] 4

О ( Д

•

( 26 )

Д а л е е , из р а зл о ж ен и й (2 4 ) — (2 6 ) с л е д у е т

° i 6'i 4

0262 4

0363 = (o i

4

02 4

0з) f 4

4

т[ (о 2а 2 4

0заз) ft 4

(026214

03 (6314

632) ffu\ 4

Н

1(02а2 4

03а\) ftt

+

2 Д а + л

4 0заз (631 4 632)) Д и 4

4

((<JI 6221 4

03 (631 4

6,12) 2 ) П и и

4

2 о 3Ь32а 2/ ( / и 4 2a3b32b2X ( /и)г] 4

О (Д ) •

(27)

224

Получим теперь разложение по степеням т разностного отношения (un+i— —ц„)/г, входящего в выражение для погрешности аппроксимации (23). При этом учтем, что в силу уравнения ( 1 ) справедливы следующие соотношения:

u' = f,

U"=*ft +ffu,

U,"=^, +2ff,t,+ /a^ш+^^+/(^),.

(28)

Тогда будем иметь

' = “я + 2 “я + е и" ^ °

=

= f + ^ V t + ffu)

+

Utt +

Vftu + Pfuu + faft + f (/„)*] +

О (+)•

Отсюда и из (27) получаем следующее разложение выражения для погреш­

ности аппроксимации (2 3 ):

^ 11 = (^ + 0 2 + а 3 - 1) / +

+ 5J" ([2 (ога2 + Оз^з) — 1] f t

+

[2 (02621+03 (631 + 63a)) — 1] //„ } +

"I

{[3 (ооЦ^ + Ojd3) — 1 J ftt

+

[ 6

(о-.а+л +

о3а3 (6 31 +

6 32)) — 2 ] fftu +

+ [3

+ а 3

( 6 3 1

+

6 32)2) —

П / 2/ ии +

+

(6 а36 3.2а2 -

1) f j t +

(6 0 3 6 3 3 6 2 1

-

!) f (f u)2} +

О (+ ).

Приравнивая нулю коэффициенты при т\ / = О, I, 2, получаем условия треть­

его порядка аппроксимации:

Oi + 0 2

+ О3 =

1 ,

0 2 % +

о3а3 =

O2621 + о3 (6 31 +

6 32) =

0,5,

aia2 +

о3а3 =

О0О2621 + а3а3 (6 31 +

6 32) = 0 2 б31

+ а3 (6 3] +

6 3 .)s = — ,

О

03632<22= O3632621 —

•

О

После проведения эквивалентных преобразований эту систему уравнений можно записать в более простом виде:

1	2	, 1	,
0 2 а 2 + о3а3 — , 0 3а2 + а3а3 =
а2 = 6 2 1 , а3 — 6 31 +	6 32 ,	Оз632а2 =	о ,	(29)
о1—1—02—03,	03^=0,

Исключим с помощью (29) из выражений (22) коэффициенты btj. Тогда получим метод

*1 = / (<„. Уя).	= / (/„ +	а2т, !/л +	“aX*i),
/		т (6, — &,) \		(30)
*3 = / ^я + аэТ, Уп + O&h +			---- 1 ,
^Я+1	о2 6 2 + о36 3,	Oi =	1 — 0 2 — 0 ;ь
= 01&1 +
g А. А Самарский, А . В. Гулин				225

который имеет третий порядок аппроксимации при условиях

о2я2 +

• ° 2а 2 "Ь а*а1 ~

• °з Ф 0, До ф 0.

(31)

Таким образом, в общем случае существует двухпараметрическое семейство трехэтапных методов Рунге— Кутта, имеющих третий порядок аппроксимации. Задавая а2 и а3 в качестве свободных параметров, получим из (31)

	\_		3 2 <h
а2 =	3	о3 =		(32)
	Д2 (аэ — Ф)		°s (аз — °2)

Кроме того, система (31) ний, определяемых условиями

Й2 = а3 =

2

с?

I ’

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

	2	3	а 3 Ф 0 ,
	з ,	о 2 + а 3 = —- ,
		4
		3	л ю бое.
II	О	о 2 = — , о3ф 0
		4

семейства реше­

(33)

(34)

1

Например, полагая а2 = — , а з= 1 , получим из (30), (32) следующий метод

третьего порядка аппроксимации:

К =	f (*п>Уп)< k* = f ( t n + - , y			n + - kij ,
h =	/( /„ +	"с, У„— т ^ + гт А а ),		(35)
	Уп+1	Уп	1 . . . . . . .
	--------------	=	~Г (&1 + 4йо +	кг).

т6

4.Методы четвертого порядка точности. Рассмотрим теперь четырехэтапный

метод

* 1 — / Уп) • ka = f (tn+ а2х, уп + bjiTfei),

	k3 = f (ln+	пзТ, уп + b3it*i +				632т*2).	(36)
	кк — f (*п “Н а4т> Уп +			k4llk l +		bi2ikn ~f- bwiks),
	*/nt-i =	*/n+T((Tift| +<T2fe2 + (T3fei+ (Tifti),
Погрешность аппроксимации			метода	(36)	равна по определению
	Ф„ * = —	~	+ Oikj +		о2к2 + а3*3 o4k4,		(37)
где функции k{,	1 = 1 , 2, 3, 4, получаются			из (36)		путем замены у„	на точное ре­
шение U n = U ( f n	) .					аппроксимации, необходимо раз­
Чтобы построить схемы четвертого порядка

ложить функции, входящие в (37), по формуле Тейлора до величин третьего по­

рядка	по т включительно	и	приравнять	нулю коэффициенты при степенях	т п,
п = 0,	1, 2, 3. Необходимо	при	этом учесть	соотношения (28) и аналогичное	вы­

ражение для ulv . Опуская выкладки, приведем систему уравнений, которой долж­ ны удовлетворять коэффициенты метода (36) для того, чтобы данный метод

226

имел четвертый порядок аппроксимации:

<Д+ Ог + °з + 04 = 1.

b-2i = о,,		Ьзх + 6 32 — а^, bn +							(38)
								* 42 -f- *43 =	0 4>
<T2^2 +	о 3а3		а4о4 — ^			>
о,а\ +	а3а:; +		а4а* =		у	,			(39)
а2а* +	ааа* +		а4а* =		у	;
(оз&зг +		0 46 1г) а2 +		о4Ь43а3 =			—	,
							о
(а36 32 +		а4Ь42) а\ +		<т4&43а| =			^	,	(40)
а3Ьз2а2а3 -f- о4*42о2а4 -|- а4й43з3а4 = — ;
								О
04*43*з2а2 = 7Г. •									(41)
Система (38) — (41) состоит			из	одиннадцати				уравнений	и содержит тринад­
цать неизвестных. Выберем в качестве					независимых параметров неизвестные а2
и а3 и выразим остальные величины через эти неизвестные.
Для этого сначала разрешим группу уравнений (39) относительно перемен­
ных <т2, cr-j. сг4. Определитель б этой системы
6	=	а2 а3 а4 (а3—а2) (а4 —<h) (а4 —а3) .							(42)

Это означает, что мы не рассматриваем значения параметров а2, а3, at, удовлет­ воряющие хотя бы одному из условий

02=0,	а3=0,	а4=0,	а2= а3, а2=о4, а3= а4.
Заметим сразу же, что а2ф 0, ст4 =И=0 согласно (41).
Предположим, что 8=4=0. Случай 6 =			0 будет рассмотрен позже.
При 8 ф 0 система (39) имеет следующее решение:
	_ 6 а3а4 — 4 ^ — 4 0 4 + 3			(43)
	*	12а2(а3 —о2)(а4 —а2)
		6а2а4 — 402 — 4а4 ~f- 3		(44)
		12о3 (а3 — а2) (о4 — а3) *
		6о2о3— 4д2 — 4о3 -f- 3		(45)
		12а4(а4— а2) (а4 — 03)
Точно так же, решая систему		(40), получим
офзз =	4а4 — 3			(46)
	24oa (о4 — о3) ’
,
<74*42 =	__2 ( 1 — 2 аа) (а4 — аз) — (3 — 4йз) (а3 — а2)			(47)
		24 (о3 — оа) (о4 — а3) а2
0 4*43 =	1 — 2 а2			(48)
	12 (а3 — а,) а3
			8*	227

Учтем теперь соотношение (41). Прежде всего заметим, что Сз¥=0. Действи­ тельно, при о3= 0 из (44) и (46) получаем

3 о4 = — , о2 = О,

4

но в силу (41) имеем а2ФО. Таким образом, уравнение (41) эквивалентно урав­

нению

(0 3 6 3 2 ) (оА з) а2— 2 4 ‘

Подставляя сюда выражения для 0 3 6 3 2 , 0 4 6 4 3 , сг3 из (44), (46), (48) и приводя подобные члены, получим уравнение а2 ( 1—<24) = 0 , из которого следует, что а4= 1 .

Таким образом, при 6=^=0 система (38) — (41) имеет следующее двухпарамет­ рическое семейство решений:

1	2 о3 — 1	|
@А~~~1I С? - 1

1 2 fl2 (as - f l , ) ( l - f l 1)

_____ —1________

° 3 1 2 о3 (а3 — а2) (1 — а3)

6й2^з— 4оо — 4^3 —}- 3

° 4 = 1 2 ( 1 — Оо) (1 — а3)	’
4Оу— о2 — 5п3	2

642 = — ------------------------ ’

24а4а2 (а3 — а2) (•— аз)

1— 2а2

^43 =	(а3	Г	•
1 2 о 4о 3		— а2)

641 = I -- 642-- 643,

631 = 02— 6з-2,

621 - о2,

Oj = 1 — о2 — О3 — о4.

Здесь, как уже отмечалось, о 3 ¥=0, а 4 ^=0, т. е. а2 =т^0,5, 6 а2а3— 4а2—4а3Ч -3 # 0 . Приведенное выше решение справедливо при 6=^=0, т. е. когда параметры а2, 0 3 ,0 4

удовлетворяют условиям
		а (Ф 0,	i= 2 , 3, 4,		агф а г,		агф а 4, а3ф а 4.
Рассмотрим		систему	(38) — (41)		при	тех	значениях параметров о2,	о3,	о4,
когда 6 = 0.	При	о2= 0	система	не	имеет	решения вследствие (41). При		о3=	0
система (40)	принимает вид
		(0 3 6 3 2 +	0 4 6 4 2 ) а2	=			+ °А » ) “а = "{г ’	(49)
			a4a«a4bi2— 1 »
						О
откуда следует, что			"Г и
			и
			а3632+ 04642 -				—	(50)
							О
228

1

Далее, система (39) при а2 = — , а3= 0 принимает вид

4 0а + 2<т4я4 = 1, 0а + 4о4а4 = — , а2 + 8 а4а4 = 2

и имеет единственное решение

а4 = 1, Oj - : 04 = '

Подставляя эти значения а4, а2, а4 и 0 2 = — в уравнения (49), (50), получаем

А	*in =	1
*42 — „ I

1 2 о3

Кроме того, из (41) имеем

1 * 4 3 1* 24а46з202 = 6а3.

Таким образом, система (38) имеет следующее семейство решений, завися­ щее от параметра сг3 ^=0 :

1

0 2	=	’	а 3 =		0 ,	04	=
	2
	2				1
0 2	=		0 4	=	—
	3				6
		1				3
*32	=		»		&42	2	’
	12а3				_
		1	—	6 а 3 | *31 1			1
*41	=						-
	_	2 ‘					1 2 о , ' ’ Ьа1 = 2 ’ 0 1 = 6 _ ° Э-
Точно так же при условии						а2= а 3 система		(38) — (41) имеет решение
	а2 —Оз —				,	04 —1,
	0 2 =		—	—	03 ,1	0 4	—
			3				6
			1
	^Э-2 =		6 а 3		*42 =		1 — З а 3 , *43 =		З а 3 ,
	и					1	1		1	1
		—	0 ,	* 3 i			*21	=		Ст1=7
	° 4 \					2			2 *
							6 а 3 ’	“

зависящее от параметра а3 =5^ 0 . При а2 = а 4 имеем решение

а2 —04 —If			«3 —„ I			02 — _ — 04> 0з — „ I
					2		о	о
032 —	1	*42	—	—	1	,	643 —	1
	_ .				6а4			-
	8	1				3		За4
*41 —	1 —	-	.	*31 =		„ I	*21 —	U
		6а4				8
01 =

зависящее от параметра с^фО.

229

О п р едел и тел ь	б	о б р а щ а ет ся в н уль ещ е				в д в у х	сл уч ая х:		при	а 4= 0 и a3=-at.
О к азы в ается ,	что	в	эти х	сл у ч а я х	си ст ем а	( 3 8 ) — (4 1 ) не		и м еет		реш ен и я . П уст ь ,
н ап р и м ер , а 4= 0.	Т о гд а		и з	п ервы х	д в у х ур авн ен и й		си стем	ы	(3 9 )	п олучи м

			3а*— 2				2 — 3а,
		0>==	6а2(а3 — а2)	,	Оз	^ 6а3(а3 — а2) ,
П р и этом п о сл ед н ее ур ав н ен и е си стем ы				(3 9 )		п р и в оди т к усл ови ю
			6 д2^з—4д2—4а3+ 3 =				0.	(51)
А н ал оги ч н о		н аходи м , что си стем а		(4 0 ), (4 1 )			р азр еш и м а отн оси тел ьн о	Ь32, &«,
bt3 тольк о при усл ови и			6 Д2 Д3— 6 й2—4а3+ 3 =				0.	(52)
И з (5 1 ),	(5 2 ) н а х о д и м		а2—0, что	н е в о зм о ж н о в си л у (4 1 ). Т очн о так				ж е д о ­
к а зы в ает ся ,	что	не с у щ ес т в у ет реш ен и й		с	аг—а4.

§3. Многошаговые разностные методы

1.Формулировка методов. Для решения задачи Коши

^ =	/ М ) .	^ > 0 . w(0) = uo	(1)
at
введем сетку
(Щ={Д=«т, гс=0, 1, . . .}
с постоянным шагом	т> 0 .	Обозначим через	yn = y(tn), f„ —

= f(tn, у„) функции, определенные на сетке сот. Линейным т-шаго-

вым разностным методом называется система		разностных	урав­
нений
аяУп+ а1Уп-1 + • • • + атУ/1-п ■—	+ &1Д- 1 +	. . ■+ bmfn-	(2)
п = т, т+ 1, ...,
где а,„ bh —числовые коэффициенты, не зависящие от п, k = 0,			1, ...

..., т, причем а„Ф0.

Уравнение (2) следует рассматривать как рекуррентное соот­

ношение, выражающее новое значение yn = y(tn)		через	найденные
ранее значения	уп- 2, ... , уп- т.
Расчет начинается с п = т, т. е. с уравнения
°о Ут + а1Ут- 1 +	• • • + атУп		b j 0.
х	= b0f m + b i / m - i +	. . . +

Отсюда видно, что для начала расчета необходимо задать т на­

чальных значений у0, уи		ym- t■Значение у0 определяется исход­
ной задачей (1),	а именно	полагают у0= и0.		Величины уи у2, ...
. . . , t/m_, можно вычислить,		например, с помощью метода				Рунге —
Кутта. В дальнейш ем будем			предполагать,	что	начальные значе­
ния г/0, уи ..., ym- t	заданы.		что в отличие	от	методов	Рунге —
Из уравнения	(2) видно,

Кутта многошаговые разностные методы допускают вычисление правых частей только в точках основной сетки cot.

230