四阶龙格-库塔法求微分方程组matlab程序

龙格-库塔法是求解微分方程数值解的一种常用方法。四阶龙格-库塔法是其中的一种，它的步骤比较多。下面，我会提供一个MATLAB程序，用来实现四阶龙格-库塔法求解微分方程组。

1. 首先，我们需要定义一些变量和初始值，包括：

- h: 步长，即每次迭代的步幅；
- t0: 初始时间；
- tn: 结束时间；
- y0: 初始状态向量，即微分方程的初值。

代码：

h = 0.01;
t0 = 0;
tn = 10;
y0 = [1;0;0;1];

2. 接下来，我们需要定义一个求导函数，即微分方程的右侧。在这个例子中，我们使用的是下面这个微分方程：

y1' = y2
y2' = -y1
y3' = y4
y4' = -y3

代码：

function dydt = derivative(t,y)
dydt = [y(2); -y(1); y(4); -y(3)];

3. 然后，我们就可以开始用四阶龙格-库塔法求解微分方程组了。具体步骤如下：

- 定义初始值；
- 定义一个空的数组来存储解；
- 定义一个循环，每次迭代计算下一个解；
- 在每次迭代中，使用四阶龙格-库塔法计算下一个解，并将其存储到数组中。

代码：

% Define initial values
t = t0;
y = y0;
i = 1;
results(i,:) = y';

% Loop over time steps
while t < tn
% Calculate next solution using RK4
k1 = h * derivative(t,y);
k2 = h * derivative(t+h/2, y+k1/2);
k3 = h * derivative(t+h/2, y+k2/2);
k4 = h * derivative(t+h, y+k3);

y = y + (k1 + 2*k2 + 2*k3 + k4) / 6;
t = t + h;
i = i + 1;
results(i,:) = y';
end

4. 最后，我们可以画出结果。在这个例子中，我们只有4个状态变量，因此可以在4个子图中分别画出它们的时间演化曲线。

代码：

% Plot results
figure;
for i = 1:4
subplot(2,2,i);
plot(0:h:tn, results(:,i));
title(['y', num2str(i)]);
end

到此为止，我们就完成了一个简单的四阶龙格-库塔法求解微分方程组的MATLAB程序。