开普勒问题的数值计算（Matlab）

贡献者： addis

预备知识　Matlab 画图，开普勒问题的运动方程

　　本文讨论开普勒问题中，若已知初始位置和初速度矢量，求轨道方程和动画的算法并给出 Matlab 代码。以下所有符号和定义沿用 “开普勒问题”。

图 1：在同一位置以同一初速度不同仰角发射的质点的椭圆轨道，具有相同的半长轴和不同的离心率，动画见这里。虚线为参考圆，例如地球表面，可见除了水平发射的质点都会落回地面。若初速度很小，那么参考圆外的运动就趋近于抛体运动（足够扁的椭圆的一端趋近于抛物线）。

图 2：在同一位置以同

45^{\circ}

仰角和不同初速度发射的质点的轨道，动画见这里。速度最大的轨道为双曲线，其他轨道为椭圆。

图 3：平方反比斥力产生的双曲线轨道，动画见这里。这就是著名的卢瑟福散射。

1. 轨道的计算

　　若 $F (r)$ 是平方反比的力（斥力为正引力为负），即

\begin{matrix} (1) & F (r) = \frac{k}{r^{2}} . \end{matrix}

其中

k

为非零实常数。由于质点的运动只和初始位置、初始速度和

k

有关，与质点的质量无关，所以我们令

m = 1

并不影响一般性。令发射位置的极坐标为

(R, α)

，初速度为

v_{0}

，发射仰角为

β

（顺时针增量为正，

β = 0

时质点延逆时针方向水平发射）。

　　所以质点机械能为动能加势能 $E = v_{0}^{2} / 2 + k / R$ ，角动量（逆时针为正）为 $L = v_{0} R \cos β$ 。于是离心率 $e$ 可由式 2 计算，半通径 $p$ 由式 3 计算。然后代入圆锥曲线的极坐标方程（式 1 ）即可得到轨道方程。然而由于发射点已经选好，实际的轨道还需要旋转一个角度 $θ_{0}$ （逆时针为正），即轨道方程变为

\begin{matrix} (2) & r (θ) = \frac{- s p}{1 + s e \cos (θ - θ_{0})} (r > 0) . \end{matrix}

其中当

k > 0

时

s = 1

，

k < 0

时

s = - 1

。要求

θ_{0}

，可以代入初始条件，即

r (α) = R

。得

\begin{matrix} (3) & θ_{0} = α \pm \arccos (\frac{- s - p / R}{e}) (e \neq 0) . \end{matrix}

这里有两个解是因为圆锥曲线和半径为

R

的圆必定有两个交点。我们需要根据发射角度来判断哪个是发射点，稍作分析易得当

\sin 2 β ⩾ 0

时上式取正号，否则取负号。

2. 运动的计算

　　我们使用 “开普勒问题的运动方程” 给出的方法求解三种圆锥曲线的含时运动。沿用其中的符号定义，令 $t = 0$ 时质点位于近日点，近日点对应极角 $θ_{0} + π$ ，令相对近日点的极角增量为 $Δ θ = θ - θ_{0} - π$ 。接下来，先计算函数 $t (Δ θ)$ ，默认逆时针旋转，然后数值求解反函数 $Δ θ (t)$ 即可。若要考虑顺时针旋转，函数值乘以 $- 1$ 即可。易得 $\cos β > 0$ 时质点绕中心天体逆时针运动， $\cos β < 0$ 时顺时针运动。

　　 抛物线：时间 $t$ 在 $(- \infty, \infty)$ 取值，式 4 就是 $t (Δ θ)$ 。

　　椭圆：时间 $t$ 在 $[- T / 2, T / 2]$ 取值，周期 $T$ 由式 5 计算。由于这是周期运动，其他任意时间都可以加减整数个周期后落到 $[- T / 2, T / 2]$ 中。由式 5 和式 6 计算 $t (Δ θ)$ 。

　　 双曲线：时间 $t$ 在 $(- \infty, \infty)$ 取值，使用式 7 到式 11 可得 $k < 0$ 和 $k > 0$ 两种情况时的 $t (Δ θ)$ 。

　　求出发射质点的初始时刻为 $t_{0} = t (α - θ_{0} - π)$ ，然后生成等间距时间格点 $t_{i}$ （ $i = 0, 1, 2, \dots$ ），用二分法解方程得到对应的 $θ_{i}$ 即可。

