Тема: методи Рунге – Кутта, кратний перерахунок методів Рунге – Кутта

Лабораторна робота 18.

Завдання.

Контрольні питання.

1. Що таке звичайне диференціальне рівняння, що таке порядок диференціального рівняння?
2. Як ставиться задача Коші для звичайного диференціального рівняння першого порядку?
3. Що означає чисельно розв’язати задачу Коші?
4. Що називають похибкою наближеного значення у_k в точці х_k? Що називають кроком інтегрування?
5. Яка рекурентна формула визначає метод Ейлера?
6. Яка рекурентна формула визначає удосконалений метод Ейлера?
7. Чи значно відрізняються похибки методу Ейлера і удосконаленого методу Ейлера?
8. Як пов’язані між собою величина похибки методу та його збіжність?

Задача. 1) Знайти чисельні розв’язки наступної задачі Коші на відрізку [0; 1] з кроками інтегрування h = 0,2 0,1 0,05 0,025 методомЕйлера.

2) Знайти удосконаленим методомЕйлера чисельні розв’язки тої ж задачі Коші з тими ж кроками інтегрування.

3) Порівняти отримані результати.

Варіант	Диференціальне рівняння	Початкова умова
	у´ = e^xcosx	y(–1) = 1
	у´ = x^2e^-x	y(1) = 5
	у´ = x^2arctgx	y(1) = 3
	у´ = e^sinxsin2x	y(1) = – 1
	у´ =	y(5) = 10
	у´ = e^x/(1+x)	y(0) = 4
	у´ = (1+x)^3/2e^-x	y(1) = 0,5
	у´ = e^x/(3+2cosx)	y(0) = 2
	у´ = e^cosx(1+x)	y(–1) = 0
	у´ = (1+x)ln^2(1+x)	y(0) = 1

(2г.)

Мета: Отримати відомості про методи Рунге – Кутта другого порядку точності,

про апостеріорні методи оцінки точності методів Рунге – Кутта

та навчитися застосовувати ці методи до конкретних задач.

Теоретичні відомості.

Нехай φ_k(h) = у_k – у(х_k) – це похибка наближеного значення у_k в точці х_k для розв’язку задачі Коші у(х) даним методом інтегрування (k = 1, 2, …, n). Тоді

Означення 1. Порядок точності методу інтегрування диференціального рівняння відносно його кроку h дорівнює натуральному числу s, якщо ch^s+1 ≤ ô φ_k(h)ô ≤ Сh^s+1 для всіх k = 1, 2, …, n і всіх достатньо малих h при деяких сталих С, с > 0.

Збіжність методу буде тим швидшою, чим більшим є порядок його точності. Так порядок точності методу Ейлера дорівнює одиниці, а порядок точності удосконаленого методу Ейлера дорівнює двом.

Означення 2. Методи Рунге – Кутта другого порядку точності визначаються рекурентними формулами

у_k+1 = y_k + p₁h f (x_k, y_k) + p₂h f (,), (1)

де

= x_k + α₂h, = y_k + β₂₁h f (x_k, y_k), (2)

за умови, що чотири параметри: p₁, p₂, α₂, β₂₁ зв’язані системою алгебраїчних рівнянь

p₁ + p₂ = 1

p_{2 ×}× α₂ = (3).

p₂ × β₂₁ = .

Очевидно, p₂ ≠ 0, α₂ ≠ 0, β₂₁ ≠ 0, α₂ = β₂₁. Наприклад, якщо p₂ = 1, то звідси p₁ = 0, α₂ = β₂₁ = , у_k+1 = у_k + h f (,), де = x_k + h, = y_k + h f (x_k, y_k).

Це рекурентна формула удосконаленого методу Ейлера (формула (3) з лабораторної роботи 17: тут у_{k + ½} = ).

Застосуємо далі апостеріорні (тобто отримані після і в результаті розрахунків) методи оцінки похибок чисельних наближень у_k до точних значень розв’язків у(х_k) в заданих точках х_k.

Нехай задані послідовність х₀, х₁, …, х_n рівновіддалених точок (х_k = х₀ + kh) з кроком інтегрування диференціального рівняння h, задача Коші у´ = f (x, y); у(х₀) = у₀ та її чисельний розв’язок у₀, у₁, …, у_n у заданих точках х₀, х₁, …, х_n, отриманий деяким методом порядку точності s.

Нехай ще задані послідовність точок , , …, з кроком h/2 так, що = х_k (k = 0, 1, …, n), та ж задача Коші у´ = f (x, y); у(х₀) = у₀ та її чисельний розв’язок , , …, у заданих точках , , …, , отриманий тим самим методом порядку точності s.

Метод подвійного перерахунку для методів Рунге – Кутта є узагальненням такого методу для інтегрування функцій і ґрунтується на двох аналогічних формулах.

1. ε = (правило Рунге) (4)

2. у(х_k) ≈ + ε = + (формула екстраполяції за Річардсоном). (5)

Тут ε – це оцінка похибки наближеного значення у_k в точці х_k, + ε – таке нове наближення до точного значення розв’язку у(х_k), порядок точності якого вже s +1. Процес можна продовжити і далі, отримуючи наближення порядку s + 2, s + 3, …. Це складає апостеріорний метод кратного перерахунку, який є узагальненням методу подвійного перерахунку.

Хід роботи.

Задача 1. Знайти чисельні розв’язки задачі Коші на відрізку [1; 2] з кроками інтегрування h = 0,1 0,05 методом Рунге – Кутта другого порядку точності по формулам , де , і оцінити похибку отриманого розв’язку методом подвійного перерахунку.

Розв’язання. Згідно з означенням2, задана рекурентна формула визначає метод Рунге – Кутта другого порядку точності за умови, що його параметри задовольняють (3). Тут , тож α₂ = , звідки β₂₁ = , р₂ = 1, р₁ = 0. Отже, , де , = y_k + h f (x_k, y_k). Це удосконалений метод Ейлера. Кроки інтегрування задамо у чарунках H1, H2: тут H1 = 0,1, H2 = 0,05. Побудуємо електронну таблицю для чисельного розв’язання даної задачі Коші з кроком інтегрування 0,1 по таким формулам. Надамо чарункам таких значень:

	A	B	C
	x	y	w1
			= $H$1*(2*A2+3*SIN(B2^2))
	= A2 + $H$1	= B2 + F2	↓
	↓	↓	↓
	D	E	F
	x1	y1	w2
	= A2 + 0,5*$H$1	= B2 + 0,5*C2	=$H$1*(2*D2+3*SIN(E2^2))
	↓	↓	↓
	↓	↓	↓

Тут у стовпці С значення w₁(h) = h f (x_k, y_k), у стовпці F значення w₂(h) = h f (,).

Формулу у чарунці F2 можна просто скопіювати з чарунки С2. В результаті дістанемо:

	A	B	C	D	E	F
	х	у	w1	х1	у1	w2
			0,2	1,05	0,1	0,213
	1,1	0,213	0,233606	1,15	0,329803	0,262567
	1,2	0,475567	0,307272	1,25	0,629203	0,365691
	1,3	0,841257	0,455029	1,35	1,068772	0,542874
	1,4	1,384131	0,56232	1,45	1,665291	0,398036
	1,5	1,782167	0,289644	1,55	1,926989	0,147683
	1,6	1,92985	0,154909	1,65	2,007304	0,097318
	1,7	2,027167	0,092905	1,75	2,07362	0,075162
	1,8	2,10233	0,072751	1,85	2,138705	0,072866
	1,9	2,175195	0,080055	1,95	2,215223	0,095675
		2,270871	0,129148	2,05	2,335445	0,188847

Як і в лабораторній роботі 1, можна скопіювати попередні формули, наприклад, у діапазон J1:O3, а потім у чарунках A3, C2, D2, F2 у $H$1 замінити 1 на 2, звідки дістанемо:

	J	K	L
	x	y	w1
			= $H$2*(2*J2+3*SIN(K2^2))
	= J2 + $H$2	= K2 + O2	↓
	↓	↓	↓
	M	N	O
	x1	y1	w2
	= J2 + 0,5*$H$2	= K2 + 0,5*L2	= $H$2*(2*M2+3*SIN(N2^2))
	↓	↓	↓
	↓	↓	↓

В результаті обчислень (деякі рядки в наступній таблиці випущені) дістанемо:

	J	K	L	M	N	O
	х	у	w1	х1	у1	w2
			0,1	1,025	0,05	0,102875
	1,05	0,102875	0,106587	1,075	0,156169	0,111158
	1,1	0,214033	0,116869	1,125	0,272467	0,123626
	1,15	0,337658	0,132065	1,175	0,403691	0,141837
	1,9	2,196925	0,040975	1,925	2,217413	0,045627
	1,95	2,242552	0,052457	1,975	2,268781	0,061468
		2,30402	0,075872	2,025	2,341956	0,095061

