Міністерство освіти і науки України

Прикарпатський національний університет імені Василя Стефаника

Факультет математики та інформатики
кафедра диференціальних рівнянь і прикладної математики
Курсова робота на тему:

Деякі скінченно-різнецеві методи розв’язування


звичайних диференціальних рівнянь

Виконав:

студент групи ПМ-41

Васьків Святослав
Перевірив:

науковий керівник:

Василишин П.Б.
Івано-Франківськ 2010


План
Вступ.

1. Чисельна ітерація рівнянь Ньютона

2. Алгоритм Бімана і Шофілда

3. Метод Рунге-Кутта

a) Метод Рунге — Кутта 4-го порядку

б) Неявні схеми Рунге-Кутта

в) Неявні інтерполяційні схеми

г) Програма Рунге-Кутта на мові С#

д)  Програма Beeman

4. Метод Адамса

5.Метод  Крилова

6. Метод Чаплигіна

Висновок

Список використаної літератури


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


1. Чисельна ітерація рівнянь Ньютона

Для спрощення запису розглянемо одновимірний рух частини і запишемо рівняння Ньютона у виді:
                           (1)

                         (2)
Метою усіх скінцеворізнецевих методів являється знаходження значень x _n+1 і  v _n+1(точка в "фазовому просторі") у момент часу t_n+1=t_n+∆t Нам вже відомо, що величину кроку ∆t  потрібно вибирати таким чином, щоб  метод інтегрування породжував приймати однакове рішення. Один із способів перевірки стійкості методу полягає в контролі величини повною енергії і забезпеченні того, щоби вона не відхилялася від початкового значення у разі, коли повна енергія зберігалась. Досить велике значення кроку приводить до не збереження повної енергії і до різних розв’язків для х_n+1 i v_n+1, тобто до таких розв’язків, які все більше відхиляються з потоком часу від  істинного розв’язку.
Суть багатьох алгоритмів, можна зрозуміти, розкладаючи

 v_n+1 ≡ v(t_n+∆t) i

x_n+1 ≡ x(t_n+∆t ) в ряд Тейлора. Запишемо



                    (3)

і

            (4)
Добре відомий метод Ейлера еквівалентний збереженню в формулі (3) членів
                      (5)

і

               (6)
Оскільки ми утримали у формулах (5-6) члени порядку ∆t, то «локальна» погрішність (погрішність на кроці) складає величину O(∆t)²

Оскільки ми від кроку до кроку погрішності накопичуються, позтому можна припускати, що «глобальна» погрішність, що є сумарною погрішністю за розглядом проміжок часу, буде величиной O(∆t). Ось ця оцінка погрішності цілком правдоподібна, оскільки число кроків, на яке розбивається часовий інтервал, пропорційна 1/∆t. Звідси випливає, порядок глобальної погрішності збільшується в ∆t разів по відношенню до локальної погрішності. Оскільки прийнято "говорити, що метод має n-й порядок аппроксимації, якщо ця локальна погрішність рівна О((∆t)ⁿ⁺¹), то метод Ейлера відноситься до методів першого порядку.

Метод Ейлера являється асиметричним, оскільки він просуває вирішення на один часовий  крок ∆t, а використовує при цьому інформацію про похідну тільки в початковій точці інтервалу. Ми вже переконалися в тому, що точність методу Ейлера обмежується і частенько породжуване його рішення нестійке. На щастя, як правило, немає необхідності використовувати більш складні алгоритми. Наприклад, швидше ми  знайшли, що проста модифікація методу (5-6), запропонована Кромером і іншими авторами, породжує стійке рішення для коливальних систем. Для закріплення повторимо алгоритм Ейлера-Кромера або наближення  по «останній точці»:



                  (7)

і

          (8)
Мабуть найбільш очевидний шлях удосконалення методу Ейлера, полягає в використанні для обрахування нового значення координати середини на відрізку швидкості. Відповідний метод середньої точки можна записати в вигляді:
                  (9)

і

     (10)
Замітимо, що якщо підставити вираз (9) для v_n+1 в (10) то отримаємо
    (11)
Звідси випливає, що метод середньої точки другого порядку точності по координаті і першого порядку по швидкості. Хоча наближення по середній точці дає точні результати для випадку постійного прискорення, взагалом він не приводить до значного покращення результату, ніж метод Ейлера. На справді оці два методи однаково погані, оскільки з кожним кроком погрішність збільшується.

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

і

           (13)
Замітио, що метод напів кроку не являє собою «самостартючим», так чи однакще фомули не дозволять порахвати v_1/2. Цю незручність модна здолати поклавши
                       (14)
Оскільки формули (12-14) можна повторювати до безмежності, то метод напів кроку отримав широке розповсюдження в навчальній літературі.

Один із найбільш відомих алгоритмів вищого порядку привласнюєтья Верле. Запишемо в ряд Тейлора для х_n-1



           (15)

Якщо скласти формули інтегрування вперед і назад (вирази (3) і (14) відповідно) то отримаємо
                (16)

або

                            (17 )
Аналогічно розв'язання розкладу в ряд Тейлора для x_n
+1
i

x_n
-1  дає
                                 (18)
Підмітимо, що звязок з алгоритом Верле (18) велика погрішність має третій порядок для координати і другий порядок для швидкості. Проте швидкість не бере участі в інтегруванні рівнянь руху. В літературі по чисельному аналізу алгоритм Верле називається «неявна симетричність різновидної схеми».

Менш відомим, про те математично еквівалентної версії алгоритму Верле являє собою схема
               (19)

і

             (20)


Видно, що схема (19-20), називається швидкісною формулою алгоритму Верле, являється самостартуючою і не приводить до накопичення погрішностей округлення. Формули (19-20) можна вивести із формул (16-19) наступним чином.. Спочатку додамо і віднімемо із рівнянь (16-19) по (1/2)х_n+1 і запишемо:


        (21)
Тут ви використаємо вираз (18). Із (17) знайдемо а_n для метода Верле:



                         (22)
Легко помітити, що підстановка (22)в вираз (21) приводить до (19). В томуж дусі перепишемо (17) для v_n+1



                              (23)
Тепер перепишемо формулу (17) для х_n+2 і підтавимо її в отриманий результат формули (23). Отримуємо:



              (24)

Тоді використовуємо вираз (17) для x_n+1, повторимо цю процедуру і поставимо x_n+1 в (24); після не важких перестановок отримаємо потрібний результат (20)

2. Алгоритм Бімана і Шофілда

Інший корисний алгоритм, в якому немає нагромаження погрішностей округлення, як в алгоритмі Верле, належить Біміану і Шофілду. Запишемо алгоритм Бімана в наступному виді:
               (25а)

і

                (25б)
Підмітимо, що точність розразунку траєктторії по схемі (25) не вища, ніж в алгоритмі Верле. Її перевага заключається в тому, що просто вона краще зберігає енергію. Однак алгоритм Бімана не самостартуючий. Алгорим Бімана і алгоритм Верле в швидкісній вормулі викоритані в програмі BEEMAN

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





                         (26а)
Передбучуване значення коодринати дозволить оприділити прискорення

Тоді, використовуючи , отримаємо скоректироване значення v_n+1 і x_n+1
коректор:  


                                (26б)
Скоректироване значення x_n+1 використовується для визначення нового передбачуваного значення а_n+1 і, значить, нових передбачених значень v_n+1 i x_n+r Ця процедура повторяється до тих пір, доки передбачення і скоррективонане значення x_n+1 відрізняються менше ніж на задану величину! Даний метод можна розробити на схемі більш високого порядку, які зв’язуються між собою не тільки  x_{n+1 ,}  x_n і v_n , але і так само значеннями v_n-1 і v_n-2. Замітимо, що метод предиктора-коректора не являється самостартуючим.
3. Метод Рунге-Кутта
Для пояснення методу Рунге-Кутта подивимось спочатку розв’язок диференціального рівняння першого порядку



                                     (27)
Метод Рунге-Кутти другого порядку для розв’язку рівняння (27) модна, використовуючи стандартні значення, записати наступним чином:

                         (28)
Сенс формул (28) полягає у наступному: В методі Ейлера допускається, що для екстаполяції в наступну точку модна використовувати нахил кривої f(x_n,y_n)в точці (x_n,y_n) так чи однакше y_n+1=y_n+f(x_n,y_n)*∆x. Однак можна повисити почність оцінки нахилу, якщо методом Ейлера повести екстраполяцію в середню точку відрізку, а потім використати центральну похідну на всьому відрізку. Звідси оцінка нахилу в методі Рунге-Котти рівна
 де
