天天看点

利用matlab实现三体问题(双星、3星、多星运动)1地月系统模拟2双星问题模型3多星问题模拟

利用matlab实现3体问题(双星、3星、多星运动)

  • 1地月系统模拟
  • 2双星问题模型
  • 3多星问题模拟
    • 2维平面
    • 3维空间

这篇文章随便水一下,写一个简单的3体问题,权当娱乐

利用matlab实现三体问题(双星、3星、多星运动)1地月系统模拟2双星问题模型3多星问题模拟

1地月系统模拟

万有引力的公式为:

F = − G M m R 2 F=-G\frac{Mm}{R^2} F=−GR2Mm​

查询地月系统的参数,可以得到:

参数 数值
引力常数 G 6.67e-11
地球质量 M 5.965e+24 Kg
近地点距离 r 363300000 m
远地点距离 R 405493000 m
近地点月球速度 Vnear 1074.82 m/s
轨道离心率 0.0549

其中近地点月球的速度通过下面公式求得:

V n e a r 2 = G M 2 R ( R + r ) r V_{near}^{2}=GM\frac{2R}{\left( R+r \right) r} Vnear2​=GM(R+r)r2R​

在已知初始条件和参数后,我们开始建立模型,对月球的运行轨道建立方程,可以得到在任意时刻月球的加速度a为:

a x = − G M x 2 + y 2 x x 2 + y 2 a_x=-\frac{GM}{x^2+y^2}\frac{x}{\sqrt{x^2+y^2}} ax​=−x2+y2GM​x2+y2

​x​

a y = − G M x 2 + y 2 y x 2 + y 2 a_y=-\frac{GM}{x^2+y^2}\frac{y}{\sqrt{x^2+y^2}} ay​=−x2+y2GM​x2+y2

​y​

其中ax为x轴方向上的加速度,ay为y轴方向上的加速度。该方程属于典型的常微分方程ODE,整理为一阶微分方程的形式:

