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

Метод Рунге-Кутта є одним з методів підвищеної точності, але має багато загального з методом Ейлера.

Нехай на відрізку потрібно знайти чисельний розв’язок рівняння

(9.21)

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

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

(9.22)

Якщо розкласти функцію в у ряд Тейлора й обмежитися членами до h₄ включно, то збільшення функції у можна представити у вигляді

(9.23)

де похідні визначаються послідовним диференціюванням з рівняння (9.4).

Замість безпосередніх обчислень по формулі (9.23) у методі Рунге-Кутта визначаються чотири числа:

(9.24)

Можна довести, що якщо числам k₁, k₂, k₃, k₄ додати відповідно вагу 1/6; 1/3; 1/3; 1/6, те середньозважене цих чисел, тобто

(9.25)

з точністю до четвертих ступенів дорівнює значенню у, обчисленому по формулі (9.23):

(9.26)

Таким чином, для кожної пари поточних значень x_і і y_і по формулах (9.24) визначаються значення

По формулі (9.26) знаходиться

(9.27)

і потім

Числа k₁, k₂, k₃, k₄ мають простий геометричний зміст. Нехай крива M₀CM₁ (рис. 9.11) являє собою розв’язок диференціального рівняння (9.21) з початковою умовою (9.22). Точка С цієї кривої лежить на прямій, рівнобіжній осі Оу і розділяючою відрізок навпіл, В и G - точки перетинання дотичної, проведеної до кривої в точці М₀, з ординатами АС і N₁M₁,. Тоді число k₁ з точністю до множника h (де ) є кутовий

коефіцієнт дотичної в точці М₀ до інтегральної кривої M₀CM₁ тобто

Точка В має координати , . Отже число k₂ з точністю до множника є кутовий коефіцієнт дотичної, проведеної до ін тегральної кривої в точці В (BF - відрізок цієї дотичної).

Рисунок 9.11 – Геометрична інтерпретація метода Рунге-Кутта

Через точку М₀ проведемо пряму, рівнобіжну відрізку BF. Тоді точка D має координати і k₃ з точністю до множника h - кутовий коефіцієнт дотичної, перевіреної до інтегральної кривої в точці D (DR₁ - відрізок цієї дотичної). Нарешті, через точку М₀ проведемо пряму, рівнобіжну DR₁, що перетне продовження M₁N₁ у точці .

Тоді k₄ з точністю до множника h є кутовий коефіцієнт дотичної, проведеної до інтегральної кривої в точці R₂.

Обчислення по методу Рунге-Кутта зручно розташовувати за алгоритмом,

1.Значення х₀ и у₀ підставляють у праву частину диференціального рівняння (9.21), визначають .

2.Отримане значення множать на крок інтегрування h, обчислюють

.

3. Змінюють значення х₀

4. визначають допоміжне значення

5. визначення похідної в точці (, )

6. визначають значення k₂

7. визначають нове допоміжне значення

8. визначення похідної в точці (, )

9. визначають значення k₃

10. визначають нове значення допоміжного у_1д

11. змінюють значення

12. визначають допоміжну похідну в точці ()

13. визначають значення k₄

14. визначають нове значення у₁ за формулою

Для визначення повторюють ітераційний процес, починаючи з першого кроку.

Схема схема алгоритму метода Рунге-Кутта представлена на рисунку 9.12.

Рисунок 9.12 – Схема алгоритму метода Рунге – Кутта

Потім всі обчислення повторюють починаючи з І кроку, доти, поки не буде пройдений весь відрізок .

Метод Рунге - Кутта має порядок точності h⁴ на усьому відрізку . Оцінка точності методу цього дуже складна. Грубу оцінку погрішності можна одержати за допомогою "подвійного прорахунку" по формулі

(9.28)

де - значення точного розв’язку рівняння (9.21) у точці наближені значення, отримані з кроком h/ 2 і h. Якщо - задана точність розв’язку, то число п (число розподілів) для визначення кроку інтегрування вибирається таким чином, щоб

