Метод Ейлера за допомогою Excel

У наступній задачі дістанемо чисельні розв’язки деякої задачі Коші спочатку методом Ейлера, потім удосконаленим методомЕйлера і порівняємо отримані результати.

Задача 1. 1) Знайти чисельні розв’язки задачі Коші у´(х) = 2у + sin(3x – y), y(0) = 1 на відрізку [0; 1] з кроками інтегрування h = 0,2 0,1 0,05 0,025 методомЕйлера.

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

Розв’язання. 1)Кроки інтегрування задамо у чарунках H1: H4 так, що H1 = 0,2, а H4 = 0,025. Побудуємо електронну таблицю для чисельного розв’язання даної задачі Коші з кроком інтегрування 0,2 методомЕйлера. Надамо чарункам таких значень:

 

	A	B	C
	x	y	Dy
			= $H$1*(2*B2+0,5*SIN(3*A2-B2))
	= A2+$H$1	= B2+C2	↓
	↓	↓	↓

 

Тут у чарунках А2:В2 координати початкової точки, у чарунці С2 підраховується кутовий коефіцієнт f (х₀; у₀), помножений на крок інтегрування, у стовпці А отримуються значення х_k з кроком 0,2 на відрізку [0; 1]. Тоді у стовпці В будуть обчислені відповідні значення у(х_k) згідно з рекурентною формулою (1). (Як і завжди, символ ↓ тут означає копіювання попередніх чарунок). В результаті отримаємо таку таблицю:

 

	A	B	C
	х	у	Dy
			0,315853
	0,2	1,315853	0,460715
	0,4	1,776568	0,656112
	0,6	2,43268	0,913941
	0,8	3,346621	1,257504
		4,604125	1,741706

 

Звичайно, стовпець А копіюємо лише до чарунки із значенням 1, в даному випадку до А7. Звідси у(0,2) ≈ 1,315853, у(1) ≈ 4,604125. Для отримання аналогічної таблиці з кроком інтегрування 0,1 можна скопіювати попередні формули, наприклад у діапазон Е1:G3, як це зробили ми, а потім у чарунках Е3 і G2 у $H$1 замінити 1 на 2. Таким чином отримуємо таку таблицю:

	E	F	G
	x	y	Dy
			= $H$2*(2*F2+0,5*SIN(3*E2-F2))
	= E2+$H$2	= B2+C2	↓
	↓	↓	↓

 

І тоді в результаті обчислень дістанемо:

 

	E	F	G
	х	у	Dy
			0,157926
	0,1	1,157926	0,193761
	0,2	1,351687	0,236194
	0,3	1,587881	0,285831
	0,4	1,873712	0,343548
	0,5	2,21726	0,410586
	0,6	2,627846	0,488745
	0,7	3,116592	0,580802
	0,8	3,697394	0,691336
	0,9	4,38873	0,828093
		5,216823	1,003441

 

	A	B	C
	х	у	Dy
			0,078963
	0,05	1,078963	0,087871
	0,1	1,166835	0,097626
	0,75	3,583389	0,33404
	0,8	3,917429	0,366778
	0,85	4,284207	0,403754
	0,9	4,687961	0,44594
	0,95	5,133901	0,494482
		5,628383	0,550564

Так само скопіюємо попередні формули у два нових вільних діапазони, (скажімо у А9:С11 та J1:L3, як це зробили ми), а потім у відповідних чарунках замінимо в $H$1 1 відповідно на 3 і 4. В результаті дістанемо:

	J	K	L
	х	у	Dy
			0,039482
	0,025	1,039482	0,041702
	0,875	4,642152	0,220832
	0,9	4,862984	0,232778
	0,925	5,095762	0,245642
	0,95	5,341404	0,259503
	0,975	5,600907	0,274432
		5,875339	0,290478

 

 

Тепер побудуємо електронну таблицю для чисельного розв’язання даної задачі Коші з кроком 0,2 удосконаленим методомЕйлера. Надамо чарункам таких значень:

 

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

 

