69314

Однокрокові методи розв’язування диференційних рівнянь

Лекция

Математика и математический анализ

Методи чисельного інтегрування диференціальних рівнянь у залежності від числа використовуваних у формулі (8.8) попередніх значень функції чи її похідної підрозділяються на однокрокові (коли використовується інформація тільки про одну попередню точку)...

Украинкский

2014-10-03

802.5 KB

6 чел.

Лекція 18. Однокрокові методи розв’язування диференційних рівнянь

Основні поняття

1. Звичайним диференціальним рівнянням (ЗДР) називається рівняння виду

 (8.1)

яке з’єднує одну незалежну змінну t, шукану функцію x та її похідні до n–го порядку. Функція x(t) незалежної змінної, яка задовольняє рівнянню (8.1), називається розв’язком цього рівняння, а задача знаходження розв’язків диференціального рівняння називається, інакше, задачею інтегрування диференціального рівняння.

2. Порядком ЗДР називається порядок старшої похідної від шуканої функції, наприклад, x” + t3x’ + x3 = 0 — це ЗДР другого порядку.

3. ЗДР називається лінійним, якщо воно має наступний вигляд:

. (8.2)

Наприклад, ЗДР x” + t3x’ + tx = e - t лінійне рівняння другого порядку, а x” + t3x’ + x3 = 0 — нелінійне.

4. Загальним інтегралом рівняння (8.1) називають функцію Ф(t,x,C1,…,Cn), що

пов’язує незалежну змінну t, шукану функцію x(t) і n констант інтегрування за допомогою рівняння:

, (8.3)

тобто розв’язок x(t) входить у (8.3) неявно, причому кількість констант інтегрування дорівнює порядку ЗДР.

5. Загальним розв’язком ЗДР називається функція

, (8.4)

яка з’єднує незалежну змінну t і n констант інтегрування Ci , тобто розв’язок x(t) визначається явно.

6. Для визначення констант інтегрування Сi задаються додаткові умови числом, яке дорівнює кількості констант інтегрування чи порядку ЗДР.

Якщо в додаткові умови підставити (8.3) і розв’язати отриману систему відносно Сi , а потім підставити в (8.3), то отримаємо частковий розв’язок чи частковий інтеграл x(t) = ψ(t).

Аналогічні процедури з загальним розв’язком (8.4) дають частковий інтеграл ψ(t,x(t)) = 0

7. Якщо всі додаткові умови задаються в одній точці x0, то сукупність ЗДР і додаткові умови називають задачею Коші для розглянутого ЗДР. У цьому випадку додаткові умови називають початковими умовами.

8. Якщо додаткові умови задаються більш ніж в одній точці, то сукупність ЗДР і додаткові умови називають крайовою задачею Коші для розглянутого ЗДР, а додаткові умови називають межовими чи крайовими умовами.

9. Таким чином, задача Коші для ЗДР (8.1) має вигляд:

де x0,x0(1),…,x0(n - 1) — задані числа. Порядок старшої похідної в початкових умовах дорівнює (n - 1) і не перевищує n, де n — порядок ЗДР.

10. Якщо старша похідна в ЗДР виражається як функція t, x і похідних більш низького порядку, то рівняння називають розв’язаним щодо старшої похідної, і воно записується в наступному вигляді: x(n) = f(t,x,x(1),…,x(n - 1))

11. Система рівнянь першого порядку виду :

  (8.5)

тобто розв’язаних щодо похідних від шуканих функцій, називається системою, що має нормальну форму Коші. Часто систему (8.5) для скорочення і зручностей виконання обчислень записують у векторній формі. Введемо наступні векторні позначення:

(8.6)

тоді задача Коші для системи (8.5) набуде вигляду:

. (8.7)

12. Задачу Коші для ЗДР n–го порядку x(n) = f(t,x,x(1),…,x(n - 1)) заміною x(t) = y1(t),
x(t) = y2’(t),…,x(n - 1)(t) = yn - 2’(t) = yn - 1(t) можна звести до задачі Коші для нормальної системи ЗДР :

13. Приведемо формулювання умов, що гарантують існування й одиничність розв’язку задачі Коші (8.7) [ 29 ].

Припустимо, що функції fi, i = 1,2,...n неперервні по всіх аргументах у замкнутій області

З неперервності функцій fi випливає їх обмеженість, тобто існування константи M>0 такої, що усюди в D виконуються нерівності |fi|≤M i = 1,2,...n.

Припустимо, крім того, що в D функції fi, задовольняють умові Ліпшиця по аргументах x1,x2,…,xn тобто

для будь–яких точок (t,x11,x12,…,x1n,) і (t,x21,x22,…,x2n,) області D.

Якщо виконано сформульовані вище припущення, то існує єдиний розв’язок

системи (8.7), визначений при |t0|≤t0 = min(a,b/M) і приймаючий при t = t0 задані початкові значення x0.

8.2. Методи Ейлера і Рунге–Кутта розв’язку задачі Коші для ЗДР

Чисельні методи дозволяють знаходити часткові розв’язки задачі Коші, наприклад, для заданих початкових умов за допомогою рекурентних співвідношень загального вигляду :

(8.8)

у які для оцінки значення функції в поточний момент часу tm використовуються вже відомі значення цієї функції xm - i і її похідних xm - i, отримані в попередні моменти часу tm - i. При цьому замість послідовності точних значень функції x(t1),x(t2)… в обрані моменти часу t1,t2,…обчислюються їх наближення u(t1),u(t2)… з деякою похибкою εm = xm - um в кожен момент часу tm.

Методи чисельного інтегрування диференціальних рівнянь у залежності від числа використовуваних у формулі (8.8) попередніх значень функції чи її похідної підрозділяються на однокрокові (коли використовується інформація тільки про одну попередню точку) і багатокрокові (коли використовується інформація про декілька попередніх точок). До однокрокових методів відносяться методи Ейлера, Хейна, Рунге–Кутта, а до багатокрокових — Адамса–Башфорта, Адамса–Мултона, ГіраБрайтона .

Методи чисельного інтегрування диференціальних рівнянь підрозділяються також на явні і неявні в залежності від значень вагових коефіцієнтів при значеннях функції xm - i і її похідних xm - i у формулі (8.8). Якщо α0 = 0, β0 = 0, тобто шукане значення функції xm чи її похідної xm в момент часу tm містяться тільки в лівій частині виразу (8.8), то мова йде про явні методи розв’язку. Якщо ж α0≠0, β0≠0, то вираз (8.8) , де шукане значення функції xm одночасно міститься в лівій і правій своїх частинах, з урахуванням рівняння (8.7) перетворюється в нелінійне алгебраїчне рівняння, властиве процедурам неявних методів розв’язку.

Порядок методу k чисельного інтегрування диференціальних рівнянь, від якого залежить точність чисельного розв’язку, визначається числом k значень функції
xn,xn - 1,…,xn - k, які використовуються в рекурентній формулі розв’язку (8.8).

При цьому для підвищення порядку однокрокових методів використовується не послідовність значень функції в попередні моменти часу tm - i, а ряд пробних оцінок значень цієї функції, отриманих на відрізку (tm - 1,tm) довжиною в один крок за вихідною інформацією про одну попередню точку.

