正面教材分享 http://blog.sciencenet.cn/u/wdlang 70%的以色列人是无神论者,不过他们都相信上帝给了他们那块土地。这个世界经不起思考

博文

如何击落美国间谍卫星?

已有 3226 次阅读 2019-5-29 12:05 |个人分类:计算方法|系统分类:科普集锦

最近帮朋友找三体系统中的周期轨道,作为基础,先写程序验证了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.gif


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



http://wap.sciencenet.cn/blog-100379-1181810.html

上一篇:计算2的任意高次幂
下一篇:Mott的处女作

5 谢力 王振亭 王安良 刘全慧 陶勇

该博文允许注册用户评论 请点击登录 评论 (7 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-7-25 02:02

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部