69315

БАГАТОКРОКОВІ МЕТОДИ РОЗВ’ЯЗУВАННЯ ДИФЕРЕНЦІЙНИХ РІВНЯНЬ

Лекция

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

В главі 8 було розглянуто однокрокові алгоритми обчислення наближеного розв’язку в точці tn + 1 з використанням інформації про розв’язувану задачу тільки на відрізку (tn,tn + 1) завдовжки в один крок. Логічно припустити, що можна підвищити точність методу...

Украинкский

2014-10-03

555 KB

7 чел.

ЛЕКЦІЯ 17. БАГАТОКРОКОВІ МЕТОДИ РОЗВ’ЯЗУВАННЯ ДИФЕРЕНЦІЙНИХ РІВНЯНЬ

Явні методи Адамса - Башфорта

В главі 8 було розглянуто однокрокові алгоритми обчислення наближеного розв’язку в точці tn + 1 з використанням інформації про розв’язувану задачу тільки на відрізку (tn,tn + 1) завдовжки в один крок. Логічно припустити, що можна підвищити точність методу, якщо використовувати інформацію про поведінку розв’язку в попередніх точках tn,tn + 1,tn - 2,…. Такі методи отримали назву багатокрокових. Загальна схема побудови багатокрокових методів, що використовуються для наближеного розв’язку задачі Коші

(9.1)

виглядає таким чином. Проінтегруємо рівняння (9.1) на відрізку [tn,tn + 1] і обчислимо значення xn + 1:

. (9.2)

Нехай відомо наближений розв’язок задачі (9.1) ui в m вузлах сітки tn - m + 1,…,tn - 1,tn. Отже, в цих точках відрізка відоме і значення Fi = f(ti,ui) правої частини диференціального рівняння (9.1) при i = n - m + 1...,n - 1,n, причому f(ti,u(ti)) буде вже функцією лише однієї змінної ti:f(ti,u(ti)) = F(ti). Замінимо в (9.2) функцію F(ti) інтерполяційним поліномом Hm(t), наприклад, поліномом Ньютона (6.29) з глави 6 і обчислимо наближене значення розв’язку un + 1:

. (9.3)

Інтегруючи, отримаємо різницеву схему для розв’язку диференціального рівняння. Порядок схеми визначається значенням залишкового члена інтерполяційного полінома. По суті, інтерполяційний багаточлен Hm(t) в формулі (9.3) використовується поза області інтерполяції, тобто в даному випадку це екстраполяційний поліном. Тому отриманий таким чином метод часто називають екстраполяційним методом Адамса—Башфорта. Отримаємо формули Адамса—Башфорта для різної кількості точок m. Задамо таблицю (tn - i,un - i),i = m - 1,
m - 2,…,0 заданих значень наближеного розв’язку, по яких можна обчислити значення
Fн - i = f(tn - i,un - i) і скласти таблицю (tn - i,Fn - i).

Візьмемо m = 2. Тоді таблиця матиме вигляд :

t

F

tn - 1

Fn - 1

tn

Fn

а відповідний поліном Ньютона запишеться в формі

Н1(t) = Fn + (t - tn)F(tn,tn - 1) = Fn + (t - tn)(Fn - Fn - 1)/h.

Підставимо Н1(t) в формулу (9.3):

.

Провівши інтегрування, отримаємо таку формулу для обчислення наближеного розв’язку:

. (9.4)

На відміну від методу Рунге—Кутта, який не вимагає попередньої побудови початку таблиці значень наближеного розв’язку, екстраполяційний метод Адамса при m = 2 вже передбачає попереднє знаходження величини u1. Визначимо похибку апроксимації отриманої формули як відхил при підстановці в різницеве рівняння значень точного розв’язку:

.

Розкладемо в ряди в околі точки tn функції, що входять у формулу для відхилу. Для точного розв’язку x(tn + h) отримаємо:

Оскільки fn - 1 = f(tn - 1,x(tn - 1)) = x(tn - h), то, розкладаючи в ряд похідну точного розв’язку, отримаємо:

Знайдемо відхил (нев’язка):

rn + 1 = - (xn + 1 - xn)/h–(xn - 1–3xn)/2 =

Таким чином, похибка формули (9.4) визначається останньою рівністю. Формули Адамса—Башфорта для інших значень числа кроків m можна отримати аналогічно. Приведемо ці формули для m = 3,4,5:

(9.5)

Трьохикроковий явний різницевий метод Адамса—Башфорта, очевидно, може застосовуватися, починаючи із значень n = 2. Похибка останньої формули визначається аналогічно попередньому випадку і дорівнює:

Звідси випливає, що головний член похибки на кроці для трьохточкової формули Адамса—Башфорта дорівнює:

.

Приведемо ще одну, більш точну, екстраполяційну формулу Адамса п'ятого порядку (m = 5):

(9.6)

Її похибка дорівнює:

Формула (9.6), очевидно, може застосовуватися, починаючи із значення n = 4. Як правило, відсутні значення функцій обчислюються в точках t1,t2,…,tm - 1 по методу Рунге—Кутти відповідного порядку, що є недоліком методу. Перевага методу Адамса полягає в тому, що на кожному кроці права частина диференціального рівняння обчислюється тільки один раз, а в методі Рунге—Кутти четвертого порядку точності на кожному кроці функція f(t,u) обчислюється чотири рази. Тому необхідно враховувати співвідношення кроків двох методів, що забезпечують задану точність.