Похибка методу k–го порядку визначається співвідношенням:

, (8.9)

де h = tm - tm - 1 — часовий крок обчислень, з яким вводиться сітка по змінній t, тобто визначається множина часових точок {tn = , n = 0,1,2,…}, у яких виконуються обчислення.

Важливість задачі Коші в багатьох галузях науки і техніки з’явилася причиною розробки для її розв’язку великої кількості як аналітичних (де це можливо), так і наближених чисельних методів. Останні істотно відрізняються за необхідним обсягом обчислень, за стійкістю розв’язку до вибору і зміни часового кроку h = tm - tm - 1, за реалізованою точністю розв’язку та ін.

Для простоти викладу основних ідей обчислювальних методів розв’язку задач з початковими умовами розглянемо спочатку випадок одного звичайного диференціального рівняння першого порядку:

(8.10)

Введемо по змінній t рівномірну сітку з кроком h, тобто визначимо множину точок {tn = , n = 0,1,2,…}, яку будемо називати сіткою. Позначимо через xnточний розв’язок задачі (8.10), а через un — наближений розв’язок в момент часу tn.

8.2.1. Метод Ейлера

Є найпростішим методом розв’язку задачі Коші, але має невисоку точність, тому на практиці його використовують досить рідко. Однак на прикладі методу Ейлера зручно пояснити способи дослідження чисельних методів. Розкладемо точний розв’язок xn + 1 = x(tn + h) задачі (8.10) в околі вузла сітки tn в ряд Тейлора, отримаємо :

 (8.11)

Залишимо у формулі (8.11) тільки два перших члени, тоді замість точного розв’язку xn + 1 виходить його наближене значення un + 1, що знаходиться за формулою:

, n = 0,1,2,...  (8.12)

Оскільки значення x0 відоме, то, послідовно застосовуючи формулу (8.12), знаходимо наближений розв’язок в наступних точках.

Формулу (8.12) можна отримати геометрично (рис. 8.1). Знаючи значення розв’язку x(t0) у точці t0, відоме по початковій умові (8.10), можна обчислити значення похідної x(t0) = f(t0,x0) в цій точці. Запишемо рівняння дотичної до графіка розв’язку x = x(t) в точці (t0,x0):

.

При досить малому кроці τордината u1 = x0 + hf(t0,x0) цієї дотичної буде мало відрізнятися від ординати x(t1) розв’язку задачі (8.10). Отже, точка (t1,u1) може бути наближено прийнята за нову початкову точку. Через цю точку знову проведемо дотичну і т.д. У підсумку отримаємо формулу (8.12).

Рис. 8.1. Геометрична інтерпретація обчислень по методу Ейлера

8.2.2. Апріорна оцінка похибки

Припустимо, що задача (8.10) розв’язується на відрізку t[t0,T]. Знайдемо похибку наближеного розв’язку, яку позначимо через

n = 0,1,2,... (8.13)

Віднімемо з рівності (8.12) рівність (8.11), отримаємо:

.  (8.14)

Формула (8.14) визначає залежність між похибками наближеного розв’язку на двох сусідніх кроках обчислень. Очевидно, що від кроку до кроку, тобто при багаторазовому застосуванні формули (8.12), можливе накладання похибок, обумовлене формулою (8.14). Отримувана за n кроків похибка називається накопиченою чи глобальною похибкою за відповідне число кроків. Оцінимо цю похибку, використовуючи формулу (8.14).

Спочатку отримаємо формулу для локальної (чи, інакше, покрокової) похибки. При un = xn маємо f(tn,un) - f(tn,un) = 0 і з виразу (8.13) отримуємо наступну формулу для локальної похибки

. (8.15)

Звідси випливає, що локальна похибка має другий порядок відносно кроку h. При подальшому розгляді для нас буде важливо, що локальна похибка методу має головний член, а саме:

 (8.16)

Переважний член асимптотичного розкладання локальної похибки (8.16)

 (8.17)

інакше, головний член локальної похибки для методу Ейлера (8.12) має другий порядок відносно τ. Для похідної х’ (t) від точного розв’язку задачі Коші (8.10) можна отримати наступний вираз:

Підставляючи останнє в формулу (8.17), отримаємо головний член локальної похибки для методу Ейлера у вигляді:

(8.18)

Припустимо, що права частина f(t,x) і її частинні похідні першого порядку по t і x рівномірно обмежені в області D{t0tT, - ∞x∞} тоді головний член локальної похибки (8.17) можна оцінити, використовувавши формулу (8.18). Покладемо, що

, .

Позначимо через L постійну Ліпшиця функції f(t,x) по змінній x. Тоді для глобальної похибки (8.14) будемо мати наступну оцінку:

.

Підставляючи в останню формулу виразу для εn,εn - 1,…,ε0, отримаємо:

.

Оскільки ε0 = 0, та остання нерівність спрощується:

.

Використовуючи нерівність eLh>1 + Lh, перетворимо останній вираз до вигляду:

 (8.19)

і, застосовуючи його для (n + 1) = T, отримаємо остаточно:

. (8.20)

Таким чином, при припущеннях, прийнятих вище, метод Ейлера збігається на будь–якому кінцевому відрізку t[t0 ,T] і є методом першого порядку точності.

З іншого боку, підставляючи у вираз (8.19) M2 = 2ηn + 1/h2з виразу (8.18), знаходимо формулу, що пов’язує глобальну і локальну похибки:

, (8.21)

з якої при малих значеннях Lh(n + 1) знаходимо:

.

Формула (8.20) може в принципі застосовуватися для вибору кроку обчислень h при заданій похибці обчислень εn + 1 з урахуванням особливостей розв’язуваного диференціального рівняння (значень L і M2).

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

8.2.3. Поліпшення точності методу Ейлера застосуванням екстраполяції Рідчарсона

В інженерній практиці очевидна рекомендація зменшувати крок обчислень (h0) для підвищення точності розрахунків при застосуванні методу Ейлера не завжди прийнятна.Однак є можливість істотно підвищити точність розрахунків при збереженні досить великих значень кроку τ, що обмежується лише міркуваннями втрати стійкості обчислень, якщо скористатися формулою (1.4) екстраполяції Річардсона, приведеної раніше в параграфі 1.3 і вже використаної в главі 7. Нагадаємо, що якщо застосовується деякий алгоритм із відомою залежністю глобальної похибки від кроку обчислень, достатньо провести розв’язок задачі з різними кроками h і qh, щоб отримати потім точний розв’язок α0 за обчисленими наближеннями F(h) і F(qh):

ao = F() +  + 0(). (8.22)

Більше того, процедуру екстраполяції Річардсона можна узагальнити, як це було зроблено в рекурентних формулах Ромб ерга ( параграф 7.9) , на випадок врахування й інших членів розкладу в ряд глобальної похибки типу:

F(h) = ao + a1 p1 + a2 p2 + …,

де p1  p2  p3 

Тоді по обчислених наближеннях F(h), F(qh) , … , F(qkh), розглянутих як початковий цикл із k = 1, рекурсивно обчислюються наступні уточнення, обумовлені поточним показником ступеня pk :

