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

                     

贡献者： addis

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

1. 轨道的计算

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

其中 $k$ 为非零实常数。由于质点的运动只和初始位置、初始速度和 $k$ 有关，与质点的质量无关，所以我们令 $m=1$ 并不影响一般性。令发射位置的极坐标为 $(R,\alpha)$，初速度为 $v_0$，发射仰角为 $\beta$（顺时针增量为正，$\beta=0$ 时质点延逆时针方向水平发射）。

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

其中当 $k>0$ 时 $s=1$，$k<0$ 时 $s=-1$。要求 $\theta_0$，可以代入初始条件，即 $r(\alpha) = R$。得 这里有两个解是因为圆锥曲线和半径为 $R$ 的圆必定有两个交点。我们需要根据发射角度来判断哪个是发射点，稍作分析易得当 $\sin2\beta \geqslant 0$ 时上式取正号，否则取负号。

2. 运动的计算

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

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

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

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

　　 求出发射质点的初始时刻为 $t_0 = t(\alpha - \theta_0 - \pi)$，然后生成等间距时间格点 $t_i$（$i=0,1,2,\dots$），用二分法解方程得到对应的 $\theta_i$ 即可。

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

% 已知初始位置、发射速度、发射方向， 求轨道以及运动方程
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'
% ==============