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

序

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

一、矩形法

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 例子

为例，时间取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)$ 形式

求解$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)}$$

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