Лекція 1.3 MATLAB у задачах прикладної математики

4. Чисельний розв’язок звичайних диференціальних рівнянь

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

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

Постановка задачі:

Потрібно знайти функцію y = f(х), що задовольняє рівнянню:

                                              (3.3)

 і приймаючу при х = х₀ задане значення y₀:

                                                (3.4)

 При цьому рішення необхідно одержати в інтервалі х₀ £ х £ х_к.. З теорії диференційних рівнянь відомо, що розв’язок y(х) задачі Коші (3.3), (3.4) існує, і є гладкою функцією, якщо права частина F(x, y) задовольняє деяким умовам гладкості.

Чисельне рішення задачі Коші методом Рунге-Кутта 4-го порядку полягає в наступному. На заданому інтервалі [х₀, х_к] вибираються вузлові точки. Значення розв’язку в нульовій точці відомо y(х₀) = у_0. У наступній точці y(х₁) визначається за формулою:

      (3.5)

 тут

(3.6)

h – крок.                                                                                        

тобто, такий варіант методу Рунге-Кутта вимагає на кожнім кроці чотириразового обчислення правої частини рівняння (3.3). Цей алгоритм реалізовано у програмі ode45. Крім цієї програми MATLAB має у своєму розпорядженні великий набір аналогічних програм, що дозволяють успішно вирішувати звичайні диференційні рівняння.

Приклад.

Вирішити задачу Коші:

,             (3.7)

 Точний розв’язок має вигляд:

Розв’яжемо задачу за допомогою програми ode45. Спочатку в М-файл записуємо праву частину рівняння (3.7), сам М-файл оформляється як файл – функція, даємо йому ім'я F:

function dydx = F(x, y)

dydx = zeros(1,1);

dydx(1) = 2*(x^2+y(1));

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

Протокол програми,

>>[X Y] = ode45 ( @ F , [0 1] , [1] ) ;

% Дескриптор @ забезпечує зв'язок з файлом – функцією правої частини

% [0 1] – інтервал на якому необхідно одержати розв’язок

% [1] – початкове значення розв’язоку

>> рlot (X,Y);

>> % Побудова графіка чисельного рішення задачі Коші (3.7)

>> hold on; gtext ( ¢ y(x) ¢)

% Команда дозволяє за допомогою мишки нанести на графік напис у(х)

>> [X  Y]

>> % Остання команда виводить таблицю чисельного рішення задачі.

Результати рішення. Графік розв’язку задачі Коші (3.7) показаний на рисунку 3.3. Чисельне рішення представлене в таблиці 3.4, де наведено лише окремі вузлові точки. У програмі ode45 за замовчуванням інтервал розбивається на 40 точок із кроком h = 1/40 = 0,025.

Рис. 3.3. Графік рішення задачі Коші

Таблиця 3.4

Табличне представлення чисельного рішення задачі Коші

Як продемонстровано у таблиці 3.4 чисельне рішення програмою ode45 є точним.