9.2. Інтерполяційні методи Адамса—Мултона

Побудова інтерполяційних mкрокових методів Адамса - Мултона може бути здійснена абсолютно аналогічно випадку екстраполяційних формул. Для цього слід використовувати інтеграл (9.3), в якому інтерполяційний поліном Hm(t) побудовано по вузлах інтерполяції tn - m + 2,…,tn,tn + 1. Приведемо інтерполяційні формули Адамса різних порядків точності:

двокроковий неявний різницевий метод Адамса—Мултона:

(9.7).

трьохкроковий неявний різницевий метод Адамса—Мултона:

(9.8)

Визначення порядку апроксимації неявних методів Адамса—Мултона виконується аналогічно визначенню нев’язку методів Адамса—Башфорта. Похибка апроксимації трьохкрокового методу Адамса—Мултона дорівнює:

Помітимо, що у всіх приведених вище інтерполяційних формулах Адамса значення F1 + n = f(tn + 1,un + 1) невідоме, оскільки ще невідоме un + 1. Отже, інтерполяційні методи Адамса визначають un + 1 тільки неявно. Так, наприклад, інтерполяційна формула другого порядку точності:

(9.9)

справді є рівнянням щодо невідомого значення un + 1. Тому інтерполяційний метод Адамса—Мултона називають неявним.

На практиці, якщо рівняння (9.1) не є «жорстким» (див. параграф 10.1), рівняння (9.9) звичайно не розв’язують, а використовують спільно явну і неявну формули, що призводить до методу прогнозу і корекції. Одним з методів прогнозу і корекції, який широко використовується, є об'єднання методів Адамса четвертого порядку точності:


(9.10)
.

Звернемо увагу, що при такій послідовності обчислень цей метод є явним. Спочатку по екстраполяційній формулі Адамса—Башфорта обчислюється значення u(n + 1)(p), яке є "прогнозом" для un + 1. Потім u(n + 1)(p) використовується для обчислення наближеного значення F(n + 1)(p), яке, в свою чергу, використовується в інтерполяційній формулі Адамса—Мултона. Таким чином, формула Адамса—Мултона "коректує" наближення, отримане формулою Адамса—Башфорта. Слід відзначити, що при такій організації обчислень похибка, яка вноситься за рахунок неточності першої формули (9.10), не змінює порядок похибки обчислень по третій формулі (9.10).

Приклад 9.1.

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

(9.11)

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

Уточнене значення по інтерполяційній формулі Адамса обчислюють, використовуючи знайдене значення . Вираз називається формулою корекції:

Для того, щоб почати обчислення по формулі прогнозу, необхідно заздалегідь обчислити три перші значення ui, i = 1,2,3. Ці значення, звичайно, знаходять методом Рунге—Кутта. В розглянутому нижче прикладі ми зробимо так само. Спочатку проводиться розрахунок по методу Рунге—Кутта (див. приклад 8.14 з п.8.5) початкового відрізка таблиці значень (ti,ui), i = 1,2,3.). Потім виконується продовження обчислень таблиці значень (ti,ui), i = 4,5...,n, спочатку по формулі Адамса—Башфорта, а потім по формулі Адамса—Мултона. Одночасно виконується продовження розв’язку методом Рунге—Кутта. Графік отриманого розв’язку приведений на рис.9.1.

Рис.9.1. Графік чисельного розв’язку задачі Коші (9.17) комбінованим методом Адамса - Башфорта і АдамсаМултона

Точний розв’язок задачі Коші (9.11) легко знайти. Він має вигляд x = - 2 + t + 3exp( - t/2). Тому можна оцінити локальну і глобальну похибки наближеного розв’язку, отриманого методом Адамса—Башфорта і Адамса—Мултона. Значення локальних похибок на кожному кроці обчислень ec і накопиченої глобальної помилки розв’язку er приведені в табл.9.1.

Таблиця 9.1. Значення похибок обчислень до прикладу 9.1

t

ec

er

0

0

0

0

1

0,13

0

2,4 × 10 - 8

2

0,25

0

4,4 × 10 - 8

3

0,38

0

6,2 × 10 - 8

4

0,5

- 6,8 × 10 - 8

7,8 × 10 - 8

5

0,63

- 6,5 × 10 - 8

9,2 × 10 - 8

7

0,75

- 6,0 × 10 - 8

1,0 × 10 - 7

8

0,88

- 5,7 × 10 - 8

1,1 × 10 - 7

9

1,0

- 5,3 × 10 - 8

1,2 × 10 - 7

10

1,1

- 5,0 × 10 - 8

1,3 × 10 - 7

11

1,3

- 4,7 × 10 - 8

1,3 × 10 - 7

12

1,4

- 4,4 × 10 - 8

1,4 × 10 - 7

13

1,5

- 4,1 × 10 - 8

1,4 × 10 - 7

14

1,6

- 3,9 × 10 - 8

1,4 × 10 - 7

