经典四阶龙格-库塔方法（Runge-Kutta Method）的 MATLAB 实现

经典四阶龙格-库塔方法是一种求解常微分方程的数值解法，具有较高的精度和稳定性。以下是一个基于 MATLAB 的四阶龙格-库塔方法的实现示例。

1. 算法原理

对于常微分方程：
\(\frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0\)
四阶龙格-库塔方法通过以下步骤计算 ( y ) 在 ( t ) 的下一个值：

1. \(k_1 = h \cdot f(t_n, y_n)\)
2. \(k_2 = h \cdot f(t_n + \frac{h}{2}, y_n + \frac{k_1}{2})\)
3. \(k_3 = h \cdot f(t_n + \frac{h}{2}, y_n + \frac{k_2}{2})\)
4. \(k_4 = h \cdot f(t_n + h, y_n + k_3)\)
5. \(y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)\)

2. MATLAB 实现

function [t, y] = runge_kutta_4th_order(ode_func, t0, y0, tf, h)
    % 输入：
    %   ode_func - 常微分方程的函数句柄
    %   t0 - 初始时间
    %   y0 - 初始条件
    %   tf - 终止时间
    %   h - 步长
    % 输出：
    %   t - 时间向量
    %   y - 解向量

    t = t0;
    y = y0;
    t_values = [t0];
    y_values = [y0];

    while t < tf
        k1 = h * ode_func(t, y);
        k2 = h * ode_func(t + h/2, y + k1/2);
        k3 = h * ode_func(t + h/2, y + k2/2);
        k4 = h * ode_func(t + h, y + k3);

        y = y + (1/6) * (k1 + 2*k2 + 2*k3 + k4);
        t = t + h;

        t_values = [t_values, t];
        y_values = [y_values, y];
    end

    t = t_values;
    y = y_values;
end

3. 使用示例

假设我们有一个常微分方程：
\(\frac{dy}{dt} = -2y + 4e^{-t}\)
初始条件为 ( y(0) = 1 )，我们希望在 ( t = 2 ) 时得到解，步长为 0.1。

% 定义常微分方程
ode_func = @(t, y) -2*y + 4*exp(-t);

% 初始条件
t0 = 0;
y0 = 1;
tf = 2;
h = 0.1;

% 调用四阶龙格-库塔方法
[t, y] = runge_kutta_4th_order(ode_func, t0, y0, tf, h);

% 绘制结果
figure;
plot(t, y, 'b-');
xlabel('Time');
ylabel('Solution y(t)');
title('Solution of ODE using 4th Order Runge-Kutta Method');

总结

MATLAB 代码实现了经典四阶龙格-库塔方法，可以用于求解各种常微分方程。通过调整步长和初始条件，可以得到不同精度的数值解。

参考代码 根据经典四阶龙格库塔算法方程公式编写的龙格库塔程序