Fk + 1( ) = Fk( ) + . (8.22,а)

Якщо Fk + 1 - Fk  , де — задана похибка обчислень, то Fk + 1 приймається як розв’язок. Відповідна схема обчислень набуває наступного вигляду:

Оскільки для методу Ейлера:

F(h) = ao + a1 + a2 2 + a3 3 + a4 4 + a5 5 + …, (8.23)

то при q = 2 перший рядок схеми обчислень стає рівним ∆/1, /3, /7, /13, /31,….

Приклад 8.1.

Розв’язати методом Ейлера з екстраполяцією Річардсона наступне рівняння:x = t2 - x2, t0 = 0, x0 = 1 Потрібно з кроком h = 0,1 знайти значення х(0,2).

Скориставшись формулою (8.12), записуємо un + 1 = un + hun,…, послідовно обчислюємо:

,
, .

Отже,u(0,2) = 0,82.

Тепер проведемо додатковий розрахунок із кроком h = 0.2:

, , чи .

Застосовуючи формулу екстраполяції Річардсона (8.22), знаходимо:

Приклад 8.2.

При розв’язку методом Ейлера лінійного диференціального рівняння x’ = - x при x(0) = 1

з кроком h = 0,25∙2 - i отримаємо для моменту часу t = 1 наступні чотири вихідні розв’язки Fi1 = X(1) (перший стовпець табл. 8.1.), який рекурсивно обробляємо, використовуючи формулу (8.22,а), і отримані дані розміщені в інших стовпцях цієї таблиці.

Таблиця 8.1. Результати обчислень за рекурентною формулою Ричардсона

X(1)

(/1)∙103

(/3)∙106

(/7)∙109

F01 = 0,316406

27203

F11 = 0,343609

F12 = 0,370812

12465

- 758

F21 = 0,356074

F22 = 0,368539

F23 = 0,367781

5981

- 168

12

F31 = 0,362055

F32 = 0,368036

F33 = 0,367868

F34 = 0,36788

Можна прийняти за розв’язок значення F34 = 0,36788 , тому що |F33 - F34| = 0,000012 при локальній помилці ≤1/210 - 6. Точне ж значення розв’язку рівняння дорівнює x(1) = e - 1 = 0.367879.

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

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

Базуються на використанні скороченого ряду Тейлора, але не вимагають обчислення похідних високих порядків. Для пояснення розкладемо точний розв’язок xi + 1 = x(ti + h) задачі (8.10) в околі вузла сітки ti в ряд Тейлора p–го порядку, отримаємо:

 8.24)

Основна ідея побудови явних методів Рунге–Кутта p–го порядку складається в заміні частини ряду 

функцією φP(xi,ti,h), яка не потребує визначення похідних від x’[t] = f(x,t) При цьому

 (8.25)

де С — константа, яка не залежить від кроку h. Узагальнений запис методу Рунге–Кутта наступний:

, (8.26)

де , (8.27)

bi — константи, а ki мають наступний рекурсивний вигляд

, (8.28)

Часто формули (8.27), (8.28) представляють наступною таблицею Бутчера:

Методи, визначені формулами (8.26), (8.27) і (8.28) називаються явними методами Рунге–Кутта через те, що в цих методах наступні наближені значення знаходять прямим розрахунком за формулами (8.27) і (8.28). Методи Рунге–Кутта різних порядків відрізняються кількістю членів в виразі (8.27). Розглянемо спочатку методи другого порядку і відповідні їм досить прості обчислювальні алгоритми.

8.3.1. Явний метод Рунге–Кутта другого порядку точності

Розглянемо метод другого порядку, для якого

(8.29)

де k1 = f(ti,ui); k2 = (ti + c2h,ui + a21k1h).

Коефіцієнти b1,b2,c1,a21 знаходять, прирівнюючи рівняння (8.29) і скорочений ряд Тейлора xi + 1 (8.24) з урахуванням члена з другої похідної, приведеного раніше в (8.18):

(8.30)

Для цього спочатку розкладають вираз для  у ряд Тейлора, як функцію двох змінних:

(8.31)

Підстановкою отриманого виразу в рівняння (8.26) з урахуванням значення для k1 знаходимо:

Після приведення подібних членів отримуємо:

(8.32)

Порівнюючи аналогічні члени в рівняннях (8.30) і (8.32), визначаємо, що для еквівалентності цих двох рівнянь необхідне виконання наступних умов:

(8.33)

Ці три умови містять чотири невідомі константи, тому для них немає єдиного розв’язку, однак, задаючи апріорне значення однієї з умов, ми можемо визначити три інші. Тому існує сукупність методів Рунге–Кутта другого порядку.

Припустимо, що ми вибираємо значення константи , тоді із системи (8.33):

(8.34)

Найбільш відомі наступні представники методів другого порядку:

8.3.1.1. Метод Хейна

Він реалізується при виборі b2 = 1/2. Тоді, відповідно до рівнянь (8.34) і (8.29):

(8.35)

де k1 = f(ti,ui); k2 = (ti + h,ui + hk1).

Неважко бачити, що k1 відповідає значенню похідної (нахилу дотичної до кривої розв’язку) на початку часового інтервалу, а k2 — значення цієї похідної наприкінці інтервалу.

8.3.1.2. Метод полігону

Реалізується при виборі b2 = 1. Тоді, згідно рівняння (8.29) b1 = 0, c2 = a21 = 1/2, і формула обчислень набуває вигляду:

(8.36)

де k1 = f(ti,ui); k2 = (ti + 1/2h,ui + 1/2hk1).

Цей же метод можна реалізувати трохи інакше, а саме: спочатку за формулою (8.37) обчислюється методом Ейлера проміжне значення:

 (8.37)

а потім, за формулою (8.38), обчислюють значення наближеного розв’язку в точці ti + 1:

(8.38)

Реалізація методу (8.36) у вигляді двох етапів (8.37) і (8.38) називається методом прогнозу–корекції, оскільки на першому етапі (8.36) наближене значення прогнозується з невисокою точністю O(h), а на другому етапі це прогнозоване значення виправляється, так що результуюча похибка має другий порядок по h.

8.3.1.3. Метод Ралстона–Рабіновича

Припускає вибір значення b2 = 2/3, при якому мінімізується локальна похибка методів другого порядку. При b1 = 1/3, c2 = a21 = 3/4:

(8.39)

де k1 = f(ti,ui); k2 = (ti + 3/4h,ui + 3/4hk1).

8.3.2. Явний метод Рунге–Кутта третього порядку точності

При розв’язку звичайних диференціальних рівнянь часто використовуються методи третього і четвертого порядку точності. Їх можна отримати тим же чином, яким були отримані формули для методу другого порядку в п. 8.3.1.

Для n = 3 у розкладу Тейлора (8.24) використовується додатковий член, що містить похідну від змінної x, або другу похідну від правої частини розв’язуваного диференціального рівняння x’ = f(x,t):

Замість рівнянь (8.30) отримуємо систему з шести рівнянь з вісьмома невідомими. Отже, значення двох з них повинні бути обрані апріорі для знаходження інших параметрів. Один з можливих варіантів методів третього порядку виглядає так:

(8.40)

де k1 = f(ti,ui); k2 = f(ti + 1/2h,ui + 1/2hk1); k3 = f(ti + h,ui - hk1 + 2hk2) 

Для будь–якого варіанта таких методів третього порядку локальна похибка O(h4) визначається першим неврахованим членом розкладу в ряд Тейлора, а глобальна похибка O(h3) визначається порядком методу. Як і у випадку методу Ейлера, можливе застосування екстраполяції Річардсона для поліпшення розв’язку. Однак тепер замість виразу (8.23) справедливим буде такий запис:

F(h) = ao + a22 + a33 + a44 + a55 + …,

У результаті при q = 2 перший рядок схеми обчислень по Річардсону містить виправлення
/3, /7, /15,…

8.3.3. Явний метод Рунге–Кутта четвертого порядку точності

Найбільше поширення отримали методи Рунге–Кутта четвертого порядку. Якщо застосувати підхід, описаний при виборі формул методу другого порядку, то зустрінемося з величезною кількістю можливих варіантів цих методів. Один з них, називаний класичним, описує наступну формулу наближення:

(8.41)

де

Можлива графічна інтерпретація формули (8.41), представлена на рис. 8.2.

Рис. 8.2. Графічна інтерпретація методу Рунге–Кутта 4–го порядку для одномірного випадку.

З рис.8.2 видно, як з початкової точки Pn виконується перший півкрок по дотичній з нахилом k1 і знаходиться точка а, у якій визначається нахил дотичної до кривої розв’язку k2. Потім, знову з початкової точки Pn, виконується наступний півкрок по прямій, але вже з нахилом k2. Знаходиться нова точка в і новий нахил дотичної до кривої розв’язку k3 Тепер вже по прямій з нахилом k3 з початкової точки робиться повний крок, що приводить до знаходження нової точки с і нового значення нахилу дотичної k4. Обчислені раніше нахили дотичної усереднюються і з отриманим значенням kp відбувається заключний повний лінійний крок, що приводить до знаходження наступної точки наближення Pn + 1.

8.3.4. Явний метод Рунге–Кутта високих порядків

Якщо потрібна більш висока точність результатів розрахунку, то можливе застосування методу Рунге–Кутта п’ятого порядку, запропонованого Бутчером:

(8.42)

де

Метод (8.42), у якому поліпшення точності досягається за рахунок істотного обсягу необхідних обчислень, застосовується все–таки рідше порівняно з класичним методом Рунге–Кутта четвертого порядку.

Приклад 8.3.

Розв’язати методом Рунге–Кутта наступне рівняння: x’ = x + t; x(0) = 1

Для кроку h = 0,2 знайти значення для х(0,2), скориставшись формулою методу Рунге–Кутта четвертого порядку (8.41):


;
;

При кроці  = 0,2 заповнюємо табл. 8.2.

Таблиця 8.2. Результати обчислень за методом Рунне–Кутта четвертого порядку

t

x

x’

k = x’

0

1

1

0,2

0,1

1,1

1,2

0,24

0,1

1,12

1,22

0,244

0,2

1,244

1,444

0,2888

Використовуючи дані останнього стовпця таблиці, знаходимо:

.

8.4. Апостеріорні оцінки похибки методів Рунге–Кутта

При виконанні практичних розрахунків необхідно забезпечити задану точність результатів обчислень. Апріорні оцінки точності для цього малокорисні. По–перше, головні члени похибки визначаються через похідні розв’язку, яка до початку розрахунку невідома (див. розділ 8.1.2, формула (8.17)). По–друге, апріорні оцінки звичайно є мажорантними і можуть у багато разів перевершувати фактичну помилку розрахунку. Тому основним практичним прийомом є апостеріорна оцінка точності. Щоб її знайти, розрахунок проводять на двох чи більшій кількості сіток, які згущуються, і застосовують правило Рунге, яке безпосередньо випливає з основного виразу екстраполяції Річардсона (8.22):

або в прийнятих у цій главі позначеннях при q = 1/2:

(8.43)

Нагадаємо, що локальна похибка методу розв’язку з порядком p визначається як:

або  (8.44)

Приклад 8.4.

Розглянемо розв’язок задачі Коші:

x’ = 2(t - 2) + t(t - 2)2 - tx, x(2) = 0 (8.45)

на відрізку t[2,4] за формулами Рунге–Кутта четвертого порядку за допомогою пакета Mathematica. Виберемо крок h = τ = 0,2, складемо програму розв’язку і визначення максимальної похибки отриманого розв’язку. Визначимо в Mathematica праву частину рівняння і обчислення членів формули Рунге–Кутта:

In = f[t_,x_]: = 2*(t-2) + t*(t-2)^2-t*x;
k1[
τ_]: = f[t,x];
k2[
τ_]: = f[t + 0.5*τ, x + 0.5*τ*k1[τ]];
k3[
τ_]: = f[t + 0.5*τ , x + 0.5*τ *k2[τ ]];
k4[
τ_]: = f[t + τ , x + τ *k3[τ]];

Введемо межі відрізку інтегрування, початкову умову і розв’яжемо задачу з кроком і 0,5τ:

In = m = 10; t0 = 2; x0 = 0; tm = 4; τ = (tm-t0)/m;
U = Array[u, {2, m + 1}, 0]; u[0,0] = t0; u[1,0] = x0;
U1 = Array[u1, {2, 2*m + 1}, 0]; u1[0,0] = t0; u1[1,0] = x0;
Ep = Array[ep, m + 1, 0]; ep[0] = 0;
J = 0; t = t0; x = x0; Do[x1 = x; t1 = t;
Do[dx1 =
τ/2*(k1[0.5*τ] + 2*k2[0.5*τ] + 2*k3[0.5*τ] + k4[0.5*τ])/6;
T = t + 0.5*
τ; x = x + dx1; j = j + 1; u1[0,j] = t; u1[1,j] = x, {k, 1, 2}];
x = x1; t = t1; dx =
τ*(k1[τ] + 2*k2[τ] + 2*k3[τ] + k4[τ])/6;
t = t +
τ; u[0,i] = t; u[1,i] = x + dx; ep[i] = Abs[u1[1, 2*i]-u[1, i]]/(2^5-1);
x = u1[1, 2*i], {i, 1, m}];

Побудуємо графіки сіткових функцій, які відповідають наближеним розв’язкам з кроками τ і 0,5τ. Щоб розмістити обидва графіки на одному малюнку, відкриємо графічний пакет і пакет, який дозволяє додати до малюнка ідентифікуючі написи:

In = <<Graphics‘MultipleListPlot‘
In = <<Graphics‘Legend

Приведемо графіки отриманих розв’язків

In = U = Array[u, {2, m + 1}, 0]; U1 = Array[u1, {2, 2*m + 1}, 0];
MultipleListPlot[Transpose[U], Transpose{U1],
PlotLegend->{“крок
τ”, “крок 0.5τ ”}, AxesLabel->{“t”, “U”}]

Рис. 8.3. Графік чисельного розв’язку методом Рунге–Кутта задачі Коші (8.45)

На графіку помітити різницю між значеннями наближених розв’язків у співпадаючих вузлах неможливо через малу похибку. Знайдемо за правилом Рунге оцінку помилки у вузлах основної сітки за формулою (8.43) і побудуємо її графік:

In = Ep = Array[ep, m + 1, 0]
ListPlot[Abs[ep], PlotStyle->PointSize[0.015], AxesLabel->{“n”, “Оцінка Рунге”}

Рис. 8.4. Графік оцінки по Рунге похибки чисельного розв’язку методом Рунге–Кутта задачі Коші (8.45).

Як випливає з отриманих оцінок, розв’язок при кроці τ = 0,5 отримано з високою точністю. Для розглянутого прикладу в Mathematica можна знайти аналітичний розв’язок задачі Коші наступним оператором:

In = Clear[t]
f[t_,x_]: = 2*(t-2) + t*(t-2)^2-t*x;
V = Dsolve[{v’[t] = = f[t, v[t]], v[2] = = 0}, v[t], t];
y[t_] = Simplify[V[[1, 1, 2]]];
Print[“y[t] = ”, y[t]]

y[t] = ( - 2 + t)2

Аналітичний розв’язок отриманий, тому можна знайти дійсну глобальну похибку розв’язку отриманого по методу Рунге–Кутта:

In = W1 = Array[w1, {2, 2*m + 1}, 0];
Do[w1[0, i] = u1[0, i]; w1[1, i] = u1[1, i]-y[u1[0, i]], {i, 0, 2*m}];
ListPlot[Transpose[W1], PlotStyle->PointSize[0.02],
AxesLabel->{“t”, “Глобальна похибка”}]

Рис. 8.5. Графік дійсної похибки чисельного розв’язку методом Рунге–Кутта задачі Коші (8.45)

Закінчуючи розгляд прикладу, приведемо чисельний розв’язок цієї ж задачі наступним оператором Mathematica. Щоб уникнути колізії використаних ідентифікаторів, введемо нові імена.

In = Vm = NDSolve[{v1’[z] = = f[z, v1[z]], v1[2] = = 0}, v1, {z, 2, 4}]
Out = {{v1→InterpolatingFunction[{{2., 4.}}, <>]}}

Знайдемо похибку розв’язку, отриманого в Mathematica із кроком, обраним за замовчуванням пакетом Mathematica.

In = y1[z_]: = v1[z]/.Vm[[1]]; Plot[y1[z]-y[z], {z, 2, 4},
AxesLabel->{“t”, “eps”}.

Рис. 8.6. Графік похибки чисельного розв’язку стандартним оператором Mathematica задачі Коші (8.45).

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

8.5. Методи з контролем похибки на кроці. Вкладені методи Рунге–Кутта-Фельберга

Для більш гнучкого керування вибором кроку інтегрування бажано мати можливість виконувати крок інтегрування і оцінювати похибку при меншій кількості значень правих частин, які обчислюються. Розглянутий вище в прикладі 8.4 метод покрокового контролю точності чисельного інтегрування диференціальних рівнянь і вибору при цьому розрахункового кроку на основі подвійного розрахунку за правилом Рунге приводить до значного збільшення обчислень. Дійсно, нехай для розв’язку задачі (8.10) використовується метод Рунге–Кутта четвертого порядку точності. Тоді виконання одного його кроку з контролем точності за правилом Рунге вимагає одинадцяти обчислень правої частини рівняння ( див. приклад 8.4). Для підвищення обчислювальної ефективності чисельного інтегрування диференціальних рівнянь і контролю покрокової похибки були розроблені методи, названі вкладеними методами Рунге–Кутта. Основна ідея складається у виконанні кроку від точки a до b методами порядку точності p і p + 1, для того щоб отримати оцінку похибки інтегрування методом порядку p. Аналогічно до формули (8.22) запишемо:

,

звідки

або в прийнятих у цій главі позначеннях при τ  1:

. (8.46)

Відповідно до виразів (8.26) і (8.27) формули методів Рунге–Кутта р–го і (р + 1)–го порядків можна записати відповідно як:

і

тому що константи bi і оцінки похідних ki для методів різних порядків, на жаль, найчастіше не збігаються. Похибка обчислень згідно (8.46) у загальному випадку оцінюється формулою:

(8.47)

тут деякі bi і ki, відсутні в методі р–го порядку, приймають нульові значення. Формулу похибки (8.46) можна спростити, якщо з можливої множини методів різних порядків підібрати такі, для яких оцінки похідних ki збігаються. Такі дослідження проведені Фельбергом для поєднання методів Рунге–Кутта різного порядку.

Наприклад, в формулах Фельберга другого порядку спільно використовуються формули Рунге–Кутта другого і третього порядків точності. Як основну формулу обрано формулу Рунге–Кутта другого порядку точності, яка визначається виразом (8.35), у якому b1 = b1 = 1/2 і k1 = f(ti,ui); k2 = f(ti + h,ui + hk1). Нагадаємо, що при виборі коефіцієнтів формули методу третього порядку формується система з шести рівнянь з вісьмома невідомими. Отже, значення двох з них повинні бути обрані апріорі для знаходження інших параметрів. Тепер вільні параметри обираються рівними чотирьом відповідним параметрам формули методу другого порядку, після чого визначаються інші чотири. У результаті отримуємо:

(8.48)

де k1 = f(ti,ui); k2 = f(ti + h,ui + hk1); k3 = f(ti + 1/2h,ui + 1/4k1h + 1/4k2h).

Фельберг вивів багато формул різних порядків. Нижче наведено формули найбільш розповсюдженого на сьогодні методу Фельберга 4(5), який поєднує методи Рунге–Кутта четвертого і п’ятого порядків і має похибку апроксимації порядку h4 і h5 відповідно:

, (8.49)

 (8.50)

де

(8.51)

Оцінка похибки, отримувана відніманням рівнянь (8.49) і (8.50), дорівнює:

(8.52)

Є й інший альтернативний вкладений метод Зонефельда, який поєднує також методи Рунге–Кутта четвертого і п’ятого порядків, але відрізняється більш простими коефіцієнтами:

(8.53)

(8.54)

де

(8.55)

Звичайна практика контролю похибки обчислень на кроці полягає в наступному. Користувач задає програмі дані про інтервал інтегрування Т і бажану загальну похибку розв’язку Е. Програма визначає припустиму похибку розв’язку на поточному кроці hn + i1 як частину загальної похибки Е, що припадає на цей крок,

(8.56)

і порівнює її з оцінкою похибки εn + 1, отриманої або з розв’язку з різними кроками (8.43) , або з розв’язку методами різних порядків (8.46). Якщо εn + 1<En + 1, то крок збільшується (наприклад, подвоюється) і, навпаки, при εn + 1>En + 1 крок зменшується (наприклад, удвічі).

Приклад 8.5.

Розв’яжемо задачу Коші, представлену виразом (8.45),

х’ = 2(t - 2) + t(t - 2)2tx; x(2) = 0,

на відрізку t[2,4] за формулами Фельберга 4(5) за допомогою пакета Mathematica.

Оберемо крок h = τ = 0,1, складемо програму розв’язку і визначення максимальної похибки отриманого розв’язку:

Введемо межі відрізка інтегрування, початкову умову і розв’яжемо задачу з кроком .

In = m = 20; t0 = 2; x0 = 0; tm = 4;  = (tm - t0)/m;
U = Array[u, {m + 1, 2}, 0]; u[0,0] = N[t0]; u[0, 1] = x0;
Ep = Array[ep, {m + 1, 2}, 0]; Array[el, {m + 1, 2}, 0];
Et = Array[et, {m + 1, 2}, 0];
ep[0, 0] = t0; ep[0, 1] = 0; et[0, 0] = t0; et[0, 1] = 0; el[0, 0] = t0; el[0, 1] = 0;
j = 0; t = t0; x = x0;

U = Array[u, {m + 1, 2}, 0]; u[0,0] = N[t0]; u[0, 1] = x0;
Ep = Array[ep, {m + 1, 2}, 0]; Array[el, {m + 1, 2}, 0];
Et = Array[et, {m + 1, 2}, 0];
ep[0, 0] = t0; ep[0, 1] = 0; et[0, 0] = t0; et[0, 1] = 0; el[0, 0] = t0; el[0, 1] = 0;
j = 0; t = t0; x = x0;

Нижче в табл..8.4 наведено перші 10 значень наближеного розв’язку , значення оцінки похибки, отриманої за методом Фельберга 4(5) і значення покрокової помилки – et, визначеної по точному аналітичному розв’язку задачі.

Таблиця 8.3. Результати обчислень за методом Фельберга 4(5)

Для того, щоб зробити свої методи оптимальними, Фельберг намагався вибрати коефіцієнти формули, яка визначає наближений розв’язок за формулою нижчого порядку, мінімізуючи локальну похибку їх результату. Як наслідок цього виходить, що оцінка похибки за формулою (8.51) може недооцінювати локальну похибку.

8.6. Чисельний розв’язок систем диференціальних рівнянь першого порядку

Виведені вище розрахункові формули Рунге–Кутта розв’язку задачі Коші для диференціального рівняння першого порядку легко переносяться на випадок чисельного розв’язку системи диференціальних рівнянь першого порядку. Для цього систему диференціальних рівнянь першого порядку достатньо записати у векторній формі.

Нехай потрібно знайти рішення задачі Коші для системи диференціальних рівнянь першого порядку, розв’язаних відносно похідних .

 (8.57)

Введемо наступні векторні позначення

 (8.58)

У використаних позначеннях задача Коші набуде вигляду :

 (8.59)

за такою ж формою, як розглянута вище задача (8.9). Формальна відмінність полягає в тому, що у відповідні співвідношення замість скалярних величин входять вектори. Тому до векторного диференціального рівняння (8.59) можна застосувати кожний з чисельних методів, які вивчалися в цій главі. Приведемо приклад векторних формул розв’язку системи диференціальних рівнянь методом Рунге–Кутта четвертого порядку для системи, записаної у вигляді (8.59):

 8.60)

Формули (8.60) за формою збігаються з формулами для одного диференціального рівняння першого порядку. Приведемо приклад їх використання.

Приклад 8.6.

Розв’яжемо методом Ейлера задачу Коші

y1 = - 0,5y1
y2 = 4 - 0,3y2 - 0,1y1

при початкових умовах y1(0) = 4; y2(0) = 6;

З міркувань точності обираємо h = 0,5. Знайдемо значення y1(2) і y2(2) , виконавши чотири кроки обчислень. Для першого кроку відповідно до методу Ейлера:

y1(0,5) = y1(0) + 0,5 [ - 0,5 y1(0) ] = 4 - 1 = 3;

y2(0,5) = y2(0) + 0,5[4 - 0,3 y2(0) - 0,1 y1(0)] = 6 + 0,9 = 6,9.

Результати розрахунку, отримані аналогічним чином на наступних ітераціях, зведені в табл.. 8.4.

Таблиця 8.4. Результати розв’язку системи диференціальних рівнянь методом Ейлера

x

y1

y2

1,0

2,25

7,715

1,5

1,6875

8,44525

2,0

1,265625

9,0940875

Проведемо повторне визначення значень y1(2) і y2(2) з кроком h = 1 і можливого застосування екстраполяції Річардсона для поліпшення отриманих результатів:

y1(1) = y1(0) + 1 [ - 0,5 y1(0) ] = 4 - 2 = 2;

y2(1) = y2(0) + 1[4 - 0,3 y2(0) - 0,1 y1(0)] = 6 + 1,8 = 7,8.

y1(2) = y1(1) + 1 [ - 0,5 y1(1) ] = 2 - 1 = 1;

y2(2) = y2(1) + 1[4 - 0,3 y2(1) - 0,1 y1(1)] = 7,8 + 1,8 = 9,26.

Після екстраполяції Річардсона отримуємо:

y1(2) = 1,265625 + ( 1,265625 - 1) = 1,531251

y2(2) = 9,0940875 + (9,0940875 - 9,26) = 8,928175.

Приклад 8.7.

Задана задача Коші

, (8.61)

а також початкові умови і крок інтегрування h = :

In = t0 = 0; x0 = {0.5, 1}; h = τ = 0.1;

Визначимо вектор–функцію правих частин рівнянь:

Задамо число кроків інтегрування і опишемо використані масиви вектора наближеного розв’язку U:

In = n = 10; Array[U, n + 1, 0];

Далі наведено програму розрахунку наближеного розв’язку з часовим кроком h = τ на n кроках.

Виведемо табл.. 8.5 перших десяти значень отриманого наближеного розв’язку задачі (8.64):

In = TA = Table[{t0 + (і-1)*τ , U[і-1][[1]], U[і-1][[2]]}, {i, 10}];
Print[“Таблиця
8.5”];
TableForm[TA, TableHeadings->{Automatic, {“t”, “
U1 ”, “U1 ”}}]

Таблиця 8.5. Результати розв’язку системи диференціальних рівнянь методом Рунне–Кутта

t

1

0

0,5

1

2

0,1

0,534027

1,16116

3

0,2

0,579456

1,34819

4

0,3

0,637922

1,56755

5

0,4

0,711651

1,82771

6

0,5

0,802855

2,13997

7

0,6

0,913105

2,51917

8

0,7

1,04313

2,98451

9

0,8

1,19313

3,5604

10

0,9

1,36313

4,27731

Одержимо розв’язок задачі (8.61) стандартним оператором Mathematica. (якщо стандартний оператор нам не відомий, тоді краще цей приклад опустити)

Отримано інтерполяційні функції, які визначають наближений розв’язок системи на заданому відрізку. За допомогою підстановки визначимо функції

Одержимо таблиці значень наближеного розв’язку системи і виведемо табл..8.7. для функції x1[t]:

In = τ = 0.1; Print[“Таблиця 8.6”];
Y = TableForm[Table[{i*
τ , x1[i*τ ][[1]]} {i, 0, 9}]
TableHeadings->{Automatic, {“t”, “
x1 ”}}]

Таблиця 8.6. Результати розв’язку системи диференціальних рівнянь стандартним оператором Mathematica

t

x1

1

0

0,5

2

0,1

0,534027

3

0,2

0,579456

4

0,3

0,63792

5

0,4

0,711647

6

0,5

0,802851

7

0,6

0,913102

8

0,7

1,04313

9

0,8

1,19313

10

0,9

1,36313

Наближені розв’язки збіглися з точністю до 10 - 5 

8.7. Стійкість методів Рунге–Кутта

Для практичних застосувань важлива оцінка поведінки розв’язку диференціального рівняння (системи рівнянь) у залежності від обраного значення кроку обчислень h. Для методів Рунге–Кутта високих порядків з їх достатнім запасом точності розв’язку бажана максимізація значень кроку  для зменшення загального числа оцінок правої частини розв’язуваного рівняння або загального часу обчислень. Однак явні методи розв’язку не дозволяють цього зробити, оскільки для них існує тверде обмеження виду [21]:

 (8.62)

де сконстанта (с < 3), λmax — максимальне власне значення матриці Якобі J = ∂f/∂xдля правої частини розв’язуваної системи диференціальних рівнянь (8.59). Якщо вирішується тільки одне рівняння типу (8.10), то λmaxвизначається значенням частинної похідної від правої частини цього рівняння, тобто

.

У випадках, коли для «жорстких» диференціальних рівнянь матриця Якобі має сильно рознесені власні значення (λmax/λmin = 106 - 109 ), обчислення повинні проводитися тільки з кроком, який визначається λmax за формулою (8.62), навіть на ділянці з повільною зміною розв’язку, обумовленою значенням λmin. В результаті загальне число кроків розв’язку неприпустимо зростає, що сильно обмежує застосування методів Рунге–Кутта для моделювання широкосмугових цифрових систем, що складають основу інструментарію інформаційних технологій.

Справедливість загальної формули (8.62) найпростіше довести, виходячи з умов стійкості розв’язку різницевих рівнянь, що і буде зроблено в главі 9. Тут же вона буде підтверджена прикладами розв’язку методами Рунге–Кутта модельного рівняння типу x’ = λx
(λ < 0), яке широко застосовується для аналізу стійкості чисельних методів розв’язку диференціальних рівнянь.

При цьому розв’язок модельного рівняння зводиться до виразу:

де u0 — початкове значення, un + 1— значення змінної на поточному кроці розв’язку, а стійкість розв’язку, коли un + 1→0 при n→∞, гарантується при виконанні умови

 (8.63)

Наприклад, для методу Ейлера (8.12) з урахуванням модельного рівняння:

звідки  для λ < 0

Приклад 8.8.

Дослідимо вплив вибору кроку h на стійкість розв’язку методом Ейлера задачі Коші x’(t) = - x(t) + 1; x(0) = 0,99.

Це рівняння має аналітичний розв’язок x(t) = 1 - et, t→∞, x(t)→1, що дозволяє оцінювати відмінність наближених розв’язків від точного. В табл. 8.7 приведені дані про п’ять ітерацій, проведених з різними кроками  = 0,5,2,01,5,0.

Таблиця 8.7. Вплив вибору кроку h на стійкість розв’язку методом Ейлера

hn

0

1

2

3

4

5

0,5

0,990

0,9950

0,9975

0,9987

0,8994

0,9997

2,01

0,990

1,0101

0,9898

1,0103

0,9897

1,0105

5,0

0,990

1,0100

0,8480

1,6400

- 1,5600

11,200

Припустимий крок h для розглянутої задачі відповідно до (8.62) дорівнює h < 2/λmax = 2 оскільки для досліджуваного рівняння λ = 1 по модулю. Перший рядок таблиці показує стійке наближення, тому що h = 0,5 < 2. Неприємності починаються з другого рядка таблиці (наближення набуває коливального характеру), коли h = 2,01 > 2. При значному перевищенні припустимого значення кроку (третій рядок таблиці) наближення розходиться.

Приклад 8.9.

Оцінимо припустиму величину кроку обчислень при розв’язку методом Ейлера задачі Коші

у’1 = - 0,5y1
у’2 = 4 - 0,3y2 - 0,1y1

при початкових умовах y1(0) = 4; y2(0) = 6 (задача розглядалася в прикладі 8.6).

Оскільки матриця Якобі для нашого прикладу дорівнює

,

то max| = 0,5 і максимальне значення кроку τ відповідно до умови (8.62) обмежене значенням
h < 4

Можна спростити оцінку максимально можливого кроку обчислень, якщо величину максимального модуля власних значень матриці Якобі оцінювати зверху величиною норми цієї матриці λmax≤||J||.

Для нашого прикладу ||J|| = |0,1 + 0,3| = |0,4|, що відповідно до умови (8.62) обмежує припустимий крок величиною h<5.

Аналіз стійкості методів Рунге–Кутта проводиться по аналогічній процедурі. Наприклад, для методу другого порядку (8.35) і модельного рівняння x’ = λx знаходимо:

де k1 = λui; k2 = λ(ui + hk1) = λui + 2.

Після відповідних підстановок і спрощень отримуємо:

звідки умова стійкості

Неважко показати, що для методу Рунге–Кутта четвертого порядку буде справедлива наступна умова вибору кроку для стійких наближень:

. (8.64)

Розв’язок нерівностей типу (8.64) для різних h представляє непросту задачу, тому розв’язок умови стійкості прийнято відображати в графічній формі. На рис. 8.7 приведені графіки функцій Mh) для явних методів Рунге–Кутта різних порядків. З графіків функцій Мh) для дійсних λ < 0 отримуємо оцінки:

