利用matlab实现三体问题（双星、3星、多星运动）

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

1地月系统模拟

万有引力的公式为：
$$F=-G\frac{Mm}{R^2}$$
查询地月系统的参数，可以得到：

参数	数值
引力常数 G	6.67e-11
地球质量 M	5.965e+24 Kg
近地点距离 r	363300000 m
远地点距离 R	405493000 m
近地点月球速度 Vnear	1074.82 m/s
轨道离心率	0.0549
其中近地点月球的速度通过下面公式求得：
$$
V_{near}^{2}=GM\frac{2R}{\left( R+r \right) r}
$$
在已知初始条件和参数后，我们开始建立模型，对月球的运行轨道建立方程，可以得到在任意时刻月球的加速度a为：
$$
a_x=-\frac{GM}{x2+y2}\frac{x}{\sqrt{x2+y2}}
$$
$$
a_y=-\frac{GM}{x2+y2}\frac{y}{\sqrt{x2+y2}}
$$
其中ax为x轴方向上的加速度，ay为y轴方向上的加速度。该方程属于典型的常微分方程ODE，整理为一阶微分方程的形式：
$$
\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+y2}\frac{x}{\sqrt{x^2+y2}}\
-\frac{GM}{x^2+y2}\frac{y}{\sqrt{x^2+y2}}\
\end{array} \right]
$$
其中上式的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轴比例，所以看上去离心率有点大。

我之后为了效果，也为了计算省事，不再设置真实的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