15

1,8

- 3,7 × 10 - 8

1,5 × 10 - 7

16

1,9

- 3,4 × 10 - 8

1,5 × 10 - 7

17

2,0

- 3,2 × 10 - 8

1,5 × 10 - 7

18

2,1

- 3, × 10 - 8

1,5 × 10 - 7

19

2,3

- 2,8 × 10 - 8

1,5 × 10 - 7

20

2,4

- 2,7 × 10 - 8

1,5 × 10 - 7

Дійсна накопичена помилка при розрахунку по методу Адамса і Рунге—Кутта практично однакові.

9.3. Лінійні багатокрокові різницеві методи

Розглянуті вище в п.п. 9.1 і 9.2 методи Адамса є окремими випадками лінійних багатокрокових різницевих методів, які в загальному випадку можна визначити такою формулою;

(9.12)

де ai,bi — числові коефіцієнти, не залежні від n, причому a0 0. Рівняння (9.12) слід розглядати як рекурентне співвідношення, яке виражає нове значення un + 1 через знайдені раніше значення un,un - 1,un - 2,…,un - m + 1.

Розрахунок починається з n = m - 1, тобто з рівняння

.

Звідси випливає, що для початку розрахунку необхідно задати m початкових значень u0,u1,u2,…,um - 1. Значення u0 визначається початковою умовою розв’язуваної задачі. Величини u0,u1,u2,…,um - 1 можна обчислити, наприклад, за допомогою методу Рунге—Кутта. Надалі припускатимемо, що початкові значення u0,u1,u2,…,um - 1 задані.

Метод (9.12) називається явним, якщо b0 = 0, і, отже, шукане значення un + 1, виражається явно через попередні значення un,un - 1,un - 2,…,un - m + 1. Інакше (тобто, коли b00) метод називається неявним. Тоді для знаходження un + 1 буде необхідно розв’язувати нелінійне рівняння (9.12).

Похибкою апроксимації на розв’язку диференціального рівняння (9.1) або нев’язкою різницевого методу (9.12) називається функція:

(9.13)

яку отримують в результаті підстановки точного розв’язку x(t) диференціальної задачі (9.1) в різницеве рівняння (9.12).

Оскільки коефіцієнти рівняння (9.12) визначені з точністю до множника, то для усунення невизначеності вважатимемо, що виконано умову:

(9.14)

що означає, що права частина різницевого рівняння (9.12) апроксимує праву частину відповідного диференціального рівняння (9.1).

З'ясуємо питання про похибку апроксимації при h0 залежно від вибору коефіцієнтів ai, bi, i = . Припускатимемо при цьому, що всі функції мають необхідну гладкість.

Розкладаючи функції xn - i = x(tn - i*h) в точці t = tn по формулі Тейлора, отримаємо:

і

.

Підставляючи ці розклади у вираз (9.13) для похибки апроксимації, матимемо:

Виділимо в кожній з сум доданки, відповідні j = 0, і згрупуємо їх: 

Скоротимо члени першої суми на h. Змінимо порядок підсумовування в подвійних сумах.:

В першій подвійній сумі виділимо доданки, відповідні індексу j = 1, і об'єднаємо їх з подібними членами:

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

(9.15)

Звідси видно, що похибка апроксимації має порядок m, якщо виконані умови:

(9.16)
(9.17)
. (9.18)

Разом з умовою нормування (9.14) рівняння (9.16), (9.17), (9.18) утворюють систему з m + 2 лінійних алгебраїчних рівнянь відносно 2(m + 1) невідомих a0,a1,…,am,b0,b1,…,bm.

Можна дещо спростити цю систему, використавши в рівнянні (9.17) умову нормування. Тоді отримаємо систему рівнянь:

, (9.19)

яка містить p рівнянь і 2m невідомих a0,a1,…,am,b0,b1,…,bm.

Коефіцієнти a0, b0 обчислюються за формулами:

(9.20)

Для того, щоб система (9.19) не була перевизначена, необхідне виконання умови p2m. Це означає, що порядок апроксимації лінійних m - крокових різницевих методів не може перевершувати 2m.

Отже, щонайвищий досяжний порядок апроксимації неявних m - крокових методів рівний 2m, а явних— 2m - 1.

Використовуючи рівняння (9.19) і (9.20), отримаємо формули для деяких лінійних неявних різницевих методів:

1. Лінійний двокроковий неявний метод четвертого порядку апроксимації:

(9.21)

2. Лінійний трьохкроковий неявний метод шостого порядку апроксимації:

(9.22)

3. Лінійний чотирикроковий неявний метод 8 - го порядку апроксимації:

(9.23)

Щоб отримати формули для явних лінійних багатокрокових різницевих методів, достатньо в рівнянні (9.18) покласти b0 = 0 і використовувати таку умову нормування:

(9.24)

Після спрощень рівняння (9.16), (9.17), (9.18) матимуть вигляд:

, (9.25)

Система (9.24), (9.25) містить p + 2 рівнянь і 2m + 1 невідомих. Для того, щоб система (9.24), (9.25) не була перевизначена, необхідне виконання умови p + 2< = 2m + 1. Ця вимога означає, що порядок апроксимації явних лінійних mкрокових різницевих методів не може перевершувати 2m - 1.

 Використовуючи рівняння (9.24), (9.25), отримаємо формули для деяких явних лінійних багатокрокових різницевих методів:

1. Лінійний двокроковий явний різницевий метод третього порядку апроксимації:

(9.26)

2. Лінійний трьохкроковий явний різницевий метод п'ятого порядку апроксимації:

(9.27)

3. Лінійний чотирикроковий явний різницевий метод 7 - го порядку апроксимації:

(9.28)

На закінчення відзначимо ще раз, що найвищий досяжний порядок апроксимації неявних m - крокових методів рівний 2m, а явних — 2m - 1.

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

Відзначимо, що, якщо в системі (9.12) відкинути останні k рівнянь, k = 1,2...,p - 1, то отримаємо умови, що забезпечують порядок апроксимації p - k.

Для методів Адамса—Мултона відповідно до формули (9.5) в загальній формулі (9.12) лінійних різницевих методів необхідно покласти:

.

Тоді умови (9.16), (9.17), (9.18) p - го порядку апроксимації набувають вигляду:

(9.29)

Звідси видно, що найвищий порядок апроксимації m - крокового методу Адамса—Мултона дорівнює m + 1, а найвищий порядок апроксимації методу Адамса—Башфорта, оскільки b0 = 0, дорівнює m:

(9.30)

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

9.4.1. Розв’язок систем рівнянь

Виведені в цьому розділі розрахункові формули розв’язку задачі Коші для диференціального рівняння першого порядку легко переносяться на чисельний розв’язок системи диференціальних рівнянь, записаної у векторній формі, аналогічно тому, як це було зроблено в розділі 8.7 для методу Рунге—Кутти. У векторних позначеннях, використаних в розділі 8.7, задача Коші має вигляд:

X' = F(t, X), X[t0] = X0. (9.31)

Тому до векторного рівняння (9.31) можна застосувати будь—який з чисельних методів, приведених в цьому розділі. Наведемо приклад векторних формул розв’язку системи диференціальних рівнянь явним методом Адамса. Введемо такі векторні позначення для сіткових функцій

(9.32)

Тоді явна векторна формула Адамса для системи диференціальних рівнянь матиме вигляд :

(9.33)

Формули (9.33) за формою співпадають з явною формулою Адамса для одного рівняння (9.6). Наведемо приклад їх застосування.

Приклад 9.2.

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

(9.34)

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

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

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

In = F[t_, x_]: = {E - (x[[1]]^2 + [[2]]^2) + 2t, 2X[[1]]2 + X[[2]]};

Обчислюватимемо наближений розв’язок по чотирикроковій явній формулі Адамса. Введемо вектор bj коефіцієнтів формули Адамса.

In = b = {55/24,59/24,37/24,9/24};

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

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

Початкові значення наближеного розв’язку в перших чотирьох точках сітки, необхідні для початку рахунку по формулі Адамса, задамо за результатами розв’язку цієї задачі методом Рунге—Кутта (див. розділ 8.7. Таблиця 8.6):

In = U[0] = {0.5,1.}; U[1] = {0.534027,1.16116}; U[3] = {0.637922,1.56755};

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

In =

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

In = TA = Table[{t0 + (i - 1)*τ,U[i - 1][[1]],U[i - 1][[2]]},{i,10}];
Print[“
Тиблиця 9.2.”]:

TableForm[TA,TableHeadings - > {Automatic,{“t”,”u1”,”u2”}}]

Отримані дані зведені у табл.9.2., які варто порівняти з табл. 6 параграфу 8.7.

Таблиця 9.2. Результати обчислення системи диференціальних рівнянь

T

u1

u2

1

0

0,5

1.

2

0,1

0,534207

1,16116

3

0,2

0,579456

1,34819

4

0,3

0,637922

1,56755

5

0,4

0,71181

1,82747

6

0,5

0,80326

2,1395

7

0,6

0,913684

2,51856

8

0,7

1,04369

2,98382

9

0,8

1,19354

3,55952

10

0,9

1,36344

4,27613

9.4.2. Розв’язок рівнянь вищих порядків

Одним з основних способів чисельного розв’язку задачі Коші для диференціальних рівнянь вищих порядків є зведення їх до відповідних задач для систем рівнянь першого порядку. Розглянемо задачу Коші для рівняння другого порядку:

(9.35)

Введемо нову змінну y = x, тоді y’ = x’’. Підставивши останні змінні в рівняння (9.35), отримаємо таку систему рівнянь, еквівалентну рівнянню (9.35):

(9.36)

Початкові умови для неї перепишуться у вигляді:

(9.37)

До задачі Коші (9.36), (9.37) можна застосувати будь—який чисельний метод з розглянутих в цій главі. Для пояснення наведемо приклад чисельного розв’язку задачі Коші для рівняння другого порядку.

Приклад 9.3.

Задана задача Коші для рівняння другого порядку:

.

Замінимо її еквівалентною задачею Коші для системи двох рівнянь першого порядку.

(9.38)

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

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

вектор правих частин рівнянь

вектор коефіцієнтів формул Мілна