, для до = 4 ; , для до = 3, , для до = 2,1. (8.60)

Рис. 8.7. Графіки функцій М(λh)

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

Оскільки λ може мати комплексні значення, то нерівність (8.64) прийнято відображати також у вигляді годографів у комплексній площині hλ. Побудуємо геометричне місце точок комплексної площини μ = , для яких |M| = 1, тобто покладемо M = e, де φ— аргумент M (Рис. 8.8). Наприклад, рівняння годографа для явного методу Рунге–Кутта четвертого порядку буде мати вигляд:

Побудова годографів стійкості по рівняннях такого вигляду можна виконати в Mathematica.

Рис. 8.8. Годографи стійкості явних методів Рунге–Кутта різних порядків

З практичної точки зору нас цікавить розташування годографів в лівій комплексній напівплощині, оскільки для фізично стійких об’єктів їх математичні моделі характеризуються від’ємними значеннями дійсних частин комплексної змінної hλ [Re() < 0]. Будь–які точки площини всередині годографа характеризують стійкі розв’язки, а поза годографом — нестійкі. Аналіз годографів стійкості методів Рунге–Кутта показує, що з ростом порядку методу область його стійкості розширюється, що є їх унікальною якістю.

ВИСНОВКИ:

1. Однокрокові методи Ейлера і Рунге–Кутта базуються на використанні скороченого ряду Тейлора для наближення функції, при цьому врахована кількість членів цього ряду k визначає їх порядок, а перший неврахований член ряду Тейлора — локальну похибку через величину кроку обчислень O(hk + 1).

