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

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

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

та початковими умовами при

,

Задача такого плану вирішується у MatLab за допомогою наступної послідовності етапів:

1. Зведення диференційного рівняння до системи диференційних рівнянь першого порядку.

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

3. Вибір потрібного солвера.

4. Створення спеціального файла-розв’язку для візуалізації результатів.

Розглянемо приклад створення файла-функції для диференційного рівняння другого порядку з заданими початковими умовами , .

Для приведення рівняння до системи диференційних рівнянь першого порядку виконується заміна , . Після підстановки у задане рівняння та отримаємо наступну систему

Далі необхідно створити файл-функцію. Він повинен мати два вхідних аргументи: змінну , за якою буде відбуватися диференціювання, та вектор, розмір якого дорівнює кількості невідомих функцій системи (у даному прикладі дві: та ). Число та порядок аргументу фіксовані, навіть якщо змінна явно не входить до системи. Вихідним аргументом файла-функції є вектор правої частини системи.

Таким чином текст файла-функції буде мати наступний вигляд

function F=dif(t, y); %об’ява функції для

%розв’язку системи

F=[y(2); -2*y(2)-10*y(1)+sin(t)]; %запис правої

%частини системи

По завершенню формування файла-функції його треба зберегти під ім’ям, яким названа функція F, тобто «dif.m».

Основні солвери для розв’язку диференційних рівнянь з заданими початковими умовами.

Найрозповсюджені солвери та області їх застосування наведені у таблиці 5.1.

Таблиця 5.1.

Назва солвера	Область застосування
ode45	Дуже розповсюджений солвер. Базується на формулах Рунге-Кутта четвертого та п’ятого порядку точності. Часто дає задовільні результати.
ode23	Також базується на формулах Рунге-Кутта, але менш низького порядку точності. Застосовується у нежорстких задачах, коли потрібно отримати розв’язок з невисоким ступенем точності.
ode113	Базується на методі Адамса-Бешфорта-Мілтона та використовується для отримання розв’язку високої точності у нежорстких задачах. Цей солвер особливо ефективний для нежорстких систем диференційних рівнянь, праві частини яких обраховуються за складними формулами.
ode15s	Оснований на методі Гіра, який допускає зміну порядку, та застосовується для вирішення жорстких систем.
ode23s	Базується на однокроковому методі Розенброка другого порядку та застосовується для вирішення жорстких систем з невеликою точністю.

 

Усі солвери намагаються знайти розв’язок з відносною точністю 10^-3. Якщо потрібно збільшити точність обрахування, необхідно застосувати додатковий параметр options, який формується функцією odeset. Наприклад, якщо потрібна точність 10^-5, то у файл-розв’язок, що буде розглянутий пізніше, потрібно на початку включити додатковий рядок

options=odeset('RelTol', 1.0e-05)

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

Створення файла-розв’язку для чисельного вирішення диференційних рівнянь.

Розглянемо створення файла-розв’язку для вирішення диференційного рівняння, для якого був створений файл-функція з використанням солвера ode45. Вхідними аргументами солвера будуть: ім’я файла-функції, що записується у апострофах 'dif' або без них, але з символом @ на початку (@dif); вектор з початковим і кінцевим значенням часу диференціювання; вектор початкових умов. Запис у М-файлі буде наступним

[T, Y]=ode45('dif', [Time], [Y0])

де масиви [Time] та [Y0] – масиви часу та початкових умов відповідно. Ці масиви можуть записуватися безпосередньо у солвері, наприклад

[T, Y]=ode45('dif', [1 15], [1 0])

або задаватися раніше. Тоді у солвері використовуються тільки ім’я цих масивів

Time=[0 15];

Y0=[1 0];

[T, Y]=ode45('dif', Time, Y0);

При необхідності зміни точності за допомогою параметра options солвер записується наступним чином

options=odeset('RelTol', 1.0e-05)

Time=[0 15];

Y0=[1 0];

[T, Y]=ode45('dif', Time, Y0, options);

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

Текст програми, що записується у файл-розв’язок має вигляд

options=odeset('RelTol', 1.0e-04); %точність 0.0001

Time=[0 15]; %масив часу

Y0=[1 0]; %масив початкових умов

[T, Y]=ode45('dif', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'r', 'LineWidth', 2); %вивід на екран

%змінної y1

hold %затримка графіка

grid %вивід стіки

plot(T, Y(:, 2), 'b--', 'LineWidth', 2); %вивід на екран

%змінної y2

xlabel('Час t, с'); %мітки осей та підпис графіка

ylabel('y1, y2');

title('Рішення диференційного рівняня другого порядку');

Після запуску файла-розв’язку відкриється графічне вікно, зображене на рис. 5.1.

Рис. 5.1. Результати розв’язку диференційного рівняння