Формули у стовпцях А, В, С майже ідентичні відповідним формулам таблиці методу Ейлера, лише чарунці В3 надана формула = B2 + F2 замість = B2+C2 у методі Ейлера. Але в цьому й врахована відмінність цих методів, у стовпці В обраховуються значення згідно з (2): у_k+1 = у_k + f (х_k + ; у_{k + ½})h, а доданок f (х_k + ; у_{k + ½})h підраховується в стовпці F. Його аргументи х_k + і у_{k + ½} = у_k + f (х_k; у_k) обраховуються відповідно у стовпцях D і Е. Тому формулу у чарунці F2 можна просто скопіювати з чарунки С2. В результаті обрахунків за попередньою таблицею дістанемо:

 

	A	B	C	D	E	F
	x	y	Dy	x1	y1	Dy1
			0,315853	0,1	1,157926	0,387522
	0,2	1,387522	0,484148	0,3	1,629596	0,585181
	0,4	1,972703	0,719274	0,5	2,33234	0,858985
	0,6	2,831688	1,046859	0,7	3,355118	1,246989
	0,8	4,078677	1,532052	0,9	4,844703	1,853903
		5,93258	2,352282	1,1	7,108721	2,905362

 

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

 

	J	K	L
	x	y	Dy
			= $H$2*(2*K2+0,5*SIN(3*J2-K2))
	= J2 + $H$2	= K2 + O2	↓
	↓	↓	↓

 

	M	N	O
	x1	y1	Dy1
	= J2 + 0,5*$H$2	= K2 + 0,5*L2	= $H$2*(2*N2+0,5*SIN(3*M2-N2))
	↓	↓	↓
	↓	↓	↓

 

І тоді в результаті обчислень дістанемо:

 

	J	K	L	M	N	O
	x	y	Dy	x1	y1	Dy1
			0,157926	0,05	1,078963	0,175743
	0,1	1,175743	0,196748	0,15	1,274116	0,218126
	0,2	1,393869	0,24312	0,25	1,515429	0,268443
	0,3	1,662312	0,297933	0,35	1,811278	0,327763
	0,4	1,990075	0,362495	0,45	2,171323	0,397662
	0,5	2,387737	0,438765	0,55	2,60712	0,480547
	0,6	2,868285	0,529838	0,65	3,133204	0,58035
	0,7	3,448634	0,640956	0,75	3,769112	0,703889
	0,8	4,152523	0,781328	0,85	4,543187	0,863032
	0,9	5,015555	0,966349	0,95	5,49873	1,076088
		6,091644	1,215832	1,05	6,69956	1,359749

 

Так само скопіюємо попередні формули у два нових вільних діапазони, а потім замінимо в $H$1 1 відповідно на 3 і 4 у відповідних чарунках. В результаті отримуємо таблиці:

 

	A	B	C	D	E	F
	x	y	Dy	x1	y1	Dy1
			0,078963	0,025	1,039482	0,083404
	0,05	1,083404	0,088249	0,075	1,127529	0,09313
	0,1	1,176535	0,09844	0,125	1,225755	0,103781
	0,15	1,280316	0,109578	0,175	1,335105	0,115401
	0,2	1,395717	0,121713	0,225	1,456574	0,128047
	0,25	1,523765	0,134906	0,275	1,591218	0,141786
	0,3	1,665551	0,149232	0,325	1,740167	0,1567
	0,35	1,822251	0,164781	0,375	1,904642	0,172889
	0,4	1,99514	0,181665	0,425	2,085972	0,190473
	0,85	4,5881	0,43649	0,875	4,806345	0,460151
	0,9	5,048251	0,487008	0,925	5,291755	0,514551
	0,95	5,562803	0,545886	0,975	5,835746	0,577854
		6,140657	0,614042	1,025	6,447678	0,650494

