四阶龙格库塔法（Matlab）

贡献者： addis

预备知识　中点法解常微分方程（组）

　　¹龙格库塔法（Runge-Kutta）是一类数值解微分方程或方程组的算法，其中较常见的是四阶龙格库塔法。这里不进行推导，仅仅给出算法。如果你不想了解算法只想用 Matlab 自带的解算器直接求解方程，见 “Matlab 解常微分方程组（ode45）”。

　　公式如下（$y_n, t_n, h$ 的定义类比式 5 ）

\begin{equation} y_{n+1} = y_n + \frac h6 (k_1 + 2k_2 + 2k_3 + k_4)~, \end{equation}

其中

\begin{equation} \begin{aligned} k_1 &= f(y_n, t_n)~, & k_2 &= f \left(y_n + h\frac{k_1}{2}, t_n + \frac h2 \right) ~,\\ k_3 &= f \left(y_n + h\frac{k_2}{2}, t_n + \frac h2 \right) ~, \qquad &k_4 &= f(y_n + hk_3, t_n + h)~. \end{aligned} \end{equation}

由以上两式，不难把该算法拓展到方程组的情况。对于 $N$ 元微分方程组

\begin{equation} \begin{cases} y'_1(t) = f_1(y_1,\dots, y_N, t)\\ y'_2(t) = f_2(y_1,\dots, y_N, t)\\ \qquad\;\; \vdots\\ y'_N(t) = f_N(y_1,\dots, y_N, t)~. \end{cases} \end{equation}

我们可以把该式记为矢量函数的形式

\begin{equation} \boldsymbol{\mathbf{y}} '(t) = \boldsymbol{\mathbf{f}} ( \boldsymbol{\mathbf{y}} , t)~, \end{equation}

其中

\begin{equation} \boldsymbol{\mathbf{y}} '(t) = \begin{pmatrix}y'_1(t)\\ \vdots\\ y'_N(t)\end{pmatrix} ~, \qquad \boldsymbol{\mathbf{y}} = \begin{pmatrix}y_1\\ \vdots\\ y_N\end{pmatrix} ~. \end{equation}

现在我们仅需要把式 1 和式 2 中的所有 $y_i$ 和 $k_i$ 都变为 $N$ 维列矢量 $ \boldsymbol{\mathbf{y}} _i$ 和 $ \boldsymbol{\mathbf{k}} _i$ 即可将微分方程拓展为微分方程组。

1. 例程：天体运动

　　到目前为止，我们每求一个微分方程的数值解都要重新写一次程序，对于一些较为复杂的算法这样做效率较低。我们这里不妨把四阶龙格库塔法写到一个单独的函数文件 odeRK4.m 中，当我们要解某个特定的方程时，只需把式 2 中的 $f(y, t)$ 作为自变量输入即可解出 $y(t)$。这样的程序就叫做微分方程的解算器，它可以求解用户自定义的微分方程。

代码 1：odeRK4.m

% 四阶龙格库塔定步长节微分方程
% f(Y, t): 求导函数
% tspan: 二元向量，起始和终止时间
% Y0: 初值（列向量）
% Nt: 时间节点数
function [Y, t] = odeRK4(f, tspan, Y0, Nt)
Nvar = numel(Y0);  % 因变量的个数
dt = (tspan(2) - tspan(1)) / (Nt-1); % 计算步长
Y = zeros(Nvar, Nt); % 预赋值
Y(:, 1) = Y0(:); % 初值
t = linspace(tspan(1), tspan(2), Nt);

for ii=1:Nt-1
    K1 = f(Y(:,ii)          , t(ii)      );
    K2 = f(Y(:,ii)+K1*dt/2  , t(ii)+dt/2 );
    K3 = f(Y(:,ii)+K2*dt/2  , t(ii)+dt/2 );
    K4 = f(Y(:,ii)+K3*dt    , t(ii)+dt   );
    Y(:,ii+1) = Y(:,ii) + dt/6 * (K1+2*K2+2*K3+K4);
