Задание 3.1. Решить задачу Коши для обыкновенного дифференциального уравнения первого порядка. Найти численное решения уравнения y′(x) = f (xy(x)) на заданном промежутке изменения аргумента при заданном начальном условии. Численное решение уравнения найти явными методами Эйлера и Рунге-Кутта, а также с помощью функции ode45 из пакета MatLab. Построить полученные различными методами решения на разных графиках, а также на одном отдельном графике сравнить между собой все решения.y′′ + 2y′ + 2y = 2x2 + 8x + 6 y0 = 1 y0′ = 4

Предыдущий
вопрос Следующий
вопрос

Данное уравнение является дифференциальным уравнением второго порядка, поэтому для использования различных методов необходимо свести его к системе уравнений первого порядка.

Заменим y' = z, тогда получим систему:
y' = z,
z' = 2x^2 + 8x + 6 - 2z - 2y.

Зададим начальные условия:
y(0) = 1,
z(0) = 4.

Определим функцию, задающую правую часть системы уравнений:

function dz = myode(x, z)
dz = [z(2); 2x^2 + 8x + 6 - 2z(2) - 2z(1)];
end

Теперь решим данную систему численно с помощью метода Эйлера и метода Рунге-Кутта:

% Метод Эйлера
h = 0.01;
x = 0:h:1;
z0 = [1; 4];
Z = zeros(2, length(x));
Z(:, 1) = z0;
for i = 1:(length(x)-1)
Z(:, i+1) = Z(:, i) + h*myode(x(i), Z(:, i));
end

% Метод Рунге-Кутта
h = 0.01;
x = 0:h:1;
z0 = [1; 4];
Z_rk = zeros(2, length(x));
Z_rk(:, 1) = z0;
for i = 1:(length(x)-1)
k1 = hmyode(x(i), Z_rk(:, i));
k2 = hmyode(x(i) + 0.5h, Z_rk(:, i) + 0.5k1);
k3 = hmyode(x(i) + 0.5h, Z_rk(:, i) + 0.5k2);
k4 = hmyode(x(i) + h, Z_rk(:, i) + k3);
Z_rk(:, i+1) = Z_rk(:, i) + (k1 + 2k2 + 2k3 + k4)/6;
end

Далее решим данную систему с помощью функции ode45 из пакета MatLab:

h = 0.01;
[x_ode45, Z_ode45] = ode45(@myode, [0 1], [1, 4]);

Наконец, построим графики решений на одном графике:

plot(x, Z(1, :), 'r', 'LineWidth', 1.5);
hold on;
plot(x, Z_rk(1, :), 'g', 'LineWidth', 1.5);
plot(x_ode45, Z_ode45(:, 1), 'b', 'LineWidth', 1.5);
legend('Метод Эйлера', 'Метод Рунге-Кутта', 'ode45');
xlabel('x');
ylabel('y(x)');
title('Решения дифференциального уравнения y'''' + 2y'' + 2y = 2x^2 + 8x + 6');
grid on;
hold off;