та

	J	K	L	M	N	O
	x	y	Dy	x1	y1	Dy1
			0,039482	0,0125	1,019741	0,04059
	0,025	1,04059	0,04175	0,0375	1,061465	0,042913
	0,05	1,083503	0,044129	0,0625	1,105568	0,045348
	0,075	1,128851	0,046621	0,0875	1,152162	0,047897
	0,1	1,176749	0,049229	0,1125	1,201363	0,050564
	0,125	1,227312	0,051956	0,1375	1,25329	0,05335
	0,15	1,280662	0,054803	0,1625	1,308064	0,056259
	0,175	1,336921	0,057776	0,1875	1,365809	0,059295
	0,2	1,396216	0,060877	0,2125	1,426655	0,062461
	0,85	4,595439	0,218654	0,8625	4,704766	0,224559
	0,875	4,819998	0,230857	0,8875	4,935427	0,237228
	0,9	5,057226	0,244032	0,9125	5,179242	0,250911
	0,925	5,308137	0,258262	0,9375	5,437268	0,265687
	0,95	5,573824	0,27362	0,9625	5,710634	0,281618
	0,975	5,855442	0,290152	0,9875	6,000518	0,298734
		6,154175	0,307866	1,0125	6,308108	0,317014

 

Такими є результати обрахунків. Оскільки похибка на кожному кроці є систематичною, то найбільшу похибку слід очікувати у кінцевій точці 1. Тому порівняємо тепер наближені значення розв’язку задачі Коші, отримані в задачі методом Ейлера та удосконаленим методомЕйлера саме у цій точці за допомогою наступної таблиці.

 

	Метод Ейлера	Удосконалений метод Ейлера
h
0,2	4,604125		5,93258
0,1	5,216823	0,612698	6,091644	0,159064
0,05	5,628383	0,41156	6,140657	0,049013
0,025	5,875339	0,246957	6,154175	0,013518

 

Тут у стовпці h використані в попередній задачі кроки інтегрування диференціального рівняння. – це наближене значення розв’язку задачі Коші в точці 1, отримане методом Ейлера або удосконаленим методомЕйлера з кроком h і взяте із відповідної попередньої таблиці. – це – , наприклад, для методу Ейлера = 0,612698 = – = 5,216823 – 4,604125. З таблиці відразу видно, що значення для удосконаленого методуЕйлера відчутно менші відповідних значень методуЕйлера, тобто збіжність удосконаленого методуЕйлера значно швидша. Якщо подивитись більш уважно, то можна помітити, що для методуЕйлера значення зменшуються приблизно пропорційно зменшенню кроку h, а для удосконаленого методуЕйлера квадрату зменшення h. Точні оцінки збіжності цих методів, апріорні та апостеріорні отримаємо у двох наступних параграфах відповідно.

 

Методи Рунге – Кутта

Природно запитати: якщо в рекурентній формулі удосконаленого методуЕйлера

у_{k + ½} = у_k + f (х_k; у_k); у_k+1 = у_k + f (х_k + ; у_{k + ½})h замість записати : у_{k + ⅓} = у_k + f (х_k; у_k); у_k+1 = у_k + f (х_k + ; у_{k + ⅓})h, то це теж призведе до значного збільшення швидкості у порівнянні з методомЕйлера? Виявляється, що ні, тобто вибір параметру α = ½, а не α = ⅓ є вирішальним. Загальна відповідь на питання про найкращі параметри заснована на наступних формулах Рунге – Кутта.

У найбільш загальному вигляді задача ставиться так. Нехай відомим є значення y_k = у(х_k) розв’язку задачі Коші (1) в довільній точці х_k, h – крок інтегрування, треба наближено обчислити у(х_k+1) – точне значення розв’язку в наступній точці х_k+1 = х_k + h. В процесі обчислень фіксовані деякі числа (параметри) p₁, …, p_q, α₂, …, α_q, β_ij 0 < j < i ≤ q. Тоді

Означення 3. Формули для обчислення наближеного значення у_k+1 розв’язку в точці х_k+1

у_k+1 = y_k + , де

w₁(h) = h f (x_k,y_k), w₂(h) = h f (х_k + α₂h, y_k + β₂₁w₁(h)),

…………………………………………………………………………………

w_q(h) = h f (x_k + α_qh, y_k + β_q1w₁(h) + β_q2w₂(h) + … + β_qq-1w_q-1(h))

називають формулами Рунге – Кутта. Відповідні ітераційні методи за цими рекурентними формулами називають методами Рунге – Кутта.