Порівняємо отримані наближені значення розв’язку задачі Коші у кінцевій точці 2 за допомогою наступної таблиці подвійного перерахунку:

	A	B	C	D
	h	y(2)	ε	y(2) уточн.
	0,1	2,270871
	0,05	2,30402	= (В17 – В16)/3	= В17 + С17

Тут у стовпці А крок інтегрування, у стовпці В наближені значення в точці 2, отримані у попередніх таблицях з відповідним кроком. У С17 підраховується за правилом Рунге (4) похибка ε, у чарунці D17 уточнене наближення для y(2) за формулою (5). Дістанемо:

	A	B	C	D
	h	y(2)	ε	y(2) уточн.
	0,1	2,270871
	0,05	2,30402	0,01105	2,31507

Ця таблиця і є відповіддю задачі 1.

Задача 2. Знайти чисельний розв’язок задачі Коші на відрізку [0; 1] методом Рунге – Кутта другого порядку точності за формулою , де з точністю 10^-4, оціненою методом кратного перерахунку.

Розв’язання. Згідно з означенням2, задана рекурентна формула визначає метод Рунге – Кутта другого порядку точності за умови, що його параметри задовольняють (3). За умовою тут α₂ = , звідки β₂₁ = , р₂ = , р₁ = . Отже,

, де , = y_k + h f (x_k, y_k).

Спочатку задамо кроки інтегрування 0,2 0,1 0,05 0,025 у чарунках H1:H4, знайдемо чисельні розв’язки для таких кроків і оцінимо їх методом кратного перерахунку. Якщо потрібна точність не буде досягнута, то крок буде зменшено. Надамо чарункам таких значень (аналогічно задачі 1):

	A	B	C
	x	y	w1
			= $H$1*(SIN(0,5*A2+2*B2^2)+1,5*B2)
	= A2 + $H$1	= B2 + 1/4*C2 + 3/4*F2	↓
	↓	↓	↓
	D	E	F
	x1	y1	w2
	= A2 + 2/3*$H$1	= B2 + 2/3*C2	= $H$1*(SIN(0,5*D2+2*E2^2)+1,5*E2)
	↓	↓	↓
	↓	↓	↓

В результаті дістанемо:

	A	B	C	D	E	F
	х	у	w1	х1	у1	w2
			0,481859	0,133333	1,32124	0,315474
	0,2	1,35707	0,287412	0,333333	1,548678	0,270875
	0,8	3,114989	1,097962	0,933333	3,846964	0,958901
		4,108655	1,290827	1,133333	4,969206	1,429217

Скопіюємо А2:F3 у будь – які вільні від інформації діапазони, а потім у відповідних чарунках у $H$1 замінимо 1 на 2, 3 або 4. За цими обрахунками створимо наступну електронну таблицю кратного перерахунку. у кінцевій точці 1 (оскільки саме в ній слід очікувати найбільшу похибку):

	H	I	J	K	L
	№ перерахунку
	h	y(1)	ε	y(1)	ε
	0,2	4,108655
	0,1	3,971733	= 1/3*(I37–I36)	= I37+J37
	0,05	4,056332	↓	↓	= 1/7*(K38–K37)
	0,025	4,051298	↓	↓	↓

Тут у стовпці H крок інтегрування, у стовпці I наближені значення в точці 1, отримані при попередніх підрахунках з відповідним кроком, у стовпці J похибка ε, підрахована за правилом Рунге (4) ε ≈ , де s – порядок точності методу. В даному разі 2^s – 1 = 3, бо s = 2. На цьому закінчився перший перерахунок. У стовпці К знаходимо уточнені наближення порядку точності 3 за формулою (5), у стовпці L їх оцінки ε знову за правилом Рунге, але тепер 2^s – 1 = 7, бо s = 3. У сукупності це другий перерахунок. Далі аналогічно проводимо третій перерахунок: із зростанням номеру перерахунку N на одиницю порядок точності методу s теж зростає на одиницю, отже тут s = 4, 2^s – 1 = 15:

	M	N	O	P
	y(1)	ε	y(1)	ε
	= K38 + L38
	↓	= 1/15*(M39 – M38)	= M39 + N39

В результаті отримаємо таку трикутну таблицю (матрицю):

	H	I	J	K	L	M	N	O
	№
	h	y(1)	ε	y(1)	ε	y(1)	ε	y(1)
	0,2	4,108655
	0,1	3,971733	-0,04564	3,926093
	0,05	4,056332	0,0282	4,084532	0,022634	4,107166
	0,025	4,051298	-0,00168	4,04962	-0,00499	4,044633	-0,00417	4,040464

Згідно з отриманими апостеріорними оцінками тут всі наближення мають якнайбільше три значущих цифри, що є недостатньою точністю за умовою задачі. Отже, знайдемо чисельний розв’язок даної задачі Коші з кроком інтегрування h = 0,0125 і додамо до таблиці кратного перерахунку отримане значення у(1). За підрахунками:

	H	I	J	K	L	M
	х	у	w1	х1	у1	w2
			0,030116	0,008333	1,020077	0,030008
	0,0125	1,030035	0,029921	0,020833	1,049983	0,029679
	0,975	3,90607	0,068214	0,983333	3,951547	0,077847
	0,9875	3,981509	0,083468	0,995833	4,037155	0,088123
		4,068469	0,086471	1,008333	4,126116	0,077411

Отже, достатньо занести отримане значення у(1) = 4,068469 у чарунку I40 поряд з відповідним значенням кроку h = 0,0125 у чарунці H40, решту формул у стовпцях J:O можна просто скопіювати:

	H	I	J	K	L	M	N	O	P	Q	R
	№
	h	y(1)	ε	y(1)	ε	y(1)	ε	y(1)	ε	y(1)	ε
	0,2	4,108655
	0,1	3,971733	*	*
	0,05	4,056332	↓	↓	*	*
	0,025	4,051298	↓	↓	↓	↓	*	*
	0,0125	4,068469	↓	↓	↓	↓	↓	↓	=(O40-O39)/31	=O40+P40

Тут символ * означає, що у відповідній чарунці знаходиться та ж сама формула, що й у попередній таблиці. Для N = 4 s = 5, 2^s – 1 = 31, тому Р40 = (O40 – O39)/31. В результаті отримаємо таку трикутну таблицю:

	H	I	J	K	L
	№
	h	y(1)	ε	y(1)	ε
	0,2	4,108655
	0,1	3,971733	-0,04564	3,926093
	0,05	4,056332	0,0282	4,084532	0,022634
	0,025	4,051298	-0,00168	4,04962	-0,00499
	0,0125	4,068469	0,005724	4,074192	0,00351

	M	N	O	P	Q
	y(1)	ε	y(1)	ε	y(1)
	4,107166
	4,044633	-0,00417	4,040464
	4,077703	0,002205	4,079907	0,001272	4,08118

На четвертому перерахунку наближення досягло вже чотирьох значущих цифр, проте точності 10^-4, яку вимагає умова задачі, ще не досягнуто. Тому знайдемо ще й чисельний розв’язок даної задачі Коші з кроком інтегрування h = 0,0125 і отримане значення у(1) додамо до таблиці кратного перерахунку. За підрахунками:

	A	B	C	D	E	F
	х	у	w1	х1	у1	w2
	0,9625	3,846446	0,02997	0,966667	3,866426	0,030876
	0,96875	3,877096	0,03158	0,972917	3,898149	0,033341
	0,975	3,909996	0,0345	0,979167	3,932996	0,036938
	0,98125	3,946325	0,038371	0,985417	3,971905	0,040913
	0,9875	3,986602	0,042126	0,991667	4,014686	0,043685
	0,99375	4,029897	0,04403	0,997917	4,05925	0,043635
		4,073631	0,042962	1,004167	4,102272	0,040882

Отже, отримане значення у(1) = 4,068469 заносимо до таблиці кратного перерахунку. В результаті отримаємо таку таблицю:

	H	I	J	K	L
	№
	h	y(1)	ε	y(1)	ε
	0,2	4,108655
	0,1	3,971733	-0,04564	3,926093
	0,05	4,056332	0,0282	4,084532	0,022634
	0,025	4,051298	-0,00168	4,04962	-0,00499
	0,0125	4,068469	0,005724	4,074192	0,00351
	0,00625	4,073631	0,001721	4,075352	0,000166

	M	N	O	P	Q	R	S
	y(1)	ε	y(1)	ε	y(1)	ε	y(1)
	4,107166
	4,044633	-0,00417	4,040464
	4,077703	0,002205	4,079907	0,001272	4,08118
	4,075518	-0,00015	4,075372	-0,00015	4,075226	-9,45052E-05	4,075131

Як бачимо, тепер необхідна точність досягнута вже на другому перерахунку. На п’ятому перерахунку досягнуто ще більшої точності.

Відповідь: y(1) = 4,075131 ± 9,45052E-05.

Дата добавления: 2014-10-31