python 4阶龙哥库塔 四阶龙格库塔法计算题

---

序

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

---

一、矩形法

---

1.1 原理

设微分方程

$\dot y=f(y) \tag{1.1}$

求$y$。

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

$\Delta y = f(y)\Delta t \tag{1.2}$

易得

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

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

---

1.2 例子

以
$\begin{aligned} \dot y &= e^{-y} \\ y_0 &= 0 \tag{1.4} \end{aligned}$

为例，时间取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{aligned} 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{aligned} \tag{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) \tag{2.2}$

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

$\begin{aligned} 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{aligned} \tag{2.3}$

一个简单的例子是

$\begin{aligned} &\dot y = \sin(t) \\ &y_0 = 0 \end{aligned} \tag{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{aligned} 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{aligned} \tag{2.5}$

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