Треба мати на увазі, що обидва файл-функція і файл-розв’язок повинні знаходитися в одній директорії. Файл-розв’язок не буде працювати без файла-функції. Після створення та збереження файл-функцію не потрібно запускати на виконання на відміну від файла-розв’язку.

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

на проміжку часу [a, 100] для заданих початкових умов , , де =0.001.

Файл-функція для заданої системи рівнянь буде мати вигляд

function F=dif(t, y) %об’ява функції для

%розв’язку системи

F=[y(2); -1/t^2]; %запис правої

%частини системи

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

options=odeset('RelTol', 1.0e-03); %точність 0.001

Time=[0.001 100]; %масив часу

Y0=[log(0.001) 1/0.001]; %масив початкових умов

[T, Y]=ode45('dif', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'k: ', 'LineWidth', 2); %вивід на екран

%змінної y1

hold %затримка графіка

grid %вивід стіки

xlabel('Час t, с'); %мітки осей та підпис графіка

ylabel('y1');

title('Рішення диференційного рівняня другого порядку');

options=odeset('RelTol', 1.0e-04); %точність 0.0001

[T, Y]=ode45('dif', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'r--', 'LineWidth', 2); %вивід на екран

%змінної y1

options=odeset('RelTol', 1.0e-05); %точність 0.00001

[T, Y]=ode45('dif', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'g-', 'LineWidth', 2); %вивід на екран

%змінної y1

options=odeset('RelTol', 1.0e-06); %точність 0.000001

[T, Y]=ode45('dif', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'c-.', 'LineWidth', 2); %вивід на екран

%змінної y1

options=odeset('RelTol', 1.0e-08); %точність 0.00000001

[T, Y]=ode45('dif', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'b--', 'LineWidth', 2); %вивід на екран

%змінної y1

legend('10^{-3}', '10^{-4}', '10^{-5}', '10^{-6}', '10^{-8}', 4) %вивід легенди

Після запуску файла-функції на екрані з’явиться графічне вікно

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

 

Як видно з рисунка, при точності більшої за 10^-5 графіки співпадають, а при точності менше 10^-5 результати мають дуже велику розбіжність. Таким чином для даної системи диференційних рівнянь достатньою точністю для отримання точного рішення є 10^-6.

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

на проміжку [0 10] для початкових умов , , .

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

.

Створюємо файл функцію для отриманої системи з ім’ям «dif2.m»

function F=dif2(t, y); %об’ява функції для

%розв’язку системи

F=[y(2); y(3); -5*y(3)+2*y(2)+y(1)-14+cos(2*t-5)]; %запис правої частини системи

Після цього зберігається файл-функція та створюється наступний файл-розв’язок

options=odeset('RelTol', 1.0e-03); %точність 0.001

Time=[0 10]; %масив часу

Y0=[1 2 0]; %масив початкових умов

[T, Y]=ode45('dif2', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'k: ', 'LineWidth', 2); %вивід на екран

%змінної y1

Hold %затримка графіка

Grid %вивід стіки

xlabel('Час t, с'); %мітки осей та підпис графіка

ylabel('y1');

title('Рішення диференційного рівняння третього порядку');

options=odeset('RelTol', 1.0e-04); %точність 0.0001

[T, Y]=ode45('dif2', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'r--', 'LineWidth', 2); %вивід на екран

%змінної y1

options=odeset('RelTol', 1.0e-05); %точність 0.00001

[T, Y]=ode45('dif2', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'b-', 'LineWidth', 2); %вивід на екран

%змінної y1

legend('10^{-3}', '10^{-4}', '10^{-5}', 3) %вивід легенди

Результат розв’язків диференційного рівняння для різних значень точності показаний на рис. 5.3. Як видно з рисунка, для заданого диференційного рівняння третього порядку достатньою точністю є 10^-3. При збільшенні точності графіки співпадають.

Якщо потрібно вивести інші змінні диференційного рівняння, наприклад та , то файл-розв’язок для одного варіанту точності, що дорівнює 10^-3 буде мати наступний вигляд

options=odeset('RelTol', 1.0e-03); %точність 0.001

Time=[0 10]; %масив часу

Y0=[1 2 0]; %масив початкових умов

[T, Y]=ode45('dif2', Time, Y0, options); %солвер

plot(T, Y(:, 1), 'k: ', 'LineWidth', 2); %вивід на екран y1

hold %затримка графіка

grid %вивід стіки

xlabel('Час t, с'); %мітки осей та підпис графіка

ylabel('y1');

title('Рішення диференційного рівняння третього порядку');

plot(T, Y(:, 2), 'b--', 'LineWidth', 2); %вивід на екран y2

plot(T, Y(:, 3), 'k-.', 'LineWidth', 2); %вивід на екран y3

legend('y(1)', 'y(2)', 'y(3)', 3) %вивід легенди

Рис. 5.3. Графіки розв’язків диференційного рівняння третього порядку для різних значень точності

Рис. 5.4. Графіки усіх змінних розв’язку диференційного рівняння третього порядку