写一段Matlab代码,实现四阶龙格库塔来解三元微分方程,并画出原函数图和解函数图
时间: 2024-02-28 08:54:47 浏览: 63
好的,下面是一个简单的实现,以一个三元常微分方程为例:
$$
\begin{cases}
y'_1 = y_2 \\
y'_2 = y_3 \\
y'_3 = -2y_3 - 3y_2 - 4y_1 + \sin(t)
\end{cases}
$$
```matlab
% 定义右侧函数
f = @(t, y) [y(2); y(3); -2*y(3) - 3*y(2) - 4*y(1) + sin(t)];
% 定义初始状态和时间区间
y0 = [0; 0; 0];
tspan = [0, 10];
% 求解
h = 0.01;
[t, y] = runge_kutta_4(f, y0, tspan, h);
% 画图比较
figure;
hold on;
% 画原函数图
t_orig = linspace(tspan(1), tspan(2), 1000);
y_orig = zeros(3, length(t_orig));
for i = 1:length(t_orig)
y_orig(:, i) = [sin(t_orig(i)); cos(t_orig(i)); -2*cos(t_orig(i)) - 3*sin(t_orig(i)) - 4*sin(t_orig(i))];
end
plot(t_orig, y_orig(1, :), 'r', 'LineWidth', 2);
plot(t_orig, y_orig(2, :), 'g', 'LineWidth', 2);
plot(t_orig, y_orig(3, :), 'b', 'LineWidth', 2);
% 画解函数图
plot(t, y(1, :), 'r--', 'LineWidth', 1);
plot(t, y(2, :), 'g--', 'LineWidth', 1);
plot(t, y(3, :), 'b--', 'LineWidth', 1);
legend('y_1(t)', 'y_2(t)', 'y_3(t)', 'y_{1, approx}(t)', 'y_{2, approx}(t)', 'y_{3, approx}(t)');
xlabel('t');
title('Solution of the ODE using Runge-Kutta 4');
```
其中,右侧函数 `f` 表示三元常微分方程的右侧函数,`y0` 表示初始状态,`tspan` 表示时间区间。求解过程使用了之前实现的 `runge_kutta_4` 函数。最后,我们画出了原函数图和解函数图进行比较。
阅读全文