Тоді векторні формули прогнозу і корекції для системи рівнянь (9.38) матимуть вигляд:

(9.39)

Початкові значення наближеного розв’язку в перших чотирьох точках сітки, необхідні для початку розрахунку по формулі прогнозу Мілна, знайдемо методом Рунге - Кутта (див. розділ 8.4.)

Задамо початкові умови:

In = t0 = 0; X0 = {4, - 4}; h = 0.1;

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

In = F[t_,x_]: = {X[[2]], - (9X[[1]] + 6X[[2]])};

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

In = n = 10; Array[U,n + 1,0]; Array[V,n + 1,0]; bp = {8/3, - 4/3,8/3};
bc = {1/3,4/3,1/3}

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

In = U[0] = x0; V[0] = x1; t = t0; D0[K1 = F[t,U[i]]; K2 = F[t + 0.5τ,U[i] + 0.5 τK1];
K3 = F[t + 0.5τ,U[i] + 0.5 τK2]; K4 = F[t + τ,U[i] + τK3];
U[i + 1] = U[i] + τ[K1 + 2K2 + 2K3 + K4]; V[i + 1] = U[i + 1]; t = t + c, {I,0,3}];

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

In = TR = Table[{t0 + (i - 1)*τ,U[i - 1][[1]],U[i - 1][[2]],{i,4}}];
Print["Таблица 9.3."];
TableForm[TR,TableHeadings - >{Automatic,{"t","u1","u2"}}]

Отримані дані представлені у табл.9.3.

Таблиця 9.3. Перші чотирі значения для застосування методу Мілна

t

u1

u2

1

0

4

- 4

2

0,1

3,55575

- 4,74055

3

0,2

3,07311

- 4,82859

4

0,3

2,6018

- 4,5526

Тепер можна продовжити розв’язок методом Мілна

In = Do[V[i + 1] = U[i - 3] + hhbp[[j]]F[t0 + h(i + 1 - j),U[i + 1 - j]];U[i + 1] = V[i + 1];
U[i + 1] = U[i - 1] + hbc[[j + 1]]F[t0 + h(i + 1 - j),U[i + 1 - j]]{i,3,n}];
TM = Table[{t0 + (i - 1)*h,U[i - 1][[1]],U[i - 1][[2]]}{i,n}];
Print["Таблиця 9.4."];
TableForm[TM,TableHeadings - >{Automatic{"t","u1","u2"}}]

Для даної задачі можна знайти точний аналітичний розв’язок:

In = EQ = Dsolve[{x'[z] = = y[z],y'[z] + 9x[z] + 6y[z] = = 0,x[0] = = 4,у[0] = = 4}
{
x[z],у[z]},z]
Out = =

Позначимо отримані аналітичні вирази відповідними функціями:

In = x1[z_] = EQ[[1,1,2]];x2[z_] = EQ[[1,2,2]];

Використовуючи отриманий точний розв’язок, можна знайти дійсну накопичену (глобальну) помилку розв’язку методом Мілна. Нижче наведено табл. 9.4, в якій наведено значення накопиченої помилки для кожної шуканої функції.

Таблиця 9.4. Значення накопиченої помилки для кожної шуканої функції

T

er1

er2

1

0

0

0

2

0,1

- 0,00018

0,00069

3

0,2

- 0,00024

0,00095

4

0,3

- 0,00024

0,00098

5

0,4

0,0004

- 0,0015

6

0,5

- 1,1Ч10 - 6

5,1Ч10 - 6

7

0,6

0,00056

- 0,0022

8

0,7

0,000019

- 0,00012

9

0,8

0,0005

- 0,002

10

0,9

- 0,000047

0,00011

11

1,0

0,00039

- 0,0016

Дійсна накопичена помилка при розрахунку методом Мілна er1 і er2 за абсолютним значенням більша від оцінки покрокової помилки, що є цілком природно.

9.5. Стійкість методів Адамса—Башфорта і Адамса—Мултона

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

(9.40)

де λ — довільне комплексне число. Властивості різних методів аналізують на прикладі модельного рівняння (9.40). Для того, щоб це рівняння дійсно моделювало початкову систему, необхідно розглядати його при всіх таких λ, які є власними значеннями матриці системи.

Стосовно тестової задачі (9.40), використаної раніше, багатокрокові формули Адамса - Башфорта типу (9.4), (9.5) і Адамса—Мултона (9.7),(9.7) та інших легко зводяться до різницевих рівнянь вигляду:

(9.41)

де ci = ai - hλbi.

Тому задача стійкості багатокрокових методів (9.12) зводиться до аналізу стійкості різницевого рівняння (9.41), методи для якого добре пророблені [8]. Стійкість різницевого лінійного однорідного рівняння з постійними коефіцієнтами типу

, (9.42)

якому задовольняє послідовність {xj}, де xj = ezj (z0, e0) визначається коренями його характеристичного рівняння:

, (9.43)

які повинні задовольняти умову

(9.44)

Тільки в цьому випадку розв’язок різницевого рівняння (9.42) рівний