　　以下 Matlab 程序支持对同一发射位置画出多组不同初速度和仰角的轨道。相比于使用解微分方程的方法（子节 1 ）计算轨道和运动，该程序计算的轨道和运动都可以达到 13 位有效数字以上。然而缺点是只能计算理想的开普勒问题，无法拓展到非平方反比力，也无法拓展到三体乃至多体运动。程序会在当前目录中生成 Nt 张图片，完成后使用 “用 Matlab 制作 gif 动画” 中给出的代码可以制作成 gif 动画。

代码 1：kepler.m

% 已知初始位置、发射速度、发射方向， 求轨道以及运动方程
function kepler
% === 参数设置 ===
k = -1; % -GMm, (m=1)
R = 1; alpha = 0; % 初始位置（极坐标）
v0 = 1.1; % 初速度（支持行矢量或标量）
beta = linspace(0, 4*pi/11, 6); % 发射仰角（支持行矢量或标量）
t_max = 8.94827; Nt = 100; % 模拟时间和步数
axis_param = 'auto'; % 坐标范围 [xmin,xmax,ymin,ymax] 或 'auto'
% ==============

% 匹配 v0 和 beta 的尺寸
N = max([numel(R), numel(alpha), numel(v0), numel(beta)]);
if isscalar(R), R = R*ones(1, N); end
if isscalar(alpha), alpha = alpha*ones(1, N); end
if isscalar(v0), v0 = v0*ones(1, N); end
if isscalar(beta), beta = beta*ones(1, N); end

if any(abs(cos(beta))<eps)
    error('不支持纯径向运动！');
end

% 画参考圆
close all; figure;
if all(R(:) == R(1))
    tmp = linspace(0, 2*pi, 500);
    plot(R(1)*cos(tmp), R(1)*sin(tmp), 'k--'); hold on;
end
scatter(0, 0); hold on; axis equal;
xlabel x; ylabel y;
N = numel(v0); t_cell = cell(1,N); th_cell = t_cell;
p_list = zeros(1,N); e_list = p_list; th0_list = p_list;
for i = 1:N
    [p_list(i), e_list(i), th0_list(i), t_cell{i}, th_cell{i}] = ...
        kepler1(k, R(i), alpha(i), v0(i), beta(i), t_max, Nt); hold on;
    axis(axis_param);
end

% 生成动画
h = nan(1,N); % point handles
for it = 1:10000000
    exit = true;
    for i = 1:N
        s = sign(k);
        if it > size(th_cell{i})
            h(i) = nan; continue;
        end
        exit = false;
        th = th_cell{i}(it);
        r = -s*p_list(i) / (1 + s*e_list(i)*cos(th - th0_list(i)));
        if (~isnan(h(i))), delete(h(i)); end
        h(i) = scatter(r*cos(th), r*sin(th), 'r^');
    end
    if exit
        break;
    end
    title(['t = ' num2str(t_cell{1}(it) - t_cell{1}(1), '%.2f')]);
    saveas(gcf, [num2str(it) '.png']);
end
end

% 计算并画出一条轨道
% 轨道方程为 r = -sign(k) * p / (1 + sign(k)*e*cos(th - th0));
function [p, e, th0, t, th_t] = kepler1(k, R, alpha, v0, beta, t_max, Nt)
E = v0^2/2 + k/R; % 轨道总机械能
L = v0*cos(beta)*R; % 角动量
e = sqrt(1 + 2*E*L^2/k^2); % 离心率
p = L^2/abs(k); % 半通径

% 画图
dth = pi/1000; Nth = 500;
if e < 1 % 椭圆
    th = linspace(0, 2*pi, Nth);
elseif e == 1 % 抛物线
    th = linspace(dth, 2*pi-dth, Nth);
else % 双曲线
    gamma = acos(1/e); % 渐进张角 atan(b/a)
    if k < 0
        th = linspace(gamma+dth, 2*pi-gamma-dth, Nth);
    else
        th = linspace(pi-gamma+dth, pi+gamma-dth, Nth);
    end
end

% 计算圆锥曲线需要旋转的角度 th0
s = sign(k);
if e == 0 % 圆
    th0 = 0;
else
    if sin(2*beta) >= 0
        th0 = alpha + acos((-s - p/R)/e);
    else
        th0 = alpha - acos((-s - p/R)/e);
    end
end
th0 = real(mod(real(th0)+pi, 2*pi) - pi);
scatter(R*cos(alpha), R*sin(alpha));
th = th + th0;