[ x y u v ] ′ = [ u v a x a y ] = [ u v − G M x 2 + y 2 x x 2 + y 2 − G M x 2 + y 2 y x 2 + y 2 ] \left[ \begin{array}{c} x\\ y\\ u\\ v\\ \end{array} \right] ^{'}=\left[ \begin{array}{c} u\\ v\\ ax\\ ay\\ \end{array} \right] =\left[ \begin{array}{c} u\\ v\\ -\frac{GM}{x^2+y^2}\frac{x}{\sqrt{x^2+y^2}}\\ -\frac{GM}{x^2+y^2}\frac{y}{\sqrt{x^2+y^2}}\\ \end{array} \right] ⎣⎢⎢⎡​xyuv​⎦⎥⎥⎤​′=⎣⎢⎢⎡​uvaxay​⎦⎥⎥⎤​=⎣⎢⎢⎢⎡​uv−x2+y2GM​x2+y2

​x​−x2+y2GM​x2+y2

​y​​⎦⎥⎥⎥⎤​

其中上式的uv为x方向上的速度和y方向上的速度。为了提高精度,将方程带入到4阶RK求解器中,即可求解得到月球的轨迹。

matlab代码如下:

clear
clc
close all
%初始化
x0=363300*10^3;%初始近地位置
y0=0;
u0=0;
v0=1074.82;%初始近地速度
%设置计算时间
h=0.5*60*60;%步长为0.5个小时
t=0:h:27.5*24*60*60;%计算总时间维27.5天

y=ODE_RK4_hyh(t,h,[x0;y0;u0;v0]);

plot(y(1,:),y(2,:))%绘图

a=(max(y(1,:))-min(y(1,:)))/2;%椭圆半长轴
b=(max(y(2,:))-min(y(2,:)))/2;%椭圆半短轴
e=sqrt(a^2-b^2)/a;%椭圆轨道离心率

function y=ODE_RK4_hyh(x,h,y0)
%4阶RK方法
%h间隔为常数的算法
y=zeros(size(y0,1),size(x,2));
y(:,1)=y0;
for ii=1:length(x)-1
    yn=y(:,ii);
    xn=x(ii);
    K1=Fdydx(xn,yn);
    K2=Fdydx(xn+h/2,yn+h/2*K1);
    K3=Fdydx(xn+h/2,yn+h/2*K2);
    K4=Fdydx(xn+h,yn+h*K3);
    y(:,ii+1)=yn+h/6*(K1+2*K2+2*K3+K4);
end
end

function dydx=Fdydx(x,y)
%将原方程整理为dy/dx=F(y,x)的形式
G=6.67E-11;%引力常数
M=5.965E24;%地球质量
ax=-G*M*y(1)/(y(1).^2+y(2).^2)^1.5;
ay=-G*M*y(2)/(y(1).^2+y(2).^2)^1.5/1;
dydx=[y(3);y(4);ax;ay];
end
           

验证得到的周期大概小于27.5天(这个我没细算)。验证得到的偏心率为0.0549,和百度百科上完全一致,证明该方法完全可以用在定量计算上。下图为平平无奇的月球轨道,我没有用axis equal调整x轴和y轴比例,所以看上去离心率有点大。

利用matlab实现三体问题(双星、3星、多星运动)1地月系统模拟2双星问题模型3多星问题模拟

我之后为了效果,也为了计算省事,不再设置真实的G、M等参数。

2双星问题模型

上一章由于地球质量远远大于月球质量,所以假设地球不动,只建立月球的模型。

但是当两个星球质量比较接近时,就不能忽略两个星球间的相互作用了。

所以对于双星问题,原先的常微分方程需要更改为:

[ m 1 x 1 y 1 u 1 v 1 m 2 x 2 y 2 u 2 v 2 ] ′ = [ 0 u 1 v 1 a x 1 a y 1 0 u 2 v 2 a x 2 a y 2 ] = [ 0 u 1 v 1 − G m 2 ( x 1 − x 2 ) ( ( x 1 − x 2 ) 2 + ( x 1 − x 2 ) 2 ) 1.5 − G m 2 ( y 1 − y 2 ) ( ( x 1 − x 2 ) 2 + ( x 1 − x 2 ) 2 ) 1.5 0 u 2 v 2 − G m 1 ( x 2 − x 1 ) ( ( x 2 − x 1 ) 2 + ( x 2 − x 1 ) 2 ) 1.5 − G m 1 ( y 2 − y 1 ) ( ( x 2 − x 1 ) 2 + ( x 2 − x 1 ) 2 ) 1.5 ] \left[ \begin{array}{c} m_{1}\\ x_{1}\\ y_{1}\\ u_{1}\\ v_{1}\\ m_{2}\\ x_{2}\\ y_{2}\\ u_{2}\\ v_{2}\\ \end{array} \right] ^{'}=\left[ \begin{array}{c} 0\\ u_{1}\\ v_{1}\\ ax_{1}\\ ay_{1}\\ 0\\ u_{2}\\ v_{2}\\ ax_{2}\\ ay_{2}\\ \end{array} \right] =\left[ \begin{array}{c} 0\\ u_{1}\\ v_{1}\\ -\frac{Gm_{2}(x_{1}-x_{2})}{((x_{1}-x_{2})^2+(x_{1}-x_{2})^2)^{1.5}}\\ -\frac{Gm_{2}(y_{1}-y_{2})}{((x_{1}-x_{2})^2+(x_{1}-x_{2})^2)^{1.5}}\\ 0\\ u_{2}\\ v_{2}\\ -\frac{Gm_{1}(x_{2}-x_{1})}{((x_{2}-x_{1})^2+(x_{2}-x_{1})^2)^{1.5}}\\ -\frac{Gm_{1}(y_{2}-y_{1})}{((x_{2}-x_{1})^2+(x_{2}-x_{1})^2)^{1.5}}\\ \end{array} \right] ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​m1​x1​y1​u1​v1​m2​x2​y2​u2​v2​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​′=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​0u1​v1​ax1​ay1​0u2​v2​ax2​ay2​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​0u1​v1​−((x1​−x2​)2+(x1​−x2​)2)1.5Gm2​(x1​−x2​)​−((x1​−x2​)2+(x1​−x2​)2)1.5Gm2​(y1​−y2​)​0u2​v2​−((x2​−x1​)2+(x2​−x1​)2)1.5Gm1​(x2​−x1​)​−((x2​−x1​)2+(x2​−x1​)2)1.5Gm1​(y2​−y1​)​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​

每一个行星有5个信息,分别为质量、位置xy、速度uv。其中质量不发生改变,位置的导数为速度,速度的导数为加速度。

在matlab编程中,由于在计算引力过程中,两个行星还可以穷举,但是多个行星之后,就要考虑循环,第一层循环遍历每一个行星,第二层循环计算每一个行星与其它行星之间的引力之和。这种双循环在matlab中可能会导致计算速度太慢,所以我利用matlab里的meshgrid方法,避免了使用循环,这也相当于向量加速。

代码如下:

clear
clc
close all

%计算双星系统,定性计算,没有代入实际物理参数
mxyuv1=rand(5,1);%行星1,所有参数都随机
mxyuv2=rand(5,1);%行星2,所有参数都随机
h=1e-5;%步长
t=0:h:1;%时间

y=ODE_RK4_hyh(t,h,[mxyuv1;mxyuv2]);

plot(y(2,:),y(3,:),y(7,:),y(8,:))%绘轨迹图

function y=ODE_RK4_hyh(x,h,y0)
%4阶RK方法
%h间隔为常数的算法
y=zeros(size(y0,1),size(x,2));
y(:,1)=y0;
for ii=1:length(x)-1
    yn=y(:,ii);
    xn=x(ii);
    K1=Fdydx(xn,yn);
    K2=Fdydx(xn+h/2,yn+h/2*K1);
    K3=Fdydx(xn+h/2,yn+h/2*K2);
    K4=Fdydx(xn+h,yn+h*K3);
    y(:,ii+1)=yn+h/6*(K1+2*K2+2*K3+K4);
end
end

function dydx=Fdydx(x,y)
%将原方程整理为dy/dx=F(y,x)的形式
N=numel(y)/5;%N个球体
G=1;%引力系数,这里只做展示,不做定量计算
%计算两个星球之间的引力
m=y(1:5:end);
x0=y(2:5:end);
y0=y(3:5:end);
[x1,x2]=meshgrid(x0,x0);
[y1,y2]=meshgrid(y0,y0);
[m1,m2]=meshgrid(m,m);
dx=x1-x2;
dy=y1-y2;
ax=-G.*dx./(dx.^2+dy.^2).^(1.5).*m2;
ay=-G.*dy./(dx.^2+dy.^2).^(1.5).*m2;
for k=1:N
    ax(k,k)=0;
    ay(k,k)=0;
end

%建立dydx
dydx=zeros(size(y));
%1 质量不变
dydx(1:5:N*5-4)=0;
%2 x导数
dydx(2:5:N*5-3)=y(4:5:N*5-1);
%3 y导数
dydx(3:5:N*5-2)=y(5:5:N*5-0);
%4 x加速度
dydx(4:5:N*5-1)=sum(ax,1);
%3 y加速度
dydx(5:5:N*5-0)=sum(ay,1);
end
           

双星的轨迹如下(参数时随机的,所以不一样很正常):

利用matlab实现三体问题(双星、3星、多星运动)1地月系统模拟2双星问题模型3多星问题模拟

可以看到,两个行星在相互吸引后在不到一个周期内又快速互相甩出,当然之后可能还会再次吸引再次重复这个周期,但是我这里时间设置比较少。

3多星问题模拟

2维平面

代码还是以上一章的代码为主,星球的增加参考上一章节的即可,然后对最终演示做了亿点优化。

clear
clc
close all

%计算双星系统,不停循环,只显示一部分点
%初始条件
Mn=3;%更改星星数量

h=5e-5;%步长
Nn=5000;%线段长度
%t=0:h:1;

% y0=rand(Mn*5,1);%随机

% y0=[100,0   ,0  ,0  ,0  ,...
%     1  ,1.1,0.05  ,2.1,9  ,...
%     1  ,1.1,-0.05  ,-2.1  ,9  ]';%3星,双星加恒星

% y0=[100,0   ,0  ,0  ,0  ,...
%     0.1,0.2 ,0  ,0  ,22 ,...
%     0.5,0.35,0  ,0  ,17 ,...
%     1  ,0.6 ,0  ,0  ,13 ,...
%     0.4  ,0.9 ,0  ,0  ,10]';%5星,恒星加4个行星

y0=[11,0    ,0    ,4.5  ,4  ,...
    11,0.5 ,-0.132,-2  ,-1.8 ,...
    11,-0.5,0.132 ,-2  ,-1.8]';%3星,8字形运动,没调试好


y1=zeros(size(y0,1),Nn);


T=1;t=T*h;
y1(:,T)=y0;

while true
    T=T+1;
    t=T*h;
    %计算下一步
    if T<=Nn
        y=ODE_RK4_hyh2(t,h,y1(:,T-1));
        y1(:,T)=y;
    else
        y2=y1;
        y=ODE_RK4_hyh2(t,h,y1(:,end));
        y1=[y2(:,2:end),y];
    end
    

    
    %绘图
    if mod(T,100)==1
        ax=figure(1);
        ax.Color=[0 0 0.1];
        clf
        x_medium=0;y_medium=0;
        if T<=Nn
            %初始还没有超过线条点数
            hold on
            for k=1:Mn
                %计算颜色
                color_hsv=colormap(hsv);
                color_N=length(color_hsv);
                F_color=color_hsv(round(color_N/Mn*k),:);
                F_color=F_color*0.6+[1,1,1]*0.4;
                cdata=[linspace(0,F_color(1),T+1)',linspace(0,F_color(2),T+1)',linspace(0.1,F_color(3),T+1)'];
                cdata=reshape(cdata,T+1,1,3);
                %绘制线条
                patch([y1(k*5-5+2,1:T),NaN],[y1(k*5-5+3,1:T),NaN],1:T+1,'EdgeColor','interp','Marker','none',...
                    'MarkerFaceColor','flat','CData',cdata,'LineWidth',1);
                x_medium=x_medium+y1(k*5-5+2,T);
                y_medium=y_medium+y1(k*5-5+3,T);
            end
            hold off
            x_medium=x_medium/Mn;
            y_medium=y_medium/Mn;
        else
            %超过线条点数后
            hold on
            for k=1:Mn
                %计算颜色
                color_hsv=colormap(hsv);
                color_N=length(color_hsv);
                F_color=color_hsv(round(color_N/Mn*k),:);
                F_color=F_color*0.6+[1,1,1]*0.4;
                cdata=[linspace(0,F_color(1),Nn+1)',linspace(0,F_color(2),Nn+1)',linspace(0.1,F_color(3),Nn+1)'];
                cdata=reshape(cdata,Nn+1,1,3);
                %绘制线条
                patch([y1(k*5-5+2,:),NaN],[y1(k*5-5+3,:),NaN],1:Nn+1,'EdgeColor','interp','Marker','none',...
                    'MarkerFaceColor','flat','CData',cdata,'LineWidth',1);
                x_medium=x_medium+y1(k*5-5+2,end);
                y_medium=y_medium+y1(k*5-5+3,end);
            end
            hold off
            x_medium=x_medium/Mn;
            y_medium=y_medium/Mn;
        end
        xlim([x_medium-1,x_medium+1]);
        ylim([y_medium-1,y_medium+1]);
        
        
        %W_deta=Gravity_Energy(y)/W
        axis off
        pause(0.1)
        
    end
end



function y1=ODE_RK4_hyh2(x0,h,y0)
%4阶RK方法
%只向后计算1步h
%y1=zeros(size(y0,1),1);
K1=Fdydx(x0,y0);
K2=Fdydx(x0+h/2,y0+h/2*K1);
K3=Fdydx(x0+h/2,y0+h/2*K2);
K4=Fdydx(x0+h,y0+h*K3);
y1=y0+h/6*(K1+2*K2+2*K3+K4);
end

function dydx=Fdydx(x,y)
%将原方程整理为dy/dx=F(y,x)的形式
N=numel(y)/5;%N个球体
G=1;%引力系数,这里只做展示,不做定量计算
%计算两个星球之间的引力
m=y(1:5:end);
x0=y(2:5:end);
y0=y(3:5:end);
[x1,x2]=meshgrid(x0,x0);
[y1,y2]=meshgrid(y0,y0);
[m1,m2]=meshgrid(m,m);
dx=x1-x2;
dy=y1-y2;
ax=-G.*dx./(dx.^2+dy.^2).^(1.5).*m2;
ay=-G.*dy./(dx.^2+dy.^2).^(1.5).*m2;
for k=1:N
    ax(k,k)=0;
    ay(k,k)=0;
end

%建立dydx
dydx=zeros(size(y));
%1 质量不变
dydx(1:5:N*5-4)=0;
%2 x导数
dydx(2:5:N*5-3)=y(4:5:N*5-1);
%3 y导数
dydx(3:5:N*5-2)=y(5:5:N*5-0);
%4 x加速度
dydx(4:5:N*5-1)=sum(ax,1);
%3 y加速度
dydx(5:5:N*5-0)=sum(ay,1);

end

           

在文章最开头,我已经演示了8字形双星的运动(感觉轨道有误差,结果越修正越歪)。下面演示一个5星运动(恒星加行星)

利用matlab实现三体问题(双星、3星、多星运动)1地月系统模拟2双星问题模型3多星问题模拟

代码本身还是有很多bug的,比如

1遇到两个行星特别近的时候,就会出现差分步长太大导致的数值误差,表现为两个行星突然朝相反的方向快速飞出。如果设计为变步长的差分算法,或者设定一个碰撞机制,可以避免这个问题。

2视角问题,遇到多个行星的移动、拉近、远离等视角自动跟踪的方法还不是太满意,有时候表现为晃来晃去的。如果能够做到交互式视角的话,可以避免这个问题,但是交互界面需要GUI或者appdesigner,比较复杂。

3维空间

在方程上,增加了z方向的位置,以及对应的速度w。每一个行星的信息为7个。计算量上有所增加,但是具体算法没有什么变化。但是3维的视角问题更为严重,所以我暂时还没想好怎么演示,如果之后写关于交互的话,可能会把这一块补上。

继续阅读