при довільних константах ei наближається до нуля при n→∞. Цей же розв’язок при n→∞ зменшується до деякого кінцевого значення , якщо серед коренів (9.44) характеристичного рівняння є один корінь |Zi| = 1, що теж допустимо. Таким чином, для стійкості розв’язку задачі Коші необхідно, щоб корені характеристичного рівняння (9.42) лежали всередині одиничного кола {Z:|Zi|<1}, і ті з коренів, модулі яких рівні одиниці {Z:|Zi| = 1}, повинні бути простими.

Розглянемо умову стійкості методу Адамса другого порядку:

.

З урахуванням тестового рівняння x’ = λx знаходимо Fn = λun,Fn - 1 = λun - 1, тому відповідне різницеве рівняння (9.42) і його характеристичне рівняння (9.43) набувають вигляду:

,
,

звідки

або  (9.45)

Рис. 9.2. Область стійкості методу Адамса другого порядку

Отже, стійкість методу Адамса другого порядку вдвічі гірша, порівняно із стійкістю методу Ейлера першого порядку (для дійсних λ).

При достатньо великих за модулем негативних z поліном С(z) (9.42) має знак ( - 1)m. Тоді значення полінома С( - 1) повинне мати той же знак, інакше знайдеться корінь на інтервалі
( -
, - 1), який не задовольняє умову (9.44). Безпосереднє обчислення цього значення

призводить до нерівності

,(9.46)

яка накладає додаткові обмеження на коефіцієнти ci, необхідні для виконання умови стійкості (9.43).

Різницеве рівняння (9.41) з урахуванням значень коефіцієнтів ci для явних методів Адамса набуває вигляду:

,

якому відповідає характеристичне рівняння

(9.47)

У разі дійсних негативних λ для виконання умови (9.44) коефіцієнти характеристичного рівняння (9.47) повинні задовольняти нерівності (9.46), яка набуває вигляду:

(9.48)

Можна показати [26], що для явного методу Адамса m—го порядку справедливі такі співвідношення:

значення bi( - 1)m + iнегативні для всіх i;

значення  росте із збільшенням m.

Тому для виконання нерівності (9.44) необхідне виконання умови:

(9.49)

де bi — вагові коефіцієнти, що використовуються при знаходженні середнього значення похідної на декількох кроках обчислень. Наприклад, для явного методу Адамса третього порядку (9.5):


,

тому його умова стійкості :

(9.50)

Аналогічно можна показати, що для явного методу Адамса четвертого порядку

(9.51)

Це зменшення кроку для забезпечення стійкості розв’язку є своєрідною платою за зменшення об'єму обчислень. Таким чином, в методі Рунге—Кутта четвертого порядку крок може бути вибраний майже в 10 разів більшим, ніж у відповідному методі Адамса, що ставить під сумнів переваги останнього з точки зору врахування лише об'єму необхідних обчислень. З подальшим збільшенням порядку методу Адамса граничне значення продовжує зменшуватися (рис.9.3). Нагадаємо, що у випадку методів Рунге—Кутта, навпаки, збільшення порядку методу веде до розширення області стійкості методу.

Рис. 9.3. Області стійкості явних методів Адамса

Годографи стійкості, приведені на рис. 9.3, отримані на підставі виразу (9.47), який за допомогою конформного перетворення z = (U + 1)/(U - 1) зводиться до вигляду:

r(U) + hσ(U) = 0, (9.52)

де r(U) = і σ(U) = .

Про стійкість методу розв’язку певного порядку для вибраного кроку обчислень судять по положенню зображаючої точки на комплексній площині hλ відносно відповідної кривої годографа. Якщо точка потрапляє у внутрішню замкнуту область, то метод стійкий, якщо в зовнішню, — ні. Слід зазначити, що процедури прогнозу - корекції типу (9.10) не покращують властивостей стійкості використаних в них методів, якщо рівняння корекції не розв'язується як нелінійне алгебраїчне рівняння відносно un + 1.

Неявні багатокрокові методи Адамса - Мултона, наприклад, методи (9.7) і (9.8), які передбачають розв’язок нелінійних рівнянь, мають значно більшу стійкість, властиву в цілому неявним методам.

Рис. 9.4. Область стійкості неявного методу Адамса—Мултона 3—го порядку

Можна показати, що при hλ<0 для методу 3—го порядку умова стійкості h|λ|<6 (рис.9.4), а для методу четвертого порядку — h|λ|<6. Докладніше стійкість неявних методів розв’язку диференціальних рівнянь розглядається в наступному розділі.

ВИСНОВКИ:

1. Багатокрокові методи Адамса базуються на використанні інтерполяційних поліномів, що проходять через m попередніх точок розв’язку, для знаходження наступного наближення розв’язку диференціального рівняння.

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

3. Якщо розв'язано проблему початку розрахунку, то на кожному подальшому кроці необхідно лише один раз обчислити праву частину розв’язуваного рівняння, а не кілька разів, як в методах Рунге—Кутта.

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

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

6. Області стійкості багатокрокових методів Адамса значно менші порівняно з методами Рунге - Кутта тих самих порядків, при цьому із зростанням порядку методу його стійкість погіршується.

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