% 画出轨道
r = -s * p./(1 + s*e*cos(th-th0));
x = r .* cos(th); y = r .* sin(th);
plot(x, y); axis equal;

% 含时运动
% 时间格点
rot_dir = sign(cos(beta)); % 1: 逆时针 -1: 顺时针
t0 = rot_dir*t_Dth(p, e, k, alpha-th0-pi); % 初始时间
t = linspace(t0, t0+t_max, Nt);
th_t = rot_dir*Dth_t(p, e, k, dth, t) + th0 + pi;
end

% 求 t(\Delta\theta)
% Dth 支持矢量
% 对椭圆 Dth 支持任意实数，会自动 wrap 到 (-pi, pi]
function t = t_Dth(p, e, k, Dth)
if Dth ~= pi
    Dth = mod(Dth+pi, 2*pi) - pi;
end
if e == 1 && k < 0 % 抛物线
    L = sqrt(p*abs(k));
    tmp = tan(Dth/2);
    t = L^3/(2*k^2) * (tmp + tmp.^3/3);
elseif e < 1 && k < 0 % 椭圆
    a = p/(1-e^2);
    r = p./(1 + e*cos(Dth));
    psi = sign(Dth).*real(acos((1 - r/a)/e));
    t = sqrt(a^3/(abs(k)))*(psi - e*sin(psi));
elseif e > 1 && k < 0 % 双曲线 (k<0)
    a = p/(e^2-1);
    r = p./(1 + e*cos(Dth));
    xi = sign(Dth).*real(acosh((1 + r/a)/e));
    t = sqrt(a^3/(abs(k)))*(e*sinh(xi) - xi);
elseif e > 1 && k > 0 % 双曲线 (k>0)
    % 此时 Dth = th
    a = p/(e^2-1);
    r = p./(1 - e*cos(Dth));
    xi = sign(Dth).*acosh((-r/a-1)/e);
    t = sqrt(a^3/(abs(k)))*(e*sinh(xi) + xi);
else
    error('参数错误或未实现!');
end
end

% 已知时间求解 \Delta\theta
% t 支持矢量
% 对于椭圆， 支持 t 为任意实数
function Dth = Dth_t(p, e, k, dth, t)
N = numel(t); Dth = zeros(1,N);
if e < 1 % 椭圆
    a = p/(1-e^2);
    T = 2*pi*a^1.5*sqrt(1/abs(k));
    t = mod(t + T/2, T) - T/2;
    Dth_range = [-pi, pi];
elseif e == 1 % 抛物线
    Dth_range = [-pi+dth, pi-dth];
else % 双曲线
    th1 = atan(sqrt(e^2-1)); % 渐进张角的一半
    if k < 0
        Dth_range = [-pi+th1+dth, pi-th1-dth];
    else
        Dth_range = [-th1+dth, th1-dth];
    end
end
% 二分法解方程
for i = 1:N
    fun = @(Dth) t_Dth(p, e, k, Dth) - t(i);
    Dth(i) = fzero(fun, Dth_range);
end
end

　　运行结果见图 1 。这里给出另一组参数，计算以 $45^{\circ}$ 仰角和不同初速度发射的情况，替换到以上代码中即可，结果见图 2 。

% === 参数设置 ===
k = -1; % -GMm, (m=1)
R = 1; alpha = 0; % 初始位置（极坐标）
v0 = linspace(1, 1.5, 6); % 初速度（支持行矢量或标量）
beta = pi/4; % 发射仰角（支持行矢量或标量）
t_max = 4*pi; Nt = 100; % 模拟时间和步数
axis_param = [-1.1,3.5,-1.1,4]; % 坐标范围 [xmin,xmax,ymin,ymax] 或 'auto'
% ==============

　　图 3 中卢瑟福散射的参数：

% === 参数设置 ===
k = 1; % -GMm, (m=1)
N = 8; % 轨道数
x0 = 4*ones(1,N); % 初始直角坐标 x
y0 = linspace(0.1, 1.5, N); % 初始直角坐标 y
R = sqrt(x0.^2 + y0.^2); % 初始距离
alpha = atan2(y0, x0); % 初始极角
v0 = 1.6; % 初速度（支持行矢量或标量）
beta = -pi/2 + alpha; % 发射仰角（支持行矢量或标量）
t_max = 7; Nt = 100; % 模拟时间和步数
axis_param = [-4,4,-0.5,5]; % 坐标范围 [xmin,xmax,ymin,ymax] 或 'auto'
% ==============

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