9.1.1 Модифікації методу Ейлера

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

Суть удосконаленого методу Ейлера полягає в використанні ітераційної формули виду:

,          (9.16)

де - значення аргументу х в точці , а - значення функції в точці .

Розглянемо диференціальне рівняння

          (9.17)

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

Потрібно знайти розв’язок рівняння (9.13) на відрізку .

Розіб'ємо відрізок на n рівних частин точками , де – крок інтегрування

Алгоритм методу складається з:

1. визначення похідної в точці :

2. змінна х за формулою:

3. визначення значення при

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

5. використовуємо отримане значення для визначення за формулою:

,

6. змінюємо

7. повторюємо всі кроки алгоритму, починаючи з першого.

Графічна інтерпретація методу представлена на рисунку 9.7.

    

     Рисунок 9.7 – Графічна інтерпретація удосконаленого методу Ейлера

Зауваження. Оцінка похибки в точці x_і може бути отримана за допомогою "подвійного прорахунку": розрахунок повторюють із кроком h/2 похибка більш точного значення (при кроці h/2) оцінюють приблизно в такий спосіб:

               (9.18)

де у(х) - точний розв’язок диференціального рівняння. Удосконалений метод Ейлера є більш точним у порівнянні з методом Ейлера та відноситься до методів 3-го порядку точності.

Схема алгоритму представлена на рисунку 9.8.

    

Рисунок 9.8 - Схема алгоритму удосконаленого метода Ейлера

Модифікований метод Ейлера заснований на використанні ітераційної формули виду:

               (9.19)

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

    

     Рисунок 9.9 – Графічна інтерпретація

Алгоритм методу включає наступні кроки:

1. визначення похідної в точці :

2. зміна незалежної змінної х за формулою:

3. визначення проміжного значення за формулою методу Ейлера

4. визначення проміжної похідної в точці

5. визначення середньо арифметичного значення двох похідних

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

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

Схема алгоритму метода представлена на рисунку 9.10.

    

     Рисунок 9.10 - Схема алгоритму модифікованого метода Ейлера

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

               (9.20)

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

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

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

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

          (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 на ЕОМ

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