8. Вимушене зростання об'єму обчислень через обмеження часового кроку при розв’язку жорстких задач (задач з рознесеними швидкими і повільними рухами) змушує «забракувати» методи Адамса (так само як і методи Рунге—Кутта) для моделювання широкосмугових цифрових систем збору, обробки і передачі інформації й шукати інші альтернативні підходи для ефективного розв’язку таких задач.

Вправи:

9.1. Чи апроксимують наступні різницеві схеми рівняння (9.1):

9.1.1. ;

9.1.2. ;

9.1.3. ?

Якщо так, то, з яким порядком?

9.2. Для рівняння (9.1) побудуйте різницеву схему з найвищим порядком апроксимації вигляду:
.

9.3. Складіть програму обчислення коефіцієнтів для m—крокових різницевих методів Адамса—Башфорта і Адамса—Мултона, використовуючи формули (9.29) і (9.30) відповідно.

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

9.4 1. ;

9.4.2. ;

9.4.3. .

9.5. Покажить, що необхідною і достатньою умовою апроксимації рівняння (9.1) різницевими рівняннями (9.2) є виконання рівностей:

9.6. Отримайте такі різницеві формули Мілна, визначте порядок їх точності і знайдіть вираз головного члена локальної похибки:
явні формули (предиктор)

9.6.1. ;

9.6.2. ;

неявні формули (коректор)

9.6.3. ;

9.6.4. .

9.7. Складіть програму обчислення наближеного розв’язку задачі Коші таким методом Мілна:

9.8. Для задачі, розглянутої в прикладі 2 розділу 9.4, отримаєте розв’язок предиктор - коректорним трьохкроковим методом Мілна. Побудуйте графіки наближених рішень.

9.9. Різницеве рівняння xn + m + a1xn + m - 1 + … + amxn = 0 після підстановки вектора
Xj = (xj,xj + 1,…,xj + m - 1)t зводиться до вигляду: Xn + 1 = A Xn. Побудуйте відповідну матрицю А.

9.10. Методом Адамса—Башфорта 3—го порядку розв’яжить з кроком h = 0,1 задачу Коші dt/dx = x2 за таких початкових умов: x0 = 1,0000, x1 = 1,1111, x2 = 1,2500. Знайдіть значення x3 і x4.

9.11. Розвяжіть наступну задачу Коші методом Адамса—Башфорта четвертого порядку: dx/dt = - x/t в діапазоні від t = 4 до t = 5, використовуючи крок h = 0,5 і початкові значення x(2,5) = 1,2, x(3) = 1, x(3,5) = 0,857142857 і x(4) = 0,75. Використовуйте для оцінки похибки точні значення х(4,5) = 0,66666667 і х(5) = 0,6, знайдені аналітичним розв’язком.

Відповіді, вказівки, розв’язок

9.1. Так, ні, так. Скористатися формулами вправи 9.1.
Для визначення порядку апроксимації необхідно знайти
нев’язку різницевого методу (див формулу (9.13)). Використаємо Mathematica і знайдемо нев’язку рівняння вправи 9.1.2. Введемо таблиці коефіцієнтів ai і bi лівої і правої частин рівняння з прикладу 9.1.2.

