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

Розглянемо створення файла-розв’язку для вирішення диференційного рівняння, для якого був створений файл-функція з використанням солвера 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('Рішення диференційного рівняня другого порядку');

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

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

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

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

на проміжку часу [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) %вивід легенди

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

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

Як видно з рисунка, при точності більшої за 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) %вивід легенди

Результат розв’язків диференційного рівняння для різних значень точності показаний на рис. 6.29. Як видно з рисунка, для заданого диференційного рівняння третього порядку достатньою точністю є 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) %вивід легенди

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

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

Література: [1], стор. 209-220; [2], стор. 71-78.

Завдання на СРС. Рішення рівняння Лотка-Вольтерра.

Література: [1], стор. 212-213.

Контрольні запитання

1. Яка послідовність вирішення диференціального рівняння у MatLab за допомогою М-файлів?

2. Які основні солвери існують та коли вони застосовуються?

3. Як створити файл-функцію? Послідовність запису файлу.

4. Як створити файл-розв’язок? Послідовність запису файлу.

5. Що таке точність обчислення? Як її змінити?

6. Як вивести на екран графічного розв’язку різні змінні рівняння?

Лекція 9

Поиск по сайту: