写一段Matlab代码,实现四阶龙格库塔来解二元常微分方程
时间: 2024-02-28 22:54:40 浏览: 25
好的,下面是一个简单的实现:
```matlab
function [t, y] = runge_kutta_4(f, y0, tspan, h)
% 四阶龙格库塔法求解二元常微分方程
% f: 函数句柄,表示 dy/dt = f(t, y)
% y0: 初始状态,一个包含两个元素的向量
% tspan: 时间区间,一个包含两个元素的向量
% h: 步长
t0 = tspan(1);
tf = tspan(2);
t = t0:h:tf;
n = length(t);
y = zeros(2, n);
y(:, 1) = y0;
for i = 1:(n-1)
k1 = h * f(t(i), y(:, i));
k2 = h * f(t(i) + h/2, y(:, i) + k1/2);
k3 = h * f(t(i) + h/2, y(:, i) + k2/2);
k4 = h * f(t(i) + h, y(:, i) + k3);
y(:, i+1) = y(:, i) + (k1 + 2*k2 + 2*k3 + k4)/6;
end
end
```
其中,`f` 是一个函数句柄,表示二元常微分方程 `dy/dt = f(t, y)` 中的右侧函数,`y0` 是初始状态,一个包含两个元素的向量,`tspan` 是时间区间,一个包含两个元素的向量,`h` 是步长。函数返回的 `t` 是时间向量,`y` 是状态矩阵,其中每一列是一个时刻的状态。
相关问题
写一段Matlab代码,实现四阶龙格库塔来解三元微分方程
好的,下面是一个简单的实现:
```matlab
function [t, y] = runge_kutta_4(f, y0, tspan, h)
% 四阶龙格库塔法求解三元常微分方程
% f: 函数句柄,表示 dY/dt = f(t, Y)
% y0: 初始状态,一个包含三个元素的向量
% tspan: 时间区间,一个包含两个元素的向量
% h: 步长
t0 = tspan(1);
tf = tspan(2);
t = t0:h:tf;
n = length(t);
y = zeros(3, n);
y(:, 1) = y0;
for i = 1:(n-1)
k1 = h * f(t(i), y(:, i));
k2 = h * f(t(i) + h/2, y(:, i) + k1/2);
k3 = h * f(t(i) + h/2, y(:, i) + k2/2);
k4 = h * f(t(i) + h, y(:, i) + k3);
y(:, i+1) = y(:, i) + (k1 + 2*k2 + 2*k3 + k4)/6;
end
end
```
其中,`f` 是一个函数句柄,表示三元常微分方程 `dY/dt = f(t, Y)` 中的右侧函数,`y0` 是初始状态,一个包含三个元素的向量,`tspan` 是时间区间,一个包含两个元素的向量,`h` 是步长。函数返回的 `t` 是时间向量,`y` 是状态矩阵,其中每一列是一个时刻的状态。
写一段Matlab代码,实现四阶龙格库塔来解三元微分方程,并画出原函数图和解函数图
好的,下面是一个简单的实现,以一个三元常微分方程为例:
$$
\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` 函数。最后,我们画出了原函数图和解函数图进行比较。