---

序

没有对比就没有伤害，本文先给出很多时候直接采用的矩形法，然后与四阶龙格库塔法做比较，着重说明四阶龙格库塔法。

---

一、矩形法

---

1.1 原理

设微分方程

$$\dot{y}=f(y)\text{\hspace{1em}(1.1)}$$

求$y$。

使用数值方法，离散化得每一步的增量

$$\Delta y=f(y)\Delta t\text{\hspace{1em}(1.2)}$$

易得

$$y_{n+1}=y_n+f(y_n)\Delta t\text{\hspace{1em}(1.3)}$$

实际上，这就是矩形法计算积分。当 $\Delta t\to 0$时，可以得出很高精度的$y_n$，但实际工程中未必能够取很小的$\Delta t$。

---

1.2 例子

以
$$\begin{array}{rl}\dot{y}& =e^{-y}\\ y_0& =0\end{array}\text{\hspace{1em}(1.4)}$$

为例，时间取1~10s，分别取 $\Delta t=0.1,\Delta t=1$，查看不同精度下的运算结果。式(1.4)可求出解析解为$y=\ln (t+1)$，用于比较求解精度。

%% 不同步长下的矩形法比较

dt1 = 0.1;              	% 步长1
dt2 = 1.0;              	% 步长2

t1 = 0:dt1:10;          	% 时间1
t2 = 0:dt2:10;          	% 时间2
    
y1 = ode_rect(t1, 0);      	% 精度1下计算结果
y2 = ode_rect(t2, 0);      	% 精度2下计算结果

plot(t1,log(t1+1),t1,y1,t2,y2);
legend('理论值', 'dt=0.1', 'dt=1');
grid on;grid minor;xlabel 't';ylabel 'y'

%% 导数方程
function dy=f(y)
    dy = exp(-y);
end

%% 矩形法
function y = ode_rect(t, y0)
    N = length(t);
    y = zeros(N,1);
    y(1) = y0;
    for n = 1:N - 1
        dt = t(n+1) - t(n);             % 计算步长 dt
        y(n+1) = y(n) + f(y(n)) * dt;   % 累加计算 y
    end
end

可见，当步长为0.1时，矩形法的精度较高，但步长为1时，矩形法误差大。

---

二、龙格库塔法

---

2.1 $\dot{y}=f(y)$ 形式

经过多年潜心研究，龙格库塔站在前人的肩膀上，发现了一种高精度的方法。那就是把式(1.3)的计算改为

$$\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(y_n)\\ k_2& =f(y_n+k_1\frac{h}{2})\\ k_3& =f(y_n+k_2\frac{h}{2})\\ k_4& =f(y_n+k_{3}h)\end{array}\text{\hspace{1em}(2.1)}$$

此处，采用习惯上的符号$h$，上面的例子中，$h=\Delta t$。

依旧对于式(1.4)，步长取1，分别使用矩形法和四阶龙格库塔法求解，结果如下

%% 矩形法与龙格库塔法比较
dt = 1.0;
t = 0:dt:10;
y1 = ode_rect(t, 0);       % 矩形法计算    
y2 = ode_rk(t, 0);         % 龙格库塔法计算

plot(0:0.01:10,log([0:0.01:10]+1),t,y1,t,y2);
legend('理论值', '矩形法', '龙格库塔法');
grid on;grid minor;xlabel 't';ylabel 'y'

%% 导数方程
function dy=f(y)
    dy = exp(-y);
end

%% 矩形法
function y = ode_rect(t, y0)
    N = length(t);
    y = zeros(N,1);
    y(1) = y0;
    for n = 1:N - 1
        dt = t(n+1) - t(n);             % 计算步长 dt
        y(n+1) = y(n) + f(y(n)) * dt;   % 累加计算 y
    end
end

%% 四阶龙格库塔法
function y = ode_rk(t, y0)
    N = length(t);
    y = zeros(N,1);
    y(1) = y0;
    for n = 1:N-1
        h = t(n+1) - t(n);              % 步长（即时间间隔）
        k1 = f(y(n));                   % k1
        k2 = f(y(n) + h/2*k1);          % k2
        k3 = f(y(n) + h/2*k2);          % k3
        k4 = f(y(n) + h*k3);            % k4
        y(n+1) = y(n) + h/6 * (k1 + 2*k2 + 2*k3 + k4);  % 累加计算y
    end
end

可见，四阶龙格库塔法很好地接近真实值。

---

2.2 $\dot{y}=f(t)$ 形式

设

$$\dot{y}=f(t)\text{\hspace{1em}(2.2)}$$

从理论计算微分方程的角度，式(1.1) 和式(2.2)有着截然不同的求解方式。但是使用数值方法，只不过是把 $y$变为$t$。四阶龙格库塔法求解式(2.2)的方法如下

$$\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(t_n)\\ k_2& =f(t_n+\frac{h}{2})\\ k_3& =f(t_n+\frac{h}{2})\\ k_4& =f(t_n+h)\end{array}\text{\hspace{1em}(2.3)}$$

一个简单的例子是

$$\begin{array}{rl}& \dot{y}=\sin (t)\\ & y_0=0\end{array}\text{\hspace{1em}(2.4)}$$

其解析解为$y=1-\cos (t)$。设步长分别为1,0.1，使用四阶龙格库塔法发求解如下

%% 龙格库塔法步长差异比较
dt1 = 1;                % 步长为1
dt2 = 0.1;              % 步长为0.1
T = 10;                 % 总时间10s

t1 = 0:dt1:T;           % 时间t1
y1 = ode_rk(t1, 0);     % 解y1
t2 = 0:dt2:T;           % 时间t2
y2 = ode_rk(t2, 0);     % 解y2

figure(1)
plot(0:0.01:T, 1-cos(0:0.01:T));hold on                 % 解析解计算值
plot(t1,y1, '-o', t2,y2, '--');hold off                             % 数值解计算值
legend('理论值','dt=1', 'dt=0.1');title('龙格库塔法');
grid on;grid minor;xlabel 't';ylabel 'y'

%% 导数方程
function dy=f(t)
    dy = sin(t);
end

%% 四阶龙格库塔法
function y = ode_rk(t, y0)
    N = length(t);
    y = zeros(N,1);
    y(1) = y0;
    for n = 1:N-1
        h = t(n+1) - t(n);              % 步长（即时间间隔）
        k1 = f(t(n));                   % k1
        k2 = f(t(n) + h/2);          % k2
        k3 = f(t(n) + h/2);          % k3
        k4 = f(t(n) + h);            % k4
        y(n+1) = y(n) + h/6 * (k1 + 2*k2 + 2*k3 + k4);  % 累加计算y
    end
end

从例子中可见，步长为1时，龙格库塔法依旧得到精确的结果。

---

2.3 $\dot{y}=f(y,t)$ 形式

设
$$\dot{y}=f(y,t)$$

求解$y$。不过是多了一个变量，使用四阶龙格库塔法计算方法为

$$\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(y_n,t_n)\\ k_2& =f(y_n+k_1\frac{h}{2},t_n+\frac{h}{2})\\ k_3& =f(y_n+k_2\frac{h}{2},t_n+\frac{h}{2})\\ k_4& =f(y_n+k_{3}h,t_n+h)\end{array}\text{\hspace{1em}(2.5)}$$

更多变量，以此类推，不再赘述。