НУЛП, ІКНІ, САП
Тема
оцінка
підпис
КН-24
6
МЕТОДИ ЧИСЕЛЬНОГО ІНТЕГРУВАННЯ ФУНКЦІЙ
Назаревич В.О.
№ залікової: 1408545
Чисельні методи
Викладач:
к.т.н., асистент
Мельник М. Р.
Мета: Ознайомлення з методами чисельного інтегрування диференціальних рівнянь та їх практичним застосуванням.
Завдання:
1. Ознайомитись з методами чисельного інтегрування диференціальних рівнянь та їх практичним застосуванням.
2. Одержати індивідуальне завдання.
3. Знайти розв'язки заданого рівняння методом Ейлера, Рунге-Кутта, модифікованого методу Ейлера.
4. Порівняти ефективність і точність даних методів.
Варіант – 3
Хід роботи: Для виконання лабораторного завдання створюємо програму в MATLAB для знаходження розв’язку диференціального рівнняння.
Код програми:
clc
clear all;
a = 1;
b = 2;
f = @(x, y)y-(sin(x.^0.5)-cos(x.^0.5)/2*x.^0.5);
H_ME_01 = 0.1;
X_ME_01 = a:H_ME_01:b;
Y_ME_01 = zeros(1, length(X_ME_01));
Y_ME_01(1) = 0;
for i = 1:(length(Y_ME_01) - 1)
Y_ME_01(i + 1) = Y_ME_01(i) + H_ME_01*(f(X_ME_01(i), Y_ME_01(i)));
end
H_ME_001 = 0.01;
X_ME_001 = a:H_ME_001:b;
Y_ME_001 = zeros(1, length(X_ME_001));
Y_ME_001(1) = 0;
for i = 1:(length(Y_ME_001) - 1)
Y_ME_001(i + 1) = Y_ME_001(i) + H_ME_001*(f(X_ME_001(i), Y_ME_001(i)));
end
H_ME_0001 = 0.001;
X_ME_0001 = a:H_ME_0001:b;
Y_ME_0001 = zeros(1, length(X_ME_0001));
Y_ME_0001(1) = 0;
for i = 1:(length(Y_ME_0001) - 1)
Y_ME_0001(i + 1) = Y_ME_0001(i) + H_ME_0001*(f(X_ME_0001(i), Y_ME_0001(i)));
end
H_MRK_01 = 0.1;
X_MRK_01 = a:H_MRK_01:b;
Y_MRK_01 = zeros(1, length(X_MRK_01));
Y_MRK_01(1) = 0;
for i = 1:(length(Y_MRK_01) - 1)
K0 = H_MRK_01 * f(X_MRK_01(i), Y_MRK_01(i));
K1 = H_MRK_01 * f(X_MRK_01(i) + H_MRK_01/2, Y_MRK_01(i) + K0/2);
K2 = H_MRK_01 * f(X_MRK_01(i) + H_MRK_01/2, Y_MRK_01(i) + K1/2);
K3 = H_MRK_01 * f(X_MRK_01(i) + H_MRK_01, Y_MRK_01(i) + K2);
Y_MRK_01(i+1) = Y_MRK_01(i) + (K0 + 2*K1 + 2*K2 + K3)/6 ;
end
H_MRK_001 = 0.01;
X_MRK_001 = a:H_MRK_001:b;
Y_MRK_001 = zeros(1, length(X_MRK_001));
Y_MRK_001(1) = 0;
for i = 1:(length(Y_MRK_001) - 1)
K0 = H_MRK_001 * f(X_MRK_001(i), Y_MRK_001(i));
K1 = H_MRK_001 * f(X_MRK_001(i) + H_MRK_001/2, Y_MRK_001(i) + K0/2);
K2 = H_MRK_001 * f(X_MRK_001(i) + H_MRK_001/2, Y_MRK_001(i) + K1/2);
K3 = H_MRK_001 * f(X_MRK_001(i) + H_MRK_001, Y_MRK_001(i) + K2);
Y_MRK_001(i+1) = Y_MRK_001(i) + (K0 + 2*K1 + 2*K2 + K3)/6;
end
H_MRK_0001 = 0.001;
X_MRK_0001 = a:H_MRK_0001:b;
Y_MRK_0001 = zeros(1, length(X_MRK_0001));
Y_MRK_0001(1) = 0;
for i = 1:(length(Y_MRK_0001) - 1)
K0 = H_MRK_0001 * f(X_MRK_0001(i), Y_MRK_0001(i));
K1 = H_MRK_0001 * f(X_MRK_0001(i) + H_MRK_0001/2, Y_MRK_0001(i) + K0/2);
K2 = H_MRK_0001 * f(X_MRK_0001(i) + H_MRK_0001/2, Y_MRK_0001(i) + K1/2);
K3 = H_MRK_0001 * f(X_MRK_0001(i) + H_MRK_0001, Y_MRK_0001(i) + K2);
Y_MRK_0001(i+1) = Y_MRK_0001(i) + (K0 + 2*K1 + 2*K2 + K3)/6;
end
H_MRK = 1/100000;
X_MRK = a:H_MRK:b;
Y_MRK = zeros(1, length(X_MRK));
Y_MRK(1) = 0;
for i = 1:(length(Y_MRK) - 1)
K0 = H_MRK * f(X_MRK(i), Y_MRK(i));
K1 = H_MRK * f(X_MRK(i) + H_MRK/2, Y_MRK(i) + K0/2);
K2 = H_MRK * f(X_MRK(i) + H_MRK/2, Y_MRK(i) + K1/2);
K3 = H_MRK * f(X_MRK(i) + H_MRK, Y_MRK(i) + K2);
Y_MRK(i+1) = Y_MRK(i) + (K0 + 2*K1 + 2*K2 + K3)/6 ;
end
H_MME_01 = 0.1;
X_MME_01 = a:H_MME_01:b;
Y_MME_01 = zeros(1, length(X_MME_01));
Y_MME_01(1) = 0;
for i = 1:(length(Y_MME_01) - 1)
Y_MME_01(i + 1) = Y_MME_01(i) + H_MME_01*(f(X_MME_01(i) + H_MME_01/2, Y_MME_01(i) + ((H_MME_01*f(X_MME_01(i), Y_MME_01(i)))/2)));
end
H_MME_001 = 0.01;
X_MME_001 = a:H_MME_001:b;
Y_MME_001 = zeros(1, length(X_MME_001));
Y_MME_001(1) = 0;
for i = 1:(length(Y_MME_001) - 1)
Y_MME_001(i + 1) = Y_MME_001(i) + H_MME_001*(f(X_MME_001(i) + H_MME_001/2, Y_MME_001(i) + ((H_MME_001*f(X_MME_001(i), Y_MME_001(i)))/2)));
end
H_MME_0001 = 0.001;
X_MME_0001 = a:H_MME_0001:b;
Y_MME_0001 = zeros(1, length(X_MME_0001));
Y_MME_0001(1) = 0;
for i = 1:(length(Y_MME_0001) - 1)
Y_MME_0001(i + 1) = Y_MME_0001(i) + H_MME_0001*(f(X_MME_0001(i) + H_MME_0001/2, Y_MME_0001(i) + ((H_MME_0001*f(X_MME_0001(i), Y_MME_0001(i)))/2)));
end
P_ME_01 = zeros(1, length(X_ME_01));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_ME_01))
if X_MRK(q) == X_ME_01(w)
P_ME_01(e) = abs(Y_MRK(q) - Y_ME_01(w));
e = e + 1;
end
end
end
P_MRK_01 = zeros(1, length(X_MRK_01));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MRK_01))
if X_MRK(q) == X_MRK_01(w)
P_MRK_01(e) = abs(Y_MRK(q) - Y_MRK_01(w));
e = e + 1;
end
end
end
P_MME_01 = zeros(1, length(X_MME_01));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MME_01))
if X_MRK(q) == X_MME_01(w)
P_MME_01(e) = abs(Y_MRK(q) - Y_MME_01(w));
e = e + 1;
end
end
end
V_ME_01 = zeros(1, length(P_ME_01));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_ME_01))
if X_MRK(q) == X_ME_01(w)
V_ME_01(e) =abs( P_ME_01(w)/Y_MRK(q));
e = e + 1;
end
end
end
V_ME_01(1) = 0;
V_MRK_01 = zeros(1, length(P_MRK_01));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MRK_01))
if X_MRK(q) == X_MRK_01(w)
V_MRK_01(e) = abs( P_MRK_01(w)/Y_MRK(q));
e = e + 1;
end
end
end
V_MRK_01(1) = 0;
V_MME_01 = zeros(1, length(P_MME_01));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MME_01))
if X_MRK(q) == X_MME_01(w)
V_MME_01(e) = abs( P_MME_01(w)/Y_MRK(q));
e = e + 1;
end
end
end
V_MME_01(1) = 0;
P_ME_001 = zeros(1, length(X_ME_001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_ME_001))
if X_MRK(q) == X_ME_001(w)
P_ME_001(e) = abs(Y_MRK(q) - Y_ME_001(w));
e = e + 1;
end
end
end
P_MRK_001 = zeros(1, length(X_MRK_001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MRK_001))
if X_MRK(q) == X_MRK_001(w)
P_MRK_001(e) = abs(Y_MRK(q) - Y_MRK_001(w));
e = e + 1;
end
end
end
P_MME_001 = zeros(1, length(X_MME_001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MME_001))
if X_MRK(q) == X_MME_001(w)
P_MME_001(e) = abs(Y_MRK(q) - Y_MME_001(w));
e = e + 1;
end
end
end
V_ME_001 = zeros(1, length(P_ME_001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_ME_001))
if X_MRK(q) == X_ME_001(w)
V_ME_001(e) = abs(P_ME_001(w)/Y_MRK(q));
e = e + 1;
end
end
end
V_ME_001(1) = 0;
V_MRK_001 = zeros(1, length(P_MRK_001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MRK_001))
if X_MRK(q) == X_MRK_001(w)
V_MRK_001(e) = abs( P_MRK_001(w)/Y_MRK(q));
e = e + 1;
end
end
end
V_MRK_001(1) = 0;
V_MME_001 = zeros(1, length(P_MME_001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MME_001))
if X_MRK(q) == X_MME_001(w)
V_MME_001(e) = abs( P_MME_001(w)/Y_MRK(q));
e = e + 1;
end
end
end
V_MME_001(1) = 0;
P_ME_0001 = zeros(1, length(X_ME_0001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_ME_0001))
if X_MRK(q) == X_ME_0001(w)
P_ME_0001(e) = abs(Y_MRK(q) - Y_ME_0001(w));
e = e + 1;
end
end
end
P_MRK_0001 = zeros(1, length(X_MRK_0001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MRK_0001))
if X_MRK(q) == X_MRK_0001(w)
P_MRK_0001(e) = abs(Y_MRK(q) - Y_MRK_0001(w));
e = e + 1;
end
end
end
P_MME_0001 = zeros(1, length(X_MME_0001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MME_0001))
if X_MRK(q) == X_MME_0001(w)
P_MME_0001(e) = abs(Y_MRK(q) - Y_MME_0001(w));
e = e + 1;
end
end
end
V_ME_0001 = zeros(1, length(P_ME_0001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_ME_0001))
if X_MRK(q) == X_ME_0001(w)
V_ME_0001(e) = abs( P_ME_0001(w)/Y_MRK(q));
e = e + 1;
end
end
end
V_ME_0001(1) = 0;
V_MRK_0001 = zeros(1, length(P_MRK_0001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MRK_0001))
if X_MRK(q) == X_MRK_0001(w)
V_MRK_0001(e) = abs( P_MRK_0001(w)/Y_MRK(q));
e = e + 1;
end
end
end
V_MRK_0001(1) = 0;
V_MME_0001 = zeros(1, length(P_MME_0001));
e = 1;
for q = 1:(length(X_MRK))
for w = 1:(length(X_MME_0001))
if X_MRK(q) == X_MME_0001(w)
V_MME_0001(e) = abs(P_MME_0001(w)/Y_MRK(q));
e = e + 1;
end
end
end
V_MME_0001(1) = 0;
plot(X_ME_01, Y_ME_01, ':k', X_MRK_01, Y_MRK_01, '--b', X_MME_01, Y_MME_01, '-.g', X_MRK, Y_MRK, '-r')
title('Графіки при H = 0.1');
legend('Метод Ейлера', 'Метод Рунге - Кутта', 'Модифікований метод Ейлера', 'Точний графік');
figure;
plot(X_ME_001, Y_ME_001, ':k', X_MRK_001, Y_MRK_001, '--b', X_MME_001, Y_MME_001, '-.g', X_MRK, Y_MRK, '-r')
title('Графіки при H = 0.01');
legend('Метод Ейлера', 'Метод Рунге - Кутта', 'Модифікований метод Ейлера', 'Точний графік');
figure;
plot(X_ME_0001, Y_ME_0001, ':k', X_MRK_0001, Y_MRK_0001, '--b', X_MME_0001, Y_MME_0001, '-.g', X_MRK, Y_MRK, '-r')
title('Графіки при H = 0.001');
legend('Метод Ейлера', 'Метод Рунге - Кутта', 'Модифікований метод Ейлера', 'Точний графік');
figure;
plot(X_ME_01, P_ME_01, ':k', X_MRK_01, P_MRK_01, '--b', X_MME_01, P_MME_01, '-.g')
title('Абсолютна похибка для Н = 0.1');
legend('Метод Ейлера', 'Метод Рунге - Кутта', 'Модифікований метод Ейлера');
figure;
plot(X_ME_001, P_ME_001, ':k', X_MRK_001, P_MRK_001, '--b', X_MME_001, P_MME_001, '-.g')
title('Абсолютна похибка для Н = 0.01');
legend('Метод Ейлера', 'Метод Рунге - Кутта', 'Модифікований метод Ейлера');
figure;
plot(X_ME_0001, P_ME_0001, ':k', X_MRK_0001, P_MRK_0001, '--b', X_MME_0001, P_MME_0001, '-.g')
title('Абсолютна похибка для Н = 0.001');
legend('Метод Ейлера', 'Метод Рунге - Кутта', 'Модифікований метод Ейлера');
figure;
plot(X_ME_01, V_ME_01, ':k', X_MRK_01, V_MRK_01, '--b', X_MME_01, V_MME_01, '-.g')
title('Відносна похибка для Н = 0.1');
legend('Метод Ейлера', 'Метод Рунге - Кутта', 'Модифікований метод Ейлера');
figure;
plot(X_ME_001, V_ME_001, ':k', X_MRK_001, V_MRK_001, '--b', X_MME_001, V_MME_001, '-.g')
title('Відносна похибка для Н = 0.01');
legend('Метод Ейлера', 'Метод Рунге - Кутта', 'Модифікований метод Ейлера');
figure;
plot(X_ME_0001, V_ME_0001, ':k', X_MRK_0001, V_MRK_0001, '--b', X_MME_0001, V_MME_0001, '-.g')
title('Відносна похибка для Н = 0.001');
legend('Метод Ейлера', 'Метод Рунге - Кутта', 'Модифікований метод Ейлера');
Запускаємо програму в середовищі MATLAB і отримуємо графіки розвязку диференціального рівняння, а також графіки абсолютної та відносної похибок.
Результати: Графіки для точності 0.1 наведено на Рис. 1:
Рис. 1 Графіки для точності 0.1
Графіки для точності 0.01 наведено на Рис. 2:
Рис. 2 Графіки для точності 0.01
Графіки для точності 0.001 наведено на Рис. 3:
Рис. 3 Графіки для точності 0.001
Графіки абсолютної похибки для точності 0.1 наведено на Рис. 4:
Рис. 4 Графіки абсолютної похибки для точності 0.1
Графіки абсолютної похибки для точності 0.01 наведено на Рис. 5:
Рис. 5 Графіки абсолютної похибки для точності 0.01
Графіки абсолютної похибки для точності 0.001 наведено на Рис. 6:
Рис. 6 Графіки абсолютної похибки для точності 0.001
Графіки відносної похибки для точності 0.1 наведено на Рис. 7:
Рис. 7 Графіки відносної похибки для точності 0.1
Графіки відносної похибки для точності 0.01 наведено на Рис. 8:
Рис. 8 Графіки відносної похибки для точності 0.01
Графіки відносної похибки для точності 0.001 наведено на Рис. 9:
Рис. 9 Графіки відносної похибки для точності 0.001
Висновки: В ході виконання цієї лабораторної роботи я ознайомився з методами чисельного інтегрування диференціальних рівнянь та їх практичним застосуванням.
Написав програму в середовищі MATLAB і перевірив її роботу на прикладі даного мені диференціального рівняння. Кожен метод було перевірено тричі (для точностей 0.1, 0.01, 0.001). Проаналізувавши отримані результати, я зробив висновок, що метод Рунге-Кутта є найточнішим, а метод Ейлера дає найбільшу похибку. Тому метод Рунге-Кутта (з точністю 1/100000) був використаний мною для обчислення точного графіку і знаходження похибок.