matlab龙格库塔法求解微分方程代码

Matlab中使用龙格-库塔（Runge-Kutta）方法求解微分方程是一个常见的数值计算任务。以下是一个简单的4阶龙格-库塔（RK4）算法的示例代码，用于求解一阶线性常微分方程：

function [t, y] = rk4(f, t0, y0, tfinal, h)
% f: 微分方程的函数 handle (f(t, y) -> dy/dt)
% t0: 初始时间
% y0: 初始值
% tfinal: 结束时间
% h: 时间步长

% 初始化
nsteps = floor((tfinal - t0) / h);
t = zeros(1, nsteps + 1);   % 存储时间点
y = zeros(size(y0));        % 存储每个时间步的解

% 龙格-库塔方法的核心步骤
for i = 1:nsteps
    k1 = h * f(t(i), y);
    k2 = h * f(t(i) + h/2, y + k1/2);
    k3 = h * f(t(i) + h/2, y + k2/2);
    k4 = h * f(t(i) + h, y + k3);

    % 更新解
    y = y + (k1 + 2*k2 + 2*k3 + k4)/6;

    % 记录时间
    t(i+1) = t(i) + h;
end
t(nsteps+1) = tfinal;      % 设置最后一个时间点

这个函数`rk4`接受一个描述微分方程导数的匿名函数`f`作为输入，你需要自定义这个函数来表达你的具体方程。例如，如果你有一个方程`dy/dt = f(t, y)`，那么`f`会像这样定义：`f = @(t, y) ...;`。

你可以调用这个函数并传入你的方程、初始条件和参数，如下所示：

% 假设你有一个二阶线性方程 dy'' = -k*y
k = 1;       % 权重系数
f = @(t, y) [0; -k*y(1)];

% 示例使用
tspan = [0 10];     % 时间范围
y0 = [1; 0];       % 初始状态
[t, y] = rk4(f, tspan(1), y0, tspan(2), 0.1);

运行后，`t`数组将包含时间点，`y`数组则包含了对应时间下的解。对于更复杂的微分方程或特定的需求，可能还需要根据实际情况进行适当修改或扩展。