2.  З метою усунення необхідності обчислення похідних високих порядків, що входять до членів ряду Тейлора, у методах Рунге–Кутта використовується спеціальна процедура апроксимуючих функцій, перші k членів розкладу яких у ряд Тейлора збігаються з k членами ряду Тейлора для вихідної функції.

3. Порядок методу k визначає залежність глобальної похибки розв’язку від кроку обчислень O(hk). Чим вище порядок методу, тим точніше розв’язок, але й більше обсяг обчислень.

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

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

6. Вкладені методи Рунге–Кутта–Фельберга (наприклад, РК4(5)) з автоматичною процедурою регулювання кроку обчислень є досить ефективними і розповсюдженими методами розв’язку диференціальних рівнянь, реалізованими у всіх відомих математичних пакетах.

7. Основним недоліком методів Ейлера і Рунге–Кутта, властивим, до речі, всім явним методам, є обмеження значення кроку обчислень h, перевищення якого приводить до розбіжності розв’язку. Це робить ці методи малопридатними для розв’язку «жорстких» задач з сильно рознесеними власними значеннями їх матриць Якобі, які характерні для моделювання широкосмугових цифрових інформаційних систем.

8. Для методів Рунге–Кутта з порядком k ≥ 3 кількість обчислень правої частини розв’язуваного диференціального рівняння на кожному кроці дорівнює k, що робить їх трудомісткими. Тому на практиці замість них використовують різницеві методи, описані в наступній главі, які використовують інформацію, отриману на попередніх кроках.

