最近帮朋友找三体系统中的周期轨道，作为基础，先写程序验证了shooting法对单体问题的可行性。

这里的单体问题也就是kepler问题，即单个粒子在库仑势里的运动。

在kepler框架中，有个lambert问题，即给定起点和终点，以及时间间隔，找一条连接起点和终点的真实可行的轨道。本质上就是要确定粒子在起点的初始速度，以便其在给定时间间隔之后到达指定目的地。

数值上，有两个途径值得考虑。第一个办法是最小化作用量，因为这个问题是个标准的Euler-Lagrange问题。这个办法比较简单，不过方便使用的差分格式太粗糙，所以精度不高。第二个办法是shooting，先猜测一个初始速度，从初始条件出发积分到固定时间间隔，最终到达的地方肯定跟目标有差距，于是我们需要对初始速度做调整，调整的办法可以是基于newton法。

下面是按照第二个办法写的程序（程序尚有很大改进空间）。

% 2d kepler problem

clear all; close all; clc;

T = 20;

xstart = [2;0];

xend = [0 ;2];

iter = 10;

v00 = [1;0.7];

vincrement = 0.001;

relax = 1;

% euler-lagrange equation

N2 = 1e5;

dt = T/N2;

for ss = 1 : iter

    x0 = xstart;             % initial positioN2

    v0 = v00;

    y = [x0; v0];

    ylist = zeros(4, N2+1);

    Elist = zeros(1, N2+1);

    ylist(:,1) = y;

    Elist(1,1) = 0.5*norm(v0)^2 - 1/norm(x0);

    for s = 2: (N2+1)

        ytemp = y;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K1 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        ytemp = y + 0.5*K1;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K2 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        ytemp = y + 0.5*K2;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K3 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        ytemp = y + K3;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K4 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        y = y + (K1+ 2*K2 + 2*K3 + K4)/6;

        ylist(:,s) = y;

        Elist(1,s) = 0.5*norm(y(3:4))^2 - 1/norm(y(1:2));

    end

    diff = ylist(1:2, end) - xend;

    %         h3 = figure;

    plot(ylist(1,:), ylist(2,:), xstart(1), xstart(2),'*', xend(1), xend(2),'o')

    axis equal

    pause(0.1)

    Mov(ss) = getframe;

    x0 = xstart;             % initial position

    v0 = v00 + [vincrement;0];

    y = [x0; v0];

    ylist = zeros(4, N2+1);

    Elist = zeros(1, N2+1);

    ylist(:,1) = y;

    Elist(1,1) = 0.5*norm(v0)^2 - 1/norm(x0);

    for s = 2: (N2+1)

       ytemp = y;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K1 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        ytemp = y + 0.5*K1;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K2 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        ytemp = y + 0.5*K2;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K3 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        ytemp = y + K3;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K4 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        y = y + (K1+ 2*K2 + 2*K3 + K4)/6;

        ylist(:,s) = y;

        Elist(1,s) = 0.5*norm(y(3:4))^2 - 1/norm(y(1:2));

    end

    diff1 = ylist(1:2, end) - xend;

    x0 = xstart;             % initial position

    v0 =v00+ [0; vincrement];

    y = [x0; v0];

    ylist = zeros(4, N2+1);

    Elist = zeros(1, N2+1);

    ylist(:,1) = y;

    Elist(1,1) = 0.5*norm(v0)^2 - 1/norm(x0);

    for s = 2: (N2+1)

       ytemp = y;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K1 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        ytemp = y + 0.5*K1;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K2 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        ytemp = y + 0.5*K2;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K3 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        ytemp = y + K3;

        r = sqrt(ytemp(1)^2 + ytemp(2)^2);

        K4 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

        y = y + (K1+ 2*K2 + 2*K3 + K4)/6;

        ylist(:,s) = y;

        Elist(1,s) = 0.5*norm(y(3:4))^2 - 1/norm(y(1:2));

    end

    diff2 = ylist(1:2, end) - xend;

    %  v01 =  v00 -[vincrement , 0 ; 0, vincrement] * ([diff1 - diff, diff2- diff]\diff)

    v00 = v00 - relax * ([diff1 - diff, diff2 - diff]\([vincrement , 0 ; 0, vincrement]* diff))

end

% movie(Mov)

movie2gif(Mov, 'lambert.gif','DelayTime',0.3)

下面的动图展示了一个收敛的例子（*为起始点，o为目标终点）：

这样的lambert问题在宇航和军事上的应用是显然的。比如，也许15年后，中美之间就会爆发第三次世界大战，届时中国肯定会试图击落美方的间谍卫星。间谍卫星不携带燃料，不能变轨，所以我们很容易预测其在未来某个时刻的位置，我军便需要发射一颗导弹在该时刻到达该位置。确定导弹发射位置和时间后，我们就需要解一个lambert问题（如果是弹道导弹的话）。