贡献者: addis
1龙格库塔法(Runge-Kutta)是一类数值解微分方程或方程组的算法,其中较常见的是四阶龙格库塔法。这里不进行推导,仅仅给出算法。如果你不想了解算法只想用 Matlab 自带的解算器直接求解方程,见 “Matlab 解常微分方程组(ode45)”。
公式如下($y_n, t_n, h$ 的定义类比式 5 )
到目前为止,我们每求一个微分方程的数值解都要重新写一次程序,对于一些较为复杂的算法这样做效率较低。我们这里不妨把四阶龙格库塔法写到一个单独的函数文件 odeRK4.m
中,当我们要解某个特定的方程时,只需把式 2 中的 $f(y, t)$ 作为自变量输入即可解出 $y(t)$。这样的程序就叫做微分方程的解算器,它可以求解用户自定义的微分方程。
% 四阶龙格库塔定步长节微分方程
% 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 )
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 所示。
我们先来看函数 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
,所以把第二个输出变量改为波浪线。