(9.29)

Однак крок розрахунку можна змінювати при переході від однієї точки до іншої.

Для оцінки правильності вибору кроку h використовують рівність

(9.30)

де q повинно бути дорівнює декільком сотим, у противному випадку крок h зменшують.

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

(9.31)

У цьому випадку паралельно визначаються числа

(9.33)

(9.34)

Тоді одержимо розв’язок системи диференційних рівнянь другого порядку:

Отже,

і остаточно одержуємо значення шуканих функцій у точці х = 0,6:

Приклад 9.1. Розв’язати на ЕОМ систему рівнянь виду:

При розробці програми для розв‘язку системи диференційних рівнянь, створюють: 1) підпрограму, яка містить праві частини рівнянь; 2) підпрограму, яка реалізує алгоритм метода Рунге-Кутта для обчислення значень функцій системи в одній точці; 3) підпрограму, яка дозволяє отримати розв‘язок системи на заданому відрізку дослідження. Наприклад, підпрограму, яка містить дану систему можна записати на алгоритмічній мові С у вигляді:

void systema (float *Y, float *F, float X)

{ ;

;

}

де - значення правої частини рівняння заданої системи.

Приклад 9.2. Розглянемо приклад розв’язання відомої інженерної задачі на ЕОМ за допомогою методу Рунге-Кутта.

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

Рисунок 9.13 – Графічне представлення використання пружних опор в транспортних засобах

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

Пристрій, показаний на рисунку 9.13, складається з маси т, зв'язаною з жорсткою стінкою через пружину постійної жорсткості k, демпфер з коефіцієнтом демпфірування с й пружиною з нелінійною характеристикою утворюючої відновлюючу силу, рівну добутку постійної k на зміщення в третій степені. Така пружина має симетричну нелінійну характеристику, забезпечуючи захист від ударних й вібраційних навантажень. Так як рух цієї системи описується нелінійним диференціальним рівнянням:

,

то традиційні «точні» методи не дозволяють знайти залежність зміщення х від часу. Тому для розв'язку вказаного рівняння прийдеться скористатись чисельним методом. Нехай параметри системи мають наступні значення:

k =2,0 Н/см; k* = 2, 0 Н/см³; с = 0,15 Нс/см; m =1,0 кг,

а початкові умови задані в вигляді: х(0) - 10 см, х(0)=0.

Підготуємо й реалізуємо на ЕОМ програму для моделювання руху даної механічної системи в інтервалі часу 0 t 1.

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

Щоб при видачі довжина виражалась в сантиметрах, змінимо розмірністьмаси. Тоді m=0,01 Н?с² /см.

Рисунок 9.14 – Графік представлення результатів чисельного розв’язку приклада 9.11 на ЕОМ

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

9.4 Методи прогнозу і корекції (багатоточкові методи)

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

(9.35)

Рисунок 9.15 – Схема алгоритму багатоточкового методу

Розглянемо особливості алгоритмів методів прогнозу і корекції.

Так як в методах, що розглядаються використовується інформація про декілька раніше отриманих точок, то на відміну від однокрокових методів вони не володіють властивістю „самостартування”. Тому, перед тим як застосовувати метод прогнозу і корекції, необхідно обчислювати вихідні данні за допомогою якого-небудь однокрокового методу. Часто для цього використовують метод Рунге-Кутта. Обчислення проводять наступним чином. Спочатку по формулі прогнозу і початковим значенням змінних знаходять значення . Верхній індекс (0) означає, що прогнозоване значення є одним з послідовності значень , що розташовані в порядку зростання точності. По прогнозованому значенню за допомогою приведеного вище диференціального рівняння знаходять похідну

(9.36)

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

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

. (9.37)

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

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

Якщо диференціальне рівняння про інтегровано в інтервалі значень від до , то результат прийме вигляд

(9.38)

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

Дата добавления: 2014-11-08