四阶龙格库塔法的计算例子

序

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

一、矩形法

1.1 原理

设微分方程

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

求yyy。

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

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

易得

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

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

1.2 例子

以
y˙=e−yy0=0(1.4)\begin{aligned} \dot y &= e^{-y} \\ y_0 &= 0 \tag{1.4} \end{aligned}y˙y0=e−y=0(1.4)

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

%% 不同步长下的矩形法比较dt1 = 0.1;                % 步长1
dt2 = 1.0;                 % 步长2t1 = 0:dt1:10;            % 时间1
t2 = 0:dt2:10;             % 时间2y1 = 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 - 1dt = t(n+1) - t(n);             % 计算步长 dty(n+1) = y(n) + f(y(n)) * dt;   % 累加计算 yend
end

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

二、龙格库塔法

2.1 y˙=f(y)\dot y=f(y)y˙=f(y) 形式

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

yn+1=yn+h6(k1+2k2+2k3+k4)k1=f(yn)k2=f(yn+k1h2)k3=f(yn+k2h2)k4=f(yn+k3h)(2.1)\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}yn+1k1k2k3k4=yn+6h(k1+2k2+2k3+k4)=f(yn)=f(yn+k12h)=f(yn+k22h)=f(yn+k3h)(2.1)

此处，采用习惯上的符号hhh，上面的例子中，h=Δth=\Delta th=Δ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 - 1dt = t(n+1) - t(n);             % 计算步长 dty(n+1) = y(n) + f(y(n)) * dt;   % 累加计算 yend
end%% 四阶龙格库塔法
function y = ode_rk(t, y0)N = length(t);y = zeros(N,1);y(1) = y0;for n = 1:N-1h = t(n+1) - t(n);              % 步长（即时间间隔）k1 = f(y(n));                   % k1k2 = f(y(n) + h/2*k1);          % k2k3 = f(y(n) + h/2*k2);          % k3k4 = f(y(n) + h*k3);            % k4y(n+1) = y(n) + h/6 * (k1 + 2*k2 + 2*k3 + k4);  % 累加计算yend
end

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

2.2 y˙=f(t)\dot y=f(t)y˙=f(t) 形式

设

y˙=f(t)(2.2)\dot y = f(t) \tag{2.2} y˙=f(t)(2.2)

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

yn+1=yn+h6(k1+2k2+2k3+k4)k1=f(tn)k2=f(tn+h2)k3=f(tn+h2)k4=f(tn+h)(2.3)\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}yn+1k1k2k3k4=yn+6h(k1+2k2+2k3+k4)=f(tn)=f(tn+2h)=f(tn+2h)=f(tn+h)(2.3)

一个简单的例子是

y˙=sin⁡(t)y0=0(2.4)\begin{aligned} &\dot y = \sin(t) \\ &y_0 = 0 \end{aligned} \tag{2.4}y˙=sin(t)y0=0(2.4)

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

%% 龙格库塔法步长差异比较
dt1 = 1;                % 步长为1
dt2 = 0.1;              % 步长为0.1
T = 10;                 % 总时间10st1 = 0:dt1:T;           % 时间t1
y1 = ode_rk(t1, 0);     % 解y1
t2 = 0:dt2:T;           % 时间t2
y2 = ode_rk(t2, 0);     % 解y2figure(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-1h = t(n+1) - t(n);              % 步长（即时间间隔）k1 = f(t(n));                   % k1k2 = f(t(n) + h/2);          % k2k3 = f(t(n) + h/2);          % k3k4 = f(t(n) + h);            % k4y(n+1) = y(n) + h/6 * (k1 + 2*k2 + 2*k3 + k4);  % 累加计算yend
end

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

2.3 y˙=f(y,t)\dot y=f(y, t)y˙=f(y,t) 形式

设
y˙=f(y,t)\dot y=f(y, t)y˙=f(y,t)

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

yn+1=yn+h6(k1+2k2+2k3+k4)k1=f(yn,tn)k2=f(yn+k1h2,tn+h2)k3=f(yn+k2h2,tn+h2)k4=f(yn+k3h,tn+h)(2.5)\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}yn+1k1k2k3k4=yn+6h(k1+2k2+2k3+k4)=f(yn,tn)=f(yn+k12h,tn+2h)=f(yn+k22h,tn+2h)=f(yn+k3h,tn+h)(2.5)

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