Застосування методу Рунге-Кутти до рівнянь руху Ньютона дає
                             (29)
Оскільки методи Руиге-Кутти є такими, що самостартуючими, то їх часто використовують для вираховання декількох перших кроків для несамостартуючих алгоритмів.


4. Метод Рунге — Кутта 4-го порядку

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

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



        
Тоді значення невідомої функції в точці  x_n+1 обчислюється відносно значення в попередній точці  x_n по такій формулі:




де h— крок інтегрування, а коефіцієнти  k _n  розраховуються наступним чином:

Це метод 4-го порядку, тобто похибка на кожному кроці становить O(h⁵), а сумарна похибка на кінцевому інтервалі інтегрування є величиною O(h⁴) .


Прямі методи Рунге — Кутта

Група прямих методів Рунге — Кутта є узагальненням методу Рунге — Кутти 4-го порядку. Воно задається формулами




де




Конкретний метод визначається числом s і коефіцієнтами b_i,a_ij i c_i . Ці коефіцієнти часто впорядковують в таблицю

            0

         c₂       a₂₁

         c₃       a₃₁     a₃₂


            ∙       ∙       ∙     ∙

         ∙       ∙       ∙          ∙

         ∙       ∙       ∙               ∙

         c_s       a_s1     a_s2 ∙   ∙   ∙     a_{s,s − 1}


        

                   b₁      b₂                b_{s − 1}  b_s
Для коефіцієнтів методу Рунге — Кутта мають справджуватись умови 
 

Якщо ми хочемо, щоб метод мав порядок p, то варто так само забезпечити умову   — наближення, отримане по методу Рунге — Кутти. Після багаторазового диференціювання ця умова перетвориться в систему поліноміальних рівнянь на коефіцієнти методу.
5. Неявні схеми Рунге-Кутта
Викладений тут алгоритм є неявна схема Рунге-Кутта четвертого порядку. У неї, як завжди, вбудована оцінка погрішності, яка дорівнює різниці відповідей четвертого і третього порядків точності, вычисленних по використаних точках. Іншими словами, в порівнянні з рекомендованим вище явним методом Рунге-Кутта п'ятого порядку, усі порядки в точності на одиницю менше. Нагадаємо, що в методі Рунге-Кутта п’ятого порядку використовується шість точок. У неявній схемі точок буде значно менше, оскільки для неявних схем зв'язок порядку точності з кількістю використаних точок не така, як для явних схем. Наприклад, як ми вже бачили, одноточкова неявна схема може мати другий порядок точності. У нашій неявній схемі четвертого порядку ми обійдемося всього-навсього трьома точками:









На перший погляд, слід спочатку вичислити матрицю , а потім застосовувати її потрібну кількість разів до усім f
_j. Проте, як ми знаєм (див. частину I, розділ "Системи лінійних рівнянь"), це не є правильний спосіб рішення СЛР "про запас". Правильний спосіб рішення СЛР з цією матрицею  раз і назавжди (для довільної правої частини) - це LU - розкладання матриці . При цьому матриця   розкладається на ліву нижню трикутну матрицю L і праву верхню трикутну матрицю U. Після цього рішення будь-хто СЛР з матрицею  будується за допомогою прямої підстановки (для матриці L) і потім зворотної підстановки (для матриці U). Кожна з підстановок вимагає N² операцій. У нашому випадку має сенс скомбінувати праві частини усіх СЛР так, щоб кількість застосувань матриці  (точніше, кількість СЛР, які потрібно вирішити) була мінімальною. Для цього досить ввести проміжні змінні . В результаті ми отримаємо наступний рецепт для одного кроку за часом t→t + h.

Алгоритм:

Для поточної точки обчислюваний  

Крім того, вираховуєм

Виконуємо LU - розкладання матриці А. (Це робиться один раз і назавжди для цього кроку за часом t→t + h).

За допомогою LU-розкладання обчислюємо . Обчислюваний . Обчислювано . З допомогою LU  - розкладання обчислюваний  

Обчислюємо  

Обчислюваний .

За допомогою LU розкладання обчислюємо



За допомогою LU розкладання обчислюваний



Обчислюємо значення :



Обчислюємо погрішність      

- Перетворюємо вектор Е на одне число є, що характеризує погрішність на цьому кроці.


Ця схема призводить до досить забавної поведінки "швидких" змінних. Поведінка виявляється абсолютно однаковою як у разі релаксаційного рівняння , так і у разі осциляторного рівняння . Саме, при Ch>>1 "швидкі" змінні експоненциально вибуваєть, причому швидкість вибування абсолютно не залежить від величини C:  x(t+h)=x(t)/3. Отже стійкість гарантована, але сподіватися на правильний опис еволюції "швидких" змінних не доводиться.

Повний алгоритм рішення системи ОДУ з адаптивною зміною кроку будується на основі цього рецепту таким самим чином, як будувався повний алгоритм в звичайному (явному) методі Рунге-Кутта п'ятого по-рядка. Єдина модифікація полягає в тому, що міри "1/4" і "1 /5", відповідна п'ятому порядку точність явної схеми, належить замінити на мірі "1/3" і "1/4", відповідна четвертому порядку точність неявної схеми.

Необхідно також нагадати, що для "жорсткої" системи ОДР перетворити вектор Е в одне число 6, що оцінює погрішність, набагато важче, чим для звичайної системи ОДР.
6. Неявні інтерполяційні схеми
Алгоритм неявних інтерполяційних методів дослівно повторює алгоритм звичайних (явних) інтерполяційних методів. Єдина відмінність полягає в тому, що ми не використовуємо звичайну (явну) схему Рунге-Кутта другого порядку:

Ми заміняємо її на неявну схему Рунге-Кутта другого порядку:





Нагадаємо, що для цієї схеми можна використовувати два різних рецепта:

•                      шукати точне рішення цієї системи нелінійних рівнянь методом Ньютона;

•                      обмежитися лінеаризованим рівнянням, тобто зробити тільки перший крок методу Ньютона :




У будь-якому випадку доведеться вирішувати СЛР. В більшості випадків можна обійтись другим варіантом, але в дуже нестійких завданнях доведеться вдатися до першого (до речі, при цьому не доведеться сильно міняти текст програми).

Нагадаємо, що існують варіанти інтерполяційного методу з фіксованим кроком Н і змінним порядком m, з фіксованому порядку те і змінним кроком Н і, нарешті, із змінними m і Н. Тут можна привечти тільки простий алгоритм з фіксованим H.

Алгоритм:

Вибираємо постійний крок H за часом (як правило слід брати H близько 0.2). Рухаємося по траєкторії з цим кроком. Кожен окремий крок реалізується так:

Для k = 1,2,.. виконуємо наступне:

Для чергового N_k з набору {4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128} обчислюємо крок h_k = H/N_k.

За допомогою цього кроку h_k знаходимо





тут t_n = t+nh_k, n = 0.. N_k-1; при цьому . Зрозуміло, тут не слід обчислювати зворотну матрицю, потрібно просто вирішувати СЛР.

Для кожного х_і по k  точок:




виконуємо поліноміальну інтерполяцію взалежності Y(X) в точку X = 0. Результат інтерполяції, отриманий на k -му етапі, означаємо . Оцінюємо погрішність
 як
Перетворюємо вектор на число характеризуючи погрішність. Якщо  менше замовленого  те вважаємо  і завершуємо даний крок t
→
t
+
H
.

Перша неприємна відмінність від явного методу полягає в тому, що кожен маленький крок на h_k (великий крок H складається з N_k штук таких кроків) вимагає рішення СЛР.

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

На перший погляд, цей алгоритм приведе не просто до багатьох помилок в описі динаміки "швидких" змінних, а до помилок якісних. Дійсно, якщо рівняння для "швидкої" змінної має вигляд , то при Ch
>>1 неявна схема Рунге-Кутта другого порядку замість експоненціального вибуває дає x
(
t
+
h
)≈-
x
(
t
), тобто |x
(
t
+
H
)|≈|
x
(
t
)|. На щастя, наступна інтерполяція виправляє ситуацію. Інтерполяційний алгоритм "помічає", що у міру зменшення h_k  абсолютна величина          x
^(
k
)
(
t
+
H
) зменшується. Причому це зменшення тим помітніше, чим менше h_k - В результаті остаточна відповідь  буде помітна менше за абсолютною величиною, чим x(t).

№ наближення	і	x	y
І	0 1	x0 x1
ІІ	0 1 2	x0 x1 x2
ІІІ	0 1 2 3	x0 x1 x2 x3