end
end

　　我们先来看第 1 行的函数声明，输入变量中，f 是式 4 中 $ \boldsymbol{\mathbf{f}} ( \boldsymbol{\mathbf{y}} , t)$ 的函数句柄，tspan 是一个 2×1 的列矢量，tspan(1) 是初始时间，tspan(2) 是终止时间，Y0 是一个列矢量，Y0(ii) 是第 ii 个因变量的初始值，Nt 是 $t_n$ 的个数，tspan 定义的时间区间被等分为 Nt - 1 个小区间。因变量中，Y 的行数是因变量的个数，列数是 Nt，t 是一个行矢量，由第 6 行定义，Y(ii, jj) 就是第 ii 个变量在 t(jj) 时刻的值。第 5 行把初值 Y0 赋给 Y 的第 1 列，第 8-14 行的循环根据式 1 和式 2 的矢量形式由 Y 的第 ii 列（$ \boldsymbol{\mathbf{y}} _i$）求第 ii+1 列（$ \boldsymbol{\mathbf{y}} _{i+1}$）。

　　我们先来用这个函数来计算 “天体运动的简单数值计算” 中的问题。我们令因变量 $ \boldsymbol{\mathbf{y}} $ 的四个分量依次为一阶方程组（式 4 ）

\begin{equation} \begin{cases} x' = v_x\\ y' = v_y\\ v'_x = -GMx/(x^2 + y^2)^{3/2}\\ v'_y = -GMy/(x^2 + y^2)^{3/2} \end{cases}~ \end{equation}

中的 $x, y, v_x, v_y$。程序代码如下

代码 2：keplerRK4.m

function keplerRK4
% 参数设定
GM = 1; % 万有引力常数乘以中心天体质量
x0 = 1; y0 = 0; % 初始位置
vx0 = 0; vy0 = 0.7; % 初始速度
tspan = [0; 4]; % 总时间和步数
Nt = 100; % 步数

Y0 = [x0; y0; vx0; vy0]; % 因变量初值
f = @(Y, t)fun(Y, t, GM);
[Y,~] = odeRK4(f, tspan, Y0, Nt);

% 画图
figure; hold on;
plot(Y(1,:), Y(2,:));
scatter(0, 0);
axis equal;
end

function Y1 = fun(Y, ~, GM)
% 因变量
x = Y(1); y = Y(2);
vx = Y(3); vy = Y(4);
Y1 = zeros(4,1); % 预赋值
Y1(1) = vx;
Y1(2) = vy;
temp = -GM /(x^2 + y^2)^(3/2);
Y1(3) = temp * x;
Y1(4) = temp * y;
end

　　运行结果如图 1 所示。

图 1：运行结果

　　我们先来看函数 fun（20 行），这个函数就相当于式 6 。第一个输入变量 Y 是一个列矢量，是 $ \boldsymbol{\mathbf{y}} _n$ 的值，第二个输入变量是 $t$，但由于式 6 中没有出现 $t$，我们用波浪线代替。第三个输入变量是参数 $GM$，即万有引力常数和中心天体质量之积。输出变量 Y1 是一个列矢量，是 $ \boldsymbol{\mathbf{y}} '_n$ 的值。

　　再来看主函数 KeplerRK4（第 1 行），参数设定中除了步数从 4000 变为了 100，其他都和 “天体运动的简单数值计算” 中的程序一样，然而这里运行结果却精确得多（曲线几乎闭合），可见这种算法的优越性。

　　主函数第 10 行中将 fun(Y, t, GM) 变为函数句柄 f(Y, t)，这样 GM 就可以在 “参数设定” 中设置，而不用在 fun 函数内部设置。第 11 行调用了上文中的 odeRK4 函数解方程组，由于我们在画图时不需要用 t，所以把第二个输出变量改为波浪线。

1. ^ 参考 Wikipedia 相关页面。

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。