1地月系统模拟

万有引力的公式为：
F=−GMmR2F=-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

其中近地点月球的速度通过下面公式求得：
Vnear2=GM2R(R+r)rV_{near}^{2}=GM\frac{2R}{\left( R+r \right) r} Vnear2=GM(R+r)r2R
在已知初始条件和参数后，我们开始建立模型，对月球的运行轨道建立方程，可以得到在任意时刻月球的加速度a为：
ax=−GMx2+y2xx2+y2a_x=-\frac{GM}{x^2+y^2}\frac{x}{\sqrt{x^2+y^2}} ax=−x2+y2GMx2+y2x
ay=−GMx2+y2yx2+y2a_y=-\frac{GM}{x^2+y^2}\frac{y}{\sqrt{x^2+y^2}} ay=−x2+y2GMx2+y2y
其中ax为x轴方向上的加速度，ay为y轴方向上的加速度。该方程属于典型的常微分方程ODE，整理为一阶微分方程的形式：
[xyuv]′=[uvaxay]=[uv−GMx2+y2xx2+y2−GMx2+y2yx2+y2]\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+y2GMx2+y2x−x2+y2GMx2+y2y⎦⎥⎥⎥⎤
其中上式的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)-1yn=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
endfunction 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双星问题模型

上一章由于地球质量远远大于月球质量，所以假设地球不动，只建立月球的模型。
但是当两个星球质量比较接近时，就不能忽略两个星球间的相互作用了。

所以对于双星问题，原先的常微分方程需要更改为：
[m1x1y1u1v1m2x2y2u2v2]′=[0u1v1ax1ay10u2v2ax2ay2]=[0u1v1−Gm2(x1−x2)((x1−x2)2+(x1−x2)2)1.5−Gm2(y1−y2)((x1−x2)2+(x1−x2)2)1.50u2v2−Gm1(x2−x1)((x2−x1)2+(x2−x1)2)1.5−Gm1(y2−y1)((x2−x1)2+(x2−x1)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] ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡m1x1y1u1v1m2x2y2u2v2⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤′=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡0u1v1ax1ay10u2v2ax2ay2⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡0u1v1−((x1−x2)2+(x1−x2)2)1.5Gm2(x1−x2)−((x1−x2)2+(x1−x2)2)1.5Gm2(y1−y2)0u2v2−((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)-1yn=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
endfunction 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:Nax(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

双星的轨迹如下（参数时随机的，所以不一样很正常）：

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

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 trueT=T+1;t=T*h;%计算下一步if T<=Nny=ODE_RK4_hyh2(t,h,y1(:,T-1));y1(:,T)=y;elsey2=y1;y=ODE_RK4_hyh2(t,h,y1(:,end));y1=[y2(:,2:end),y];end%绘图if mod(T,100)==1ax=figure(1);ax.Color=[0 0 0.1];clfx_medium=0;y_medium=0;if T<=Nn%初始还没有超过线条点数hold onfor 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);endhold offx_medium=x_medium/Mn;y_medium=y_medium/Mn;else%超过线条点数后hold onfor 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);endhold offx_medium=x_medium/Mn;y_medium=y_medium/Mn;endxlim([x_medium-1,x_medium+1]);ylim([y_medium-1,y_medium+1]);%W_deta=Gravity_Energy(y)/Waxis offpause(0.1)end
endfunction 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);
endfunction 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:Nax(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星运动（恒星加行星）

代码本身还是有很多bug的，比如
1遇到两个行星特别近的时候，就会出现差分步长太大导致的数值误差，表现为两个行星突然朝相反的方向快速飞出。如果设计为变步长的差分算法，或者设定一个碰撞机制，可以避免这个问题。
2视角问题，遇到多个行星的移动、拉近、远离等视角自动跟踪的方法还不是太满意，有时候表现为晃来晃去的。如果能够做到交互式视角的话，可以避免这个问题，但是交互界面需要GUI或者appdesigner，比较复杂。

3维空间

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

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

利用Matlab将图片转换成素描(简笔画)风格
题目: 利用Matlab将图片转换成素描(简笔画)风格记得曾经看过别人的网络头像,是那种类似简笔画或素描的图片,一直以来都想做一个类似的头像,但始终不得要领.今天当我看到文献[1]中的图5.28时( ...
matlab图像处理将两个目标合成一个,利用MATLAB实现医学图像处理与分析
[实例简介] 利用MATLAB实现医学图像处理与分析边缘是图像最基本的特征.所谓边缘是指图像周围像素灰度有阶跃变化或屋顶状变化的像素的集合, 它存在于目标与背景.目标与目标.区域与区域.基元与基元之 ...
matlab系统的根轨迹,实验五利用MATLAB绘制系统根轨迹
<实验五利用MATLAB绘制系统根轨迹>由会员分享,可在线阅读,更多相关<实验五利用MATLAB绘制系统根轨迹(6页珍藏版)>请在人人文库网上搜索. 1.实验五利用MAT ...
matlab对图像进行增强,利用matlab对图像进行增强处理.doc
利用matlab对图像进行增强处理.doc 郑州轻工业学院课程设计任务书题目利用MATLAB对图像进行增强处理专业.班级电子信息工程07级学号姓名主要内容.基本要求.主要参考资料等:主要内容:在图像形 ...
matlab如何截取图像的中间部分_利用matlab提取并分割RGB图像中的某一个已知像素值的图像...
已知一副RGB图像中的的像素值,利用matlab将其分割出来并以二进制图像形式显示: %extract.m clear all; I=imread('new_original.png'); figur ...
【Matlab 控制】利用 Matlab Function 绘制分段函数
利用 Switch block 利用 Matlab Function block function [mean, stdev] = fcn(vals) % #codegen% calculates a ...
幅度调制信号 matlab,《利用MATLAB实现信号的幅度调制与解调.doc
<利用MATLAB实现信号的幅度调制与解调课程设计论文姓名:姜勇学院:机电与车辆工程学院专业:电子信息工程2班学号:1665090208 安徽科技学院学年第学期 < > ...
利用Matlab优化工具箱求解旅行商最短路径问题
前面介绍了利用Matlab二元整数规划求解数独问题,对于另一个问题-旅行商问题也可以用它来求解. 旅行商问题就是找到经过所有站点的最短闭合路径,如下图为在美国地图框架内产生的200个旅行站点,而旅行商 ...
利用matlab处理点云
本文主要分享利用matlab点云工具的相关模块来处理点云,并通过点云轮廓对点云体积进行简单的估计测量. 目录利用matlab处理点云目录主要的操作流程图 2具体流程 1 点云的读入和显示 2 点 ...

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

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