In = Ta = {{3/2, - 1,1/2};Tb = {1,0,0};
Підставимо в різницеве рівняння точний розв’язок x(t) рівняння (9.1), знайдемо нев’язку при роз’язку диференціального рівняння (9.1):
;
і розкладемо її в ряд Тейлора в точці tn по h, зберігши перші три члени розкладу. Нижче іде оператор розкладу лівої частини рівняння
:

In = m = 2; L = Series[[Ta[[i + 2]]*x[tni*h]h{ф,0,m}]
Out = =

і оператор розкладу правої частини рівняння:

In = R = Series[[Tb[[i + 2]]*x[tn–i*h]/h{ф,0,m}]
Out = =

Тепер можемо знайти нев’язку:
In = rn = Simplify[ - L,R]; Print["rn = ",rn];

Таким чином, має місце апроксимація з другим порядком.

9.2. Див розв’язок вправи 9.1.2.
Введемо таблиці коефіцієнтів ai і bi лівої і правої частин рівняння:

In = m = 2; T1 = {{}; T2 = {b0,b1,b2};

 Знайдемо нев’язку при роз’язку диференціального рівняння (9.1) і розкладемо її в ряд Тейлора в точці tn по h, зберігши перші три члени розкладу. Нижче іде оператор розкладу лівої частини рівняння:

In = L = - Series[[T1[[i + 2]]*x[tni*h]/h{ф,0,m}]
Out =

 і оператор розкладу правої частини рівняння:

In = R = Series[[T2[[i + 2]]*x’[tn–i*h]{ф,0,m}]
Out =

 Тепер можемо знайти нев’язку 6

In = rn = Simplify[L + R]
Out = =

 Щоб отримати максимально можливий порядок апроксимації, необхідно bi, i = 0, 1, 2, обчислити як розв’язок наступної системи рівнянь:

In = eq1 = - 1 + b0 + b1 + b2 = = 0; eq2 = b0–b2 = = 0; eq3 = - 1 + 3b0 + 3b2 = = 0;
Solve[{eq1,eq2,eq3}{b0,b1,b2}]
Out = =

9.3. Див. вказівки п.9.5.2

9.4.1. Порядок апроксимації другий. Головний член рівний  9.4.2. Для отримання формул Ністрема третього порядку побудувати інтерполяційний багаточлен з вузлами tn - i, i = і проінтегрувати його в межах tn - 1, tn + 1.

In = T4 = Table[{tni*h,Fn - i}{i,0,2}];
Im3 = InterpolatingPolynomial[T4,t];
NIS3P = un - 1 + Simplify[Integrate[Im3{t,tnh,tn + h}]];
Print["un + 1",NIS3P]

 Головний член локальної похибки дорівнює  і визначається аналогічно методу, розглянутому в параграфі 9.2, при отриманні формули (9.9). Порядок апроксимації третій.

9.4.3.Порядок апроксимації четвертий. Головний член .

9.5. Випливає з формул (9.4) і (9.13).

9.6.1.Побудувати інтерполяційний багаточлен з вузлами

In = T = Table[{tn–i*h,Fn - i}{i, - 1,3}];

 Проінтегрувати інтерполяційний багаточлен від tn - 3 до tn + 1

In = In3 = InterpolatingPolynomial[Take,[T,4],t];
MILN3P = un - 3 + Simplify[Integrate[In3{t,tn–3h,tn + h}]];
Print["un + 1 = ",MILN3P]
.

 Порядок апроксимації четвертий. Головний член .

9.6.2. Порядок апроксимації четвертий. Головний член.

9.6.3. Проінтегрувати інтерполяційний поліном від tn - 1 до tn + 1:

In = MILN3C = un - 1 + Simplify[Integrate[In3{t,tn–ф,tn + ф}]];
Print["un + 1 = ",MILN3C]
.

 Порядок апроксимації четвертий. Головний член — .

9.7. Див. розділ 9.4, розв’язок прикладу 9.2.

9.8. Див. розділ 9.4, розв’язок прикладу 9.2.

9.9.

9.10. x3 = 1,4265 x4 = 1,6596.

PAGE  199


 

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

25623. Внезародышевые органы 43.5 KB
  Амнион временный орган обеспечивающий водную среду для развития зародыша. В эмбриогенезе человека он появляется на второй стадии гаструляции сначала как небольшой пузырек дном которого является первичная эктодерма эпибласт зародыша. Желточный мешок наиболее древний в эволюции внезародышевый орган возникший как орган депонирующий питательные вещества желток необходимые для развития зародыша.
25624. Зрительная сенсорная система. Орган зрения 61 KB
  В передней части глаза склера переходит в покрытую многослойным плоским эпителием прозрачную роговицу. Наружная фиброзная оболочка глазного яблока к которой прикрепляются наружные мышцы глаза обеспечивает защитную функцию. Внутренняя чувствительная оболочка глаза сетчатка сенсорная рецепторная часть зрительного анализатора в которой происходят под воздействием света фотохимические превращения зрительных пигментов фототрансдукция изменение биоэлектрической активности нейронов и передача информации о внешнем мире в подкорковые...
25626. Гистогенез и органогенез на 2 и 3 неделе развития 26.5 KB
  Коммутирование ограничение возможных путей развития клеток. Оно совершается последовательно: сначала преобразуются крупные участки генома детерминирующие наиболее общие свойства клеток а позднее более частные свойства. Дифференцировка это изменения в структуре клеток связанные с их функциональной специализацией обусловленные активностью определенных генов. В развивающемся организме дифференцировка сопровождается определенной организацией или размещением специализирующихся клеток что выражается в установлении определенного плана...
25627. Гистогенез и органогенез 22 KB
  4 неделя Углубление желточной складки образование желточного стебля и приподнятие зародыша в полости амниона. Замыкание нервной трубки и формирование переднего невропора к 25 сут и заднего невропора к 27 сут образование нервных ганглиев; закладка легкого желудка печени поджелудочной железы эндокринных желез аденогипофиза щитовидной и околощитовидных желез. Образование ушной и хрусталиковой плакод первичной почки мезонефроса. Образование зачатков верхних и нижних конечностей 4 пар жаберных дуг.
25628. Гладкие мышечные ткани 29.5 KB
  Стволовые клетки и клеткипредшественники в гладкой мышечной ткани на этапах эмбрионального развития пока точно не отождествлены. Поверх чехликов из базальной мембраны между миоцитами проходят эластические и ретикулярные волокна объединяющие клетки в единый тканевой комплекс. Ретикулярные волокна проникают в щели на концах миоцитов закрепляются там и передают усилие сокращения клетки всему их объединению. Поэтому после поступления нервного импульса медиатор распространяется диффузно возбуждая сразу многие клетки.
25630. Дифференцировка первичной эктодермы 39 KB
  Меньшая часть эктодермы расположенная над хордой нейроэктодерма дает начало дифференцировке нервной трубки и ганглиозной пластинки. Нейруляция процесс образования нервной трубки протекает по времени неодинаково в различных частях зародыша. Замыкание нервной трубки начинается в шейном отделе а затем распространяется кзади и несколько замедленнее в краниальном направлении где формируются мозговые пузыри. Из нервной трубки образуются нейроциты и нейроглия головного и спинного мозга сетчатки глаза и органа обоняния.