Вправи:

1. Для задачі Коші x’ = 1 + t2x2, x(0) = 0 знайти x(0,5) з чотирма значущими цифрами методом Ейлера з наступною екстраполяцією Річардсона наближень, отриманих з кроками h = 0,5, h = 0,25, h = 0,125.

2. При використанні методу Ейлера отримано наступні результати з різними кроками:

0,05

0,1

0,2

X

1,22726

1,22595

1,22345

Знайти поліпшений розв’язок за допомогою екстраполяції Річардсона.

3. Визначити діапазон можливих значень кроку обчислень h при використанні методу Ейлера для розв’язку задачі Коші: 

4. Для задачі Коші x’ = - 2t3 + 12t2 - 20t + 8,5; x(0) = 1 знайти x(1,5) при кроці h = 0,5:

а) методом Рунге–Кутта другого порядку (8.35);

б) методом полігону (8.36);

в) методом Ралстона (8.39);

г) методом Рунге–Кутта третього порядку (8.40);

д) методом Рунге–Кутта четвертого порядку (8.41).

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

5. Побудуйте вкладені формули Фельберга методу Рунге–Кутта 1(2) порядку, що задовольняють умові .

6. Протестуйте за зразком приклада 8.5 метод Ческіно:
,
де
розв’язавши цим методом задачу
Точний розв’язок тестової задачі знайдіть за допомогою системи
Mathematica. Побудуйте графіки похибок і оцінок помилки.

7. Складіть програму розв’язку задачі Коші для звичайного диференціального рівняння першого порядку вкладеним методом Мерсона з автоматичним вибором кроку, що забезпечує отримання наближеного розв’язку з заданою точністю:
,
де, , ,
,
Виконайте розв’язок наступної задачі Коші
по складеній програмі і порівняйте дійсну похибку з оцінкою похибки методу Мерсона.

8. Користуючись вкладеним методом Зонефельда (8.53)–(8.55), знайти наближений розв’язок системи рівнянь

Наведіть графіки оцінки похибки.

9. Побудуйте геометричне місце точок комплексної площини μ = hλ, для яких виконується умова стійкості явного методу Ейлера М = 1 + hλ.Відповіді, вказівки, розв’язки:

1. 0,5060 при точному значенні 0,5063.

2. 1,22861

3. 0  ≤1

4. а)3,43750; б)3,109375; в)3,2777344; г)3,21875; д)3,.21875
Точне значення 3,21875.

5. Розв’язок аналогічний п. 8.5.

6. Див. розв’язок прикладу п. 8.5.

7. Розв’язок аналогічний п. 8.5.

8. Аналогічно розв’язку прикладу 8.5 у п. 8.5.

9. Можна покласти М = q = , де μ – аргумент q. Отримаємо тоді для μ рівняння μ =  - 1 і побудуємо відповідну йому криву.

Рис.8.9. Годограф стійкості для метода Ейлера

PAGE  176

EMBED Equation.3  


 

А также другие работы, которые могут Вас заинтересовать

57829. Прийняття рішень в умовах невизначеності 527 KB
  Чим займається наука статистика Відповідь. Що таке математична статистика Відповідь. Що таке вибірка Навіщо використовують вибіркове спостереження Відповідь. Що таке ранжування Відповідь.
57830. Внутрішня будова стебла. Транспорт речовин по рослині 113 KB
  Мета уроку: Сформувати в учнів поняття про внутрішню будову стебла, зокрема кору, камбій, деревину, серцевину, їх будову і значення, про взаємозв’язок будови клітин і тканин стебла з їх фунціями...
57831. ВЛАСТИВОСТІ СТЕПЕНЯ ІЗ ЦІЛИМ ПОКАЗНИКОМ 2.17 MB
  Мета: працювати над засвоєнням учнями означення степеня з цілим показником та його властивостей; формувати вміння застосовувати означення і властивості степеня з цілим показником до обчислення значень виразів і перетворення виразів зі змінними...
57832. Перетворення виразів, що містять степінь з цілим показником 651 KB
  Число до якого ми підносимо основу степеня називається показником 3 Як називається число до якого ми підносимо основу степеня 3 Число яке ми підносимо до степеняназивається основою 4 Яке число ми отримуємо при піднесенні до...
57833. Степінь з цілим показником 79 KB
  Що називається степенем числа а з натуральним показником n 2. Парна степінь відємного числа завжди яке число 4. Непарна степінь відємного числа завжди яке число 5. Як називаються два дійсних числа сума яких дорівнює нулю 6.
57834. Штучні супутники Землі. Розвиток космонавтики 355.5 KB
  Розвивати творчу ініціативу і активність учнів, спостережливість, здатність логічно мислити, вміння робити висновки, узагальнення, розвивати пізнавальні здібності, навички працювати з науково-популярною літературою...
57835. Світ після Другої світової війни 756.5 KB
  При викладенні матеріалу застосовано дослідницькопошуковий метод і метод проектної роботи з історичним матеріалом розкрито творчий потенціал учнів. Продемонстровано різні методи...
57836. Міфи різних народів про творення світу. Слов’янські міфи про створення світу 95.5 KB
  Мета: поглибити знання про міфи народів світу; визначити схожі елементи у міфах про створення світу; визначити особливості виникнення цих міфів у Давньому Китаї та давніх скандинавів; вчити бачити звязок з уявленнями інших народів...
57837. Урок - гра: Подорож до 7 чудес світу 147 KB
  Мета: Закріпити знання учнів про алгоритм тотожних перетворень раціональних виразів способи перетворення відношення двох дробових виразів та про схеми застосування властивостей арифметичних дій під час перетворення раціональних виразів....