У цих формулах враховані практично всі можливі удосконалення методу Ейлера і треба так підібрати параметри, щоб збіжність відповідного методу була найшвидша можлива. Спочатку надамо точного змісту словам про найшвидшу збіжність методу. Величина φ_k+1(h) = у_k+1 – у(х_k+1) = y_k + – у(х_k + h) – це похибка наближеного значення у_k+1 в точці х_k+1. Якщо f (x, y) є достатньо гладкою функцією своїх аргументів, то w₁(h), w₂(h), …, w_q(h) і φ_k+1(h) – такі ж гладкі функції свого аргументу h.

Означення 4. 1.Величину φ_k(h) = у_k+1 – у(х_k+1) = y_k + – у(х_k + h) називають похибкою методу Рунге – Кутта на кроці.

2. Якщо при даних фіксованих параметрах φ_k(0) = φ_k´(0) = … = φ_k^(s)(0) = 0 для всіх достатньо гладких f (x, y), що визначають відповідну задачу Коші (1), а φ_k^(s+1)(0) ≠ 0 для деякої гладкої f (x, y), то натуральне число s називають порядком точності методу.

Іншими словами, s є порядком точності методу, якщо похідні наближення у_k+1 = y_k + завжди співпадають з похідними точного розв’язку у(х_k + h) при h = 0 до порядку s включно, а похідні порядку s + 1 для деяких задач Коші (1) уже не рівні між собою. Згідно з формулою Тейлора

φ_k(h) = hⁱ + h^s+1 = h^s+1 (0 < θ < 1). (4)

Отже, h^s+1 ≤ ôφ_k(h)ô ≤ h^s+1, тобто величина похибки методу на кроці має порядок s + 1: φ_k(h) = О(h^s+1). Але з переходом від кроку до кроку ця похибка систематично зростає і на відрізку завдовжки одиниця сумарна похибка методу є O(h^s), оскільки ітерацій з кроком інтегрування h на такому відрізку порядку . Таким чином, порядок точності методу завжди на одиницю менший, ніж порядок його похибки на кроці. Так порядок точності методу Ейлера дорівнює одиниці, бо на одному кроці похибка цього методу має порядок h², а порядок точності удосконаленого методу Ейлера дорівнює двом, оскільки на одному кроці похибка цього методу має порядок h³.

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

При q = 1 маємо лише один параметр p₁ (оскільки для β_ij 0 < j < i ≤ q, звідки при q = 1 0 < j < 1, а такого натурального j не існує). Тому у_k+1 = y_k + p₁w₁(h) = y_k + p₁h f (x_k,y_k), φ_k(h) = у_k+1 – у(х_k+1) = y_k + p₁h f (x_k,y_k) – у(х_k + h). Звідси φ_k(0) = 0; φ_k´(0) = p₁ f (x_k, y_k) – у´(х_k + h)׀_h=0 = (p₁ – 1) f (x_k,y_k), оскільки у´(х_k + h)׀_h=0 = у´(х_k) = f (x_k,y_k); φ_k´´(0) = – у´´(х_k). Очевидно, що рівність φ_k´(0) = 0 для всіх f (x_k,y_k) буде виконана лише за умови, що p₁ = 1. За такої умови у_k+1 = y_k + h f (x_k,y_k), що співпадає з (1), тобто це формула методу Ейлера. Отже,

 

 

Висновок 1. 1) Метод Ейлера – це метод Рунге – Кутта з q = 1.

2) Метод Рунге – Кутта з q = 1 має перший порядок точності тоді і тільки тоді, коли це метод Ейлера.

3) Похибка методу Ейлера на кроці має вигляд φ_k(h) = – у´´(х_k + θh)h² (0 < θ < 1).

Пункт 3) випливає безпосередньо з (4), оскільки, очевидно, φ_k´´(h) = – у´´(х_k + h), звідки φ_k´´(θh) = – у´´(х_k + θh).

Розглянемо тепер випадок q = 2. Тоді маємо чотири параметри: p₁, p₂, α₂, β₂₁. Тому у_k+1 = y_k + p₁w₁(h) + p₂w₂(h) = y_k + p₁h f (x_k, y_k) + p₂h f (, ), де = x_k + α₂h, = y_k + β₂₁h f (x_k, y_k), φ_k(h) = у_k+1 – у(х_k+1) = y_k + p₁h f (x_k, y_k) + p₂h f (, ) – у(х_k + h), звідки φ_k(0) = 0. Обчислимо похідні функції φ_k(h).

φ_k´(h) = p₁ f (x_k, y_k) + p₂ f (, ) + p₂h ( + ) – у´(х_k + h) = p₁ f (x_k, y_k) + p₂ f (, ) + p₂h ( α₂ + β₂₁ f (x_k, y_k)).

Звідси, оскільки ׀_h=0 = x_k, ׀_h=0 = y_k, у´(х_k) = f (x_k, y_k), то φ_k´(0) = p₁ f (x_k, y_k) + p₂ f (x_k, y_k) – у´(х_k) = p₁ f (x_k, y_k) + p₂ f (x_k, y_k) – f (x_k, y_k) = (p₁ + p₂ – 1) f (x_k, y_k).

Отже, рівність φ_k´(0) = 0 для всіх f (x_k,y_k) буде виконана лише за умови, що p₁ + p₂ – 1 = 0.

φ_k´´(h) = 2p₂ ( α₂ + β₂₁ f (x_k, y_k)) + p₂h ( α₂ + β₂₁ f (x_k, y_k)) – у´´(х_k + h). Звідси, оскільки ׀_h=0 = x_k, ׀_h=0 = y_k, то

φ_k´´(0) = 2p₂ ( α₂ + β₂₁ f (x_k, y_k)) – у´´(х_k).

Але у´´(х_k) = (у´(х_k))´ = f (x, y(х))ôx = x_k, y = y_k = + у´(х_k) = + f (x_k, y_k). Звідси

φ_k´´(0) = 2p₂ ( α₂ + β₂₁ f (x_k, y_k)) – ( + f (x_k, y_k)) = (2p₂α₂ – 1) + f (x_k, y_k)(2p₂β₂₁ – 1).

Отже, рівність φ_k´´(0) = 0 для всіх f (x_k, y_k) буде виконана лише за умов 2p₂α₂ – 1 = 0, 2p₂β₂₁ – 1 = 0. Рівність φ_k´(0) = 0 для всіх f (x_k,y_k) виконана за умови, що p₁ + p₂ – 1 = 0. Рівність φ_k(0) = 0 виконана автоматично завжди.

Висновок 2. 1) Метод Рунге – Кутта з q = 2 має другий порядок точності (тобто φ_k(0) = φ_k´(0) = φ_k´´(0) = 0) тоді і тільки тоді, коли

 

p₁ + p₂ = 1

p₂α₂ = (5).

p₂β₂₁ =

 

2) Якщо p₂ = 1, то звідси p₁ = 0, α₂ = β₂₁ = , у_k+1 = у_k + h f (, ), де = x_k + h, = y_k + h f (x_k, y_k). Це рекурентна формула удосконаленого методу Ейлера (3).

(Тут у_{k + ½} = ). (5) – це система трьох алгебраїчних рівнянь для визначення чотирьох невідомих. Очевидно, p₂ ≠ 0, α₂ ≠ 0, β₂₁ ≠ 0, α₂ = β₂₁. Вона має нескінчену множину розв’язків, що залежить від одного параметра, і кожний її розв’язок дає формули Рунге – Кутта другого порядку точності. Наприклад, якщо покласти p₂ = , то звідси p₁ = , α₂ = β₂₁ = 1, у_k+1 = у_k + h(f (x_k, y_k) + f (, )), де = x_k + h, = y_k + h f (x_k, y_k). Такий варіант методу Рунге – Кутта другого порядку точності називають удосконаленим методом Ейлера – Коші.

Апріорні оцінки похибок чисельних наближень у_k до точних значень розв’язків у(х_k) в точках заданої послідовності х_k, які безпосередньо випливають з (4), як правило, значно завищені. Крім того, іноді оцінити значення дуже важко або й неможливо (тоді, наприклад, коли функцію f (x_k, y_k) задано графічно чи таблично). Всі ці недоліки апріорних оцінок можна усунути, застосувавши метод подвійного перерахунку, який дозволяє отримувати апостеріорні оцінки похибок чисельних наближень так само, як це було зроблено в попередньому розділі.