Реферат: Деякі скінченно-різнецеві методи розв’язування звичайних диференціальних рівнянь

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

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

Факультет математики та інформатики

кафедра диференціальних рівнянь і прикладноїматематики

Курсова робота на тему:

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

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

Виконав:

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

Васьків Святослав

Перевірив:

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

Василишин П.Б.

Івано-Франківськ 2010

План

Вступ.

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

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

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

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

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

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

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

д) Програма Beeman

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

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

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

Висновок

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

Вступ

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

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

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

/> (1)

/> (2)

Метою усіх скінцеворізнецевих методів являється знаходження значень xn+1і v n+1(точка в «фазовому просторі») у момент часу tn+1=tn+∆tНам вже відомо, що величину кроку ∆t потрібно вибирати таким чином,щоб метод інтегрування породжував приймати однакове рішення. Один із способівперевірки стійкості методу полягає в контролі величини повною енергії ізабезпеченні того, щоби вона не відхилялася від початкового значення у разі,коли повна енергія зберігалась. Досить велике значення кроку приводить до незбереження повної енергії і до різних розв’язків для хn+1 i vn+1, тобто до таких розв’язків, які все більше відхиляютьсяз потоком часу від істинного розв’язку.

Сутьбагатьох алгоритмів, можна зрозуміти, розкладаючи

vn+1≡ v(tn+∆t)i

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

/> (3)

/> (4)

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

/> (5)

/> (6)

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

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

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

/> (7)

/> (8)

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

/> (9)

/> (10)

Замітимо, щоякщо підставити вираз (9) для vn+1в(10) то отримаємо

/> (11)

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

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

/> (12)

/> (13)

Замітио,що метод напів кроку не являє собою «самостартючим», так чи однакще фомули недозволять порахвати v1/2.Цюнезручність модна здолати поклавши

/> (14)

Оскільки формули(12-14) можна повторювати до безмежності, то метод напів кроку отримав широкерозповсюдження в навчальній літературі.

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

(15)

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

/> (16)

або

/> (17 )

Аналогічно розв'язаннярозкладу в ряд Тейлора для xn+1ixn-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)для vn+1

/> (23)

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

/> (24)

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

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

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

/> (25а)

/> (25б)

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

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

/> (26а)

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

Тоді,використовуючи />, отримаємо скоректированезначенняvn+1іxn+1

коректор: />

/> (26б)

Скоректированезначення xn+1використовується для визначення нового передбачуваного значення аn+1і, значить, нових передбачених значень vn+1i xn+rЦя процедура повторяється до тих пір, доки передбачення і скоррективонанезначення xn+1відрізняються менше ніж на задану величину! Даний метод можна розробити насхемі більш високого порядку, які зв’язуються між собою не тільки xn+1, xnі vn, але і таксамо значеннями vn-1і vn-2.Замітимо, що метод предиктора-коректора не являєтьсясамостартуючим.

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

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

/> (27)

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

/> (28)

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

/> де/>

Застосуванняметоду Рунге-Кутти до рівнянь руху Ньютона дає

/> (29)

Оскільки методи Руиге-Кутти є такими, що самостартуючими, то їх частовикористовують для вираховання декількох перших кроків для несамостартуючихалгоритмів.

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

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

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

/> />

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

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

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

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

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

де

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

/> 0

c2 a21

c3 a31 a32

∙ ∙ ∙ ∙

cs as1 as2∙ ∙ ∙ as,s − 1

b1 b2 bs− 1 bs

Для коефіцієнтів методу Рунге — Кутта мають справджуватись умови

Якщо ми хочемо, щоб метод мав порядок p, то варто так самозабезпечити умову /> — наближення, отримане по методу Рунге — Кутти. Після багаторазовогодиференціювання ця умова перетвориться в систему поліноміальних рівнянь накоефіцієнти методу.

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

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

На першийпогляд, слід спочатку вичислити матрицю />,а потім застосовувати її потрібну кількість разів до усім fj.Проте, як ми знаєм (див. частину I, розділ „Системи лінійнихрівнянь“), це не є правильний спосіб рішення СЛР „про запас“.Правильний спосіб рішення СЛР з цією матрицею /> разі назавжди (для довільної правої частини) — це LU — розкладання матриці />.При цьому матриця /> розкладаєтьсяна ліву нижню трикутну матрицю L і праву верхню трикутну матрицю U. Після цьогорішення будь-хто СЛР з матрицею /> будуєтьсяза допомогою прямої підстановки (для матриці L) і потім зворотної підстановки(для матриці U). Кожна з підстановок вимагає N2 операцій. У нашомувипадку має сенс скомбінувати праві частини усіх СЛР так, щоб кількістьзастосувань матриці /> (точніше,кількість СЛР, які потрібно вирішити) була мінімальною. Для цього досить ввестипроміжні змінні />.В результаті ми отримаємо наступний рецепт для одного кроку за часом 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,…виконуємо наступне:

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

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

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

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

виконуємо поліноміальну інтерполяцію взалежності Y(X) в точку X = 0. Результатінтерполяції, отриманий на k -му етапі,означаємо />/>. Оцінюємо погрішність

/> як />

Перетворюємо вектор />начисло />характеризуючипогрішність. Якщо /> меншезамовленого /> тевважаємо /> ізавершуємо даний крок t→t+H.

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

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

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

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

Наведемо прикладпограми Рунге-Кутта на мові С#. В програмі використовуєтьсяабстрактний клас TrungeKutta, в якому потрібноперекрити абстрактний метод F, який задає перші чаcтини рівняння.

using System;

using System.Collections.Generic;

namespace rwsh_rk4

{

abstract class TRungeKutta

{

public int N;

double t; // теперішній час

public double[] Y; // шукане число Y[0] – саме рішення, Y[i] — i-та змінна розвязку

double[] YY, Y1, Y2, Y3, Y4; // внутрішня змінна

public TRungeKutta(int N) // N – розмір системи

{

this.N = N; // зберегти розміри системи

if (N < 1)

{

this.N = -1; // якщо розмірність менше одиниці, тоустановити флаг помилки

return; // і вийти із конструктора

}

Y = new double[N]; // створення вектора розв’язку

YY = new double[N]; // внутрішній розв’язок

Y1 = new double[N];

Y2 = new double[N];

Y3 = new double[N];

Y4 = new double[N];

}

public void SetInit(double t0, double[] Y0) // встановленняпочаткових умов.

{ // t0 – початковийчас, Y0 – початкова умова

t = t0;

int i;

for (i = 0; i < N; i++)

{

Y[i] = Y0[i];

}

public double GetCurrent() // повернути даний час

{

return t;

}

abstract public void F(double t, double[] Y, ref double[] FY); //перші частини с-ми.

public void NextStep(double dt) // наступний крок методаРунге-Кутта, dt – крок по часу

{

if(dt<0)

{

return;

}

int i;

F(t, Y, ref Y1); // вирахувати Y1

for (i = 0; i < N; i++)

{

YY[i] = Y[i] + Y1[i] * (dt / 2.0);

}

F(t + dt / 2.0, YY, ref Y2);

for (i = 0; i < N; i++)

{

YY[i] = Y[i] + Y2[i] * (dt / 2.0);

}

F(t + dt / 2.0, YY, ref Y3); // вирахувати Y3

for (i = 0; i < N; i++)

{

YY[i] = Y[i] + Y3[i] * dt;

}

F(t + dt, YY, ref Y4); // вирахувати Y4

for (i = 0; i < N; i++)

{

Y[i] = Y[i] + dt / 6.0 * (Y1[i] + 2.0 * Y2[i] + 2.0 *Y3[i] + Y4[i]); // виразувати р-зок на новому кроці

}

t = t + dt; // збільшити крок

}

class TMyRK: TRungeKutta

{

public TMyRK(int aN): base(aN) { }

public override void F(double t, double[] Y, ref double[] FY)

{

FY[0] = Y[1]; // приклад математичний маятник

FY[1] = -Y[0]; // y''(t)+y(t)=0

}

class Program

{

static void Main(string[] args)

{

TMyRK RK4 = new TMyRK(2);

double[] Y0 = {0, 1}; // задаємо початкові умови y(0)=0,y'(0)=1

RK4.SetInit(0, Y0);

while (RK4.GetCurrent() < 10) // розв’язуєм до 10

{

Console.WriteLine("{0}\t{1}\t{2}",RK4.GetCurrent(), RK4.Y[0], RK4.Y[1]); // вивести t, y, y'

RK4.NextStep(0.01); // розв’язати на наступному кроці,крок інтегрування dt=0.01

}

8. Програма Beeman

У програмі Вееman моделюється осцилятор Морза з допомогою алгоритмуБімана. Оскільки цей алгоритм не самостартуючий, то для всіх — числових значеньх(, і використовується швидкісна форма алгоритму Верле.

PROGRAM Beeman Імоделювання осцилятора Морза

CALL initialfx, v, aold, dt, dt2, nmax)

CALL energy(x,v, ecum,e2cum) 1значення початкової енергії

CALLVerleg x, v, a, aold, dt, dt2) CALL energy{ x, v, ecum, e2cum) LET n = 1

DO whiie n < nmax

LET n = n + 1 1 число кроку

CALL Bceman{x, v, a, aold, dl, dl2)

Іобразування повної знергії після кожного кроку за часом CALL energY(x, v, ecum,e2cum) LOOP

CALL output{ ecum, e2cum, n) END

SUB initial(x, v,aold, dt,dt2, nmax)DECLARE DEFf LETх = 2 LET v = 0 LET aold = f(x)

INPUT prompt «крок по часу(c) = »: dt LET dt2 = dt" dt

INPUT prompt «тривалість = »: tmax LETnmax = tmax/dt END SUB

SUB Ver!et( x, v, a, aold, dt, dt2) DECLARE DEF f

LET x = x + v*dt + 0.5*ao!d*dt2 LET a = f(x)

LET v = v + 0.5"{a + ao!d)*dt END SUB

SUB Beeman(x, v, a, aold, dt, dt2) DECLARE DEF f

LET x = x + v*dt + (4*a — aold)*dt2/6

LET anew = f(x) І значенняна(п+1) -му кроці LET v = v + (2*anew + 5*a — aold)*dt/6

LET aold= a значення на (n-1)-му кроці

LET а = anew значення на n-му кроці ENDSUB

DEF f(x) LET e = exp(- x) LET f = 2*e*(e — 1) ENDDEF

SUB energY(x, v, ecum, e2cuin) LET KE = 0.5*v"v LET e = exp(- x) LET PE = e*{e — 2) LET etot = KE + PE LET ecum = ecum + etotLET e2cum = e2cum + etot*etot END SUB

SUB output{ecum,e2cuiT!,n) LETn = n+ 1 І вирахування початкового значення

LET ebar = ecum/n PRINT «середня енергія = »;ebar LET sigma2 = e2cum/n- ebar*ebar PRINT «sigma = »; sqr(sigma2) END SUB

Метод Адамса

Цей метод чисельного інтегрування розроблений Адамсом в 1855році на прохання видомого англійського алтелериста Башфора, який займавсявнутрішньою балістикою. В подальшому цей метод був забутий і знову відкритийбув норвезьким математиком Штермером. Популяризація метода Адамса і подальшейого вдосконалення пов’язане із іменем Крилова.

Запишемо рівняння першого порядку

З початковими умовами /> (1,2)

Нехай xi(i=0,1,2…)-система рівнозначних значень з крокомhi y(xi).Очевидно маємо

/> (3)

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

/> (4)

де /> або /> (4а)

Підставляю вираз (4а) в формулу (3) і враховуючи те, що />будемо мати

З відси отримуємо формулу експоляриціональну Адамса

/> (5)

Для початкового процессу потрібно чотири початкових значенняy0, y1, y2, y3, — початковий відрізок,який приділяє, виходячи із початкових умов (2), яким-небуть чисельним методом.Мажна наприклад використати метод Рунге-Кутта або розкласти в ряд Тейлора

Де i=1,2,3 (або i=-1,1,2)із відповідною зміною нумерування. Знаючи ці значення, із рівнянь (1)можна знайти значення похідних /> і скласти таблицю

/> (6)

Подальше значення yi(i=4,5…) шуканого розвязку можна крокза кроком обчислювати за формулою Адамса, поповнюючи по мірі можливості таблицюрізниць (6)

Вирахувавши перше наближення для /> по формулі

Визначити /> підрахуватикінцеві різниці

/> (7)

а потім знайти друге наближення для більш точній формулі

/> (8)

/>/>

/>Якщо /> і /> відрізняються лишень на дкілька одиниць останнього зберігаючогодесяткового розряду, то можна поставити /> а потім знайшовши /> перерахувавши кінцеві різниці (7). Після цього, потрібно знову знайти />по формулі (8) Поту цей крок h повинен бути таким, щоб цейперерахунок був зміненим.

На практиці крок h вибирають малим, щоб можнабуло знехтувати членом /> в формулі (8)

Якщо за розбіжність величин />і />суттєва, то потрібно зменшити крок h.

Звичайно крок h зменшують рівно в 2 рази.Можна показати, як в цьому випадку, маючи до деякого значення і таблицювеличин хj, yj, Yj=hy’j(j<=i) з кроком />, можна просто побудувати таблицю величин /> з кроком />

На основі формули (4) будемо мати

/> (9)

Де />Звідси, /> і /> і враховуючи, що /> заходимо

/> (10)

Аналогічно при /> із формули (9) отримаєм, що аргументу /> відповідає значення

/> (11)

Що стосується значень Yi-1iYi,то вони знаходяться в старій таблиці. Після цього складаємо початковий відрізокдля нової таблиці:

і знаходимо кінцеві різниці:

/>/> />/>

Далі таблиця будується простимспособом, подальшою модифікацією формули (5):

Для роботи на компютерах формулу Адамса (5) вигідновикористовувати в розкритому виді. Враховуючи, що

Після цього маємо: /> причому/>

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

Для спрощеня запису обмежимось розглядом диференціальнихрівнянь першого порядка

/> (1)

З початковими умовами />

Введемо спочатку ряд допоміжних формул

В силу формули Адамса отримаємо

/> (2)

Введемо позначення />

Формула (2) називається формулою похилогорядка, так як в ній використовуютьсярізниці, які знаходяться на діагоналі таблиці різниць. Враховуючи, що />

Із формули (2) будемо мати

Звідси отримуємо першу допоміжну формулу – яку ще можнаназвати перша формула ламаного рядка

/> (3)

Далі враховуючи, що /> і/> із формули (3) виводимодругу формулу – друга формула ламаного рядка

/> (4)

Якщо /> отримаємоформулу горизонтального рядка

/> (5)

Підмітимо, що формулу (5) можна отримати безпосередньо задопомогою інтегрування, в межах відxiдоxi+1<sub/>розкладанням/> за допомогою першоїінтерполяційної формули Ньютона:

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

/> Післяцього знайдемо /> і складаєрізницю />, де />/>/>.

Значення які знайшли заносимо в розділ (І) основного бланку(таблиця 1)

Схема обчислення відрізка методом послідовних наближень

№ наближення

ІІ

ІІІ

Далі переходимо до другого наближення. Для того,використовуємо дані із знаходження ламаних рядків, обчислюємо значення /> і />:

/> (7)

/> (8)

Двочленні формули отримуються відповідно із формули (5) при і=0і із формули (2) при і=1 в результаті відкидання різниць порядка вищогоніж перший.

Таким чином, отримаємо можливість знайти

/> і /> ,

в результаті чого можна порахувати

і скласти різниці

Отримані результати записуємо у таблицю в розділ 2 основногобланка

Для знаходження третього наближення застосовуємо трьохчленніформули, які отримуються із формули (2) при і=2 після відкидання різницьтретього порядку. Обчислюємо значення /> із трьохчленних формул:

/> (9)

/> (10)

/> (11)

Звідси можна знайти

і обчислити />/>/>. Після цього можна заповнити розділ ІІІ в таблиці (І) знайшовши потрібнірізниці звичайним порядком.

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

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

задовольняючи подвійну нерівність

і початковими умовами /> причому такими, що /> на /> при />. Геометрично це означає, що шукана інтегральна крива /> стискається в як завгодно малий криволінійний сектор А0ВnCn (мал. 1).

/>Якщо положити /> то максимальна абсолютна погрішність наближеного розв’язку /> буде рівна />/> ця погрішність на кожному кроці визначається безпосередньо.

Покажемо ідею метода Чаплигіна для диф. рівнянь першогопорядку

/> (1)

з початковою умовою

/> (2)

Причому будемо мати на увазі, що права частина /> непереривна і має неперервні похідні /> і /> в деякому околі початкової точки />. Метод побудований на одній лемі.

Лема Чаплигіна про інтегральні нерівності.

Нехай /> - диференціальний оператор, який відповідає диференціальномурівнянню (1), і />інтеграл рівняння (1)

/> (3)

яке задовольняє початкову умову /> і вибраний при />.

Якщо функція />задовольняючи умови:

/> (4)

то на відрізку /> виконується нерівність

/> (5)

так чи однаке функція і являється наближенимрозв’язком /> .

Аналогічно і для функції />виконуються умови:

/> (6)

то на відрізку /> має місце нерівність />, (7)

так чи однаке функція /> являється верхнім наближеним розв’язком />у.

Доведення: Достатньо доказати лиш одне із нерівностей (5) або (7). Доведемонаприклад нерівність (5). Із формул (3) і (4) маємо /> і /> Звіди

/> (8)

Де

/> (9)

Функція /> втрачає зміст при х, для якого />. В цьому випадку

В силу наведених вище умов функція р(х) визначена інеперервна на відрізку />.

Помножимо обидві частини диференціальної нерівності (8) наінтегруючий множник />будемо мати

/> (10)

Звідси інтегруючи нерівність (10) в межах від /> до />, де /> отримаєм />, або так як /> то остаточно знаходимо /> при />, що і потрібно було довести.

194-Чисельні методи

Висновок:

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

Як вже підкреслювалося, не слід віддавати перевагу одному якому-небудьалгоритму. Успіхи в комп'ютерній технології нині дозволяють легкоексперементувати з різними алгоритмами для разнообразних динамічних систем.

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

1) Х. Гулд, Я.Тобочник. Компьютерноемоделирование в физике: часть1.

2) В.А.Ильина, П.К. Силаев. Численныеметоды для физиков-теоретиков часть2. Москва, Ижевск 2004р.

3) сайт uk.wikipedia.org/wiki/Метод_Рунге_—_Кутти

4) И.С. Березин, Н.П. Жидков. Методывычислений том. 2 Москва 1959р.

5) Б.П. Демидович, И.А.Марон, З.Шувалова.Численные методы анализа. «Наука» Москва1967р.

еще рефераты

Еще работы по математике

Реферат по математике

Методы математической статистики

20 Июня 2015

Реферат по математике