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



Мы поможем в написании ваших работ!


Мы поможем в написании ваших работ!



Мы поможем в написании ваших работ!


ЗНАЕТЕ ЛИ ВЫ?

Метод Ейлера за допомогою 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 підраховується кутовий коефіцієнт f0; у0), помножений на крок інтегрування, у стовпці А отримуються значення х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 + fk + ; уk + ½)h, а доданок fk + ; уk + ½)h підраховується в стовпці F. Його аргументи хk + і уk + ½ = уk + fk; у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 + fk; уk); уk+1 = уk + fk + ; уk + ½)h замість записати : уk + = уk + fk; уk); уk+1 = уk + fk + ; уk + )h, то це теж призведе до значного збільшення швидкості у порівнянні з методомЕйлера? Виявляється, що ні, тобто вибір параметру α = ½, а не α = ⅓ є вирішальним. Загальна відповідь на питання про найкращі параметри заснована на наступних формулах Рунге – Кутта.

У найбільш загальному вигляді задача ставиться так. Нехай відомим є значення yk = у(хk) розв’язку задачі Коші (1) в довільній точці хk, h – крок інтегрування, треба наближено обчислити у(хk+1) – точне значення розв’язку в наступній точці хk+1 = хk + h. В процесі обчислень фіксовані деякі числа (параметри) p1 , … , pq, α2 , …, αq, βij 0 < j < i ≤ q. Тоді

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

уk+1 = yk + , де

w1(h) = h f(xk,yk), w2(h) = h fk + α2h, yk + β21w1(h)),

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

wq(h) = h f(xk + αqh, yk + βq1w1(h) + βq2w2(h) + … + βqq-1wq-1(h))

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

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

Означення 4.1.Величину φk(h) = уk+1 – у(хk+1) = yk + – у(х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 = yk + завжди співпадають з похідними точного розв’язку у(хk + h) при h = 0 до порядку s включно, а похідні порядку s + 1 для деяких задач Коші (1) уже не рівні між собою. Згідно з формулою Тейлора

φk(h) = hi + hs+1 = hs+1 (0 < θ < 1). (4)

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

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

При q = 1 маємо лише один параметр p1 (оскільки для βij 0 < j < i ≤ q, звідки при q = 1 0 < j < 1, а такого натурального j не існує). Тому уk+1 = yk + p1w1(h) = yk + p1h f(xk,yk), φk(h) = уk+1 – у(хk+1) = yk + p1h f(xk,yk) – у(хk + h). Звідси φk(0) = 0; φk´(0) = p1 f(xk, yk) – у´(хk + h)׀h=0 = (p1 – 1) f(xk,yk), оскільки у´(хk + h)׀h=0 = у´(хk) = f(xk,yk); φk´´(0) = – у´´(хk). Очевидно, що рівність φk´(0) = 0 для всіх f(xk,yk) буде виконана лише за умови, що p1 = 1. За такої умови уk+1 = yk + h f(xk,yk), що співпадає з (1), тобто це формула методу Ейлера. Отже,

 

 

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

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

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

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

Розглянемо тепер випадок q = 2. Тоді маємо чотири параметри: p1, p2, α2, β21. Тому уk+1 = yk + p1w1(h) + p2w2(h) = yk + p1h f(xk, yk) + p2h f( , ), де = xk + α2h, = yk + β21h f(xk, yk), φk(h) = уk+1 – у(хk+1) = yk + p1h f(xk, yk) + p2h f( , ) – у(хk + h), звідки φk(0) = 0. Обчислимо похідні функції φk(h).

φk´(h) = p1 f(xk, yk) + p2 f( , ) + p2h ( + ) – у´(хk + h) = p1 f(xk, yk) + p2 f( , ) + p2h ( α2 + β21 f(xk, yk)).

Звідси, оскільки ׀h=0 = xk, ׀h=0 = yk, у´(хk) = f(xk, yk), то φk´(0) = p1 f(xk, yk) + p2 f(xk, yk) – у´(хk) = p1 f(xk, yk) + p2 f(xk, yk) – f(xk, yk) = (p1 + p2 – 1) f(xk, yk).

Отже, рівність φk´(0) = 0 для всіх f(xk,yk) буде виконана лише за умови, що p1 + p2 – 1 = 0.

φk´´(h) = 2p2 ( α2 + β21 f(xk, yk)) + p2h ( α2 + β21 f(xk, yk)) – у´´(хk + h). Звідси, оскільки ׀h=0 = xk, ׀h=0 = yk, то

φk´´(0) = 2p2 ( α2 + β21 f(xk, yk)) – у´´(хk).

Але у´´(хk) = (у´(хk))´ = f(x, y(х))ôx = xk, y = yk = + у´(хk) = + f(xk, yk). Звідси

φk´´(0) = 2p2 ( α2 + β21 f(xk, yk)) – ( + f(xk, yk)) = (2p2α2 – 1) + f(xk, yk)(2p2β21 – 1) .

Отже, рівність φk´´(0) = 0 для всіх f(xk, yk) буде виконана лише за умов 2p2α2 – 1 = 0, 2p2β21 – 1 = 0. Рівність φk´(0) = 0 для всіх f(xk,yk) виконана за умови, що p1 + p2 – 1 = 0. Рівність φk(0) = 0 виконана автоматично завжди.

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

 

p1 + p2 = 1

p2α2 = (5).

p2β21 =

 

2) Якщо p2 = 1, то звідси p1 = 0, α2 = β21 = , уk+1 = уk + h f( , ), де = xk + h, = yk + h f(xk, yk). Це рекурентна формула удосконаленого методу Ейлера (3).

(Тут уk + ½ = ). (5) – це система трьох алгебраїчних рівнянь для визначення чотирьох невідомих. Очевидно, p2 ≠ 0, α2 ≠ 0, β21 ≠ 0, α2 = β21 . Вона має нескінчену множину розв’язків, що залежить від одного параметра, і кожний її розв’язок дає формули Рунге – Кутта другого порядку точності. Наприклад, якщо покласти p2 = , то звідси p1 = , α2 = β21 = 1, уk+1 = уk + h(f(xk, yk) + f( , )), де = xk + h, = yk + h f(xk, yk). Такий варіант методу Рунге – Кутта другого порядку точності називають удосконаленим методом Ейлера – Коші.

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

 



Последнее изменение этой страницы: 2016-04-23; Нарушение авторского права страницы; Мы поможем в написании вашей работы!

infopedia.su Все материалы представленные на сайте исключительно с целью ознакомления читателями и не преследуют коммерческих целей или нарушение авторских прав. Обратная связь - 3.238.132.225 (0.018 с.)