一. 理论知识说明
二. 面向对象编程简单介绍
三. 此次程序设计所要用的类及代码
四. 主程序
五. 不同初始参数的三体运动展示
六. 后记
一. 下面先简单介绍一下编程计算中要用到的一些理论知识
把星体看做质点
我们要描述一个质点的运动就需要算位移，速度，加速度。我们可以采用笛卡尔坐标形式的运动微分方程来计算，我书上是这样写的：
由于没有解析解，我们只能用数值迭代的方法，设置一个时间步长，在每一次迭代中更新运动参数。开始时需要人为给定星体的如下这些初始值：
R：星体的半径
M：星体的质量
X：星体在X方向上的坐标
Y：星体在Y方向上的坐标
Z：星体在Z方向上的坐标
U：星体在X方向上的速度
V：星体在Y方向上的速度
W：星体在Z方向上的速度
R每次迭代中不变化，其实当质点处理了就不需要半径，后面程序中是用来显示星体大小的；
M每次迭代中不变化；
XYZUVW6个参数需要在是需要给定的初始值，而且在每次迭代中都要被重新计算。
依据输入的初始值我们可以计算出星体之间的万有引力，然后在把两个星体之间的万有引力分解到三个方向上去，算出三个方向上的万有引力，再利用三个方向上的万有引力算出三个方向上的加速度：
Ax：X方向上的加速度
Ay：Y方向上的加速度
Az：Z方向上的加速度
比如说现在有两个星体：
Star1(M1,X1,Y1,Z1,U1,V1,W1)
Star2(M2,X2,Y2,Z2,U2,V2,W2)
那么这两个星体之间的万有引力：（公众号号里要是能用LaTeX就好了...）
其中：g是万有引力常数R12是两个星体间的距离这个表示Star1对Star2的引力，那2对1的引力就是：-f12了。这里的f12是由2指向1的矢量（因为1吸引着2嘛），我们需要把它分解到三个方向上去:
Star1的坐标：XYZ1=(X1,Y1,Z1)Star2的坐标：XYZ2=(X2,Y2,Z2)
那么f12在XYZ三个方向上的分量是：
（其实就是乘了一个方向余弦的值）
这样我们干脆直接把f12用一个向量来表示：
在三体星系中由于有三个星体所以一个星体会受到两个万有引力。三个星体所受到的力的矢量和要为0，就是要构成首尾相连的多边形。所以对于这三个星体所受到的力，我按照如下的形式来写：星体1受到的力为：f31 和 -f12（等于f21）星体2受到的力为：f12 和 -f23（等于f32）星体3受到的力为：f23 和 -f31（等于f13）这样最后这这三个星体各自所受到的力（这个力已经用向量表示了）加一块正好是0，perfect。我们再把每个星体受到的力分别除上各自的质量就得到了每个星体在三个方向上的加速度：[Ax,Ay,Az]
得到了加速度，接下来们就可以根据如下的公式计算星体的坐标和速度了:
一定是要先计算坐标啊，因为计算下一时刻的位移要的是上一时刻的位置和速度，如果你先算了速度再算位移的话，那么你位移方程中的速度就不是上一时刻的速度了。
好了，到这星体的运动参数就都能被算出来了，下面我屡一下迭代计算的顺序：
计算各个星体间的相对距离
计算万有引力
计算加速度
计算位移
计算速度
classdef ClassName < handleproperties property1 property2 %.....可以在这里直接赋值endmethods function obj = ClassName(可以输入参数，也可以不输入) %这第一个function是类的构造函数，可以用于对类属性的初始化 %构造函数只会在类被装载(创建)时执行一次 %使用obj.property1 = ....，对类的property1属性进行赋值 %... end function method1_ReturnValue = method1(obj,para1,para2,...) %这是这个类的方法 %第一个参数必须是obj，在调用时不用输入 %para1,para2,...是你实际输入的第一个,第二个参数,... %调用或修改类的属性：obj.property1 %调用类的其他方法：obj.method2 end %....还可以继续添加其他方法end 类里面可以不写方法或属性。properties后面不加东西时为默认属性，Default。properties后面还加：Constant，Dependent，Hidden，这三种。由于后面的程序中我只用到了Constant，所以这里就只说Constant。Constant：就是把这个属性声明成常量，以后不变，例如本程序中的万有引力常数，如下：properties（Constant） G=6.67*10^(-11); end三. 下面开始码类的代码我们按照迭代计算的顺序来写 类 ：1· 相对距离可以作为一个类它的属性就一个：R，表示两个星体之间的距离；它不需要方法，写上构造函数就行。classdef RelativeDistance properties R%两个星体之间的相对距离 end methods function obj = RelativeDistance(XYZ1 , XYZ2) %XYZ1第一个星体的坐标向量 %XYZ2第二个星体的坐标向量 obj.R = sqrt( sum( (XYZ1 - XYZ2).^2 ) ); end endend2· 万有引力可以作为一个类它有一个默认属性：f 表示万有引力；一个常值属性：G表示万有引力常数；它同样也不需要任何方法，写一个构造函数就行。classdef RelativeDistance properties R%两个星体之间的相对距离 end methods function obj = RelativeDistance(XYZ1 , XYZ2) %XYZ1第一个星体的坐标向量 %XYZ2第二个星体的坐标向量 obj.R = sqrt( sum( (XYZ1 - XYZ2).^2 ) ); end endend3· 把星体作为一个类经过上面的关于这次编程所需要的理论的介绍，可以知道与一个星体的运动有关的属性有如下这些：
此外我们还需要添加用于控制球体显示的属性：H_star，它在这里是一个surface对象；还有球体颜色的属性H_Star_Color。
它的方法除了构造函数（用于初始化星体运动参数）外，还需要显示星体的方法Display 和 更新计算星体的运动参数的MotionParameter方法。
在贴上星体类的代码前我需要先介绍一下我的参数输入界面（通过这个界面手动输入你想要的初值）：
(这个小界面的代码我就不放了)
点击创建星体1，把星体一的这些参数放到结构体Star1中，并保存到当前文件夹下的ParameterFolder文件夹中；创建星体2和星体3同理。点击完成创建关闭这个输入参数的界面。（上面那些都是默认的初始值）
注意我在内部设置的坐标系范围为正负5000。所以XYZ值最好不要超过这个范围。
那么接下来可以编写Star类了：
classdef Star < handle properties R = 10 ; %星体的半径都为10 M%星体的质量 X%星体的横坐标 Y%星体的纵坐标 Z%星体的竖坐标 XYZ%星体的坐标向量 U%星体的横坐标方向速度 V%星体的纵坐标方向速度 W%星体的竖坐标方向速度 Ax%星体的横坐标方向加速度 Ay%星体的纵坐标方向加速度 Az%星体的竖坐标方向加速度 H_Star = [];%球体对象 H_Star_Color%球体对象的颜色 end methods function obj = Star(i) % 第i个星体 load('ParameterFolder\Star1.mat') load('ParameterFolder\Star2.mat') load('ParameterFolder\Star3.mat') obj.M = str2double( eval(['Star' , num2str(i) , '.M' ]) ); obj.M = obj.M * 10^19; obj.X = str2double( eval(['Star' , num2str(i) , '.X' ]) ); obj.Y = str2double(eval(['Star' , num2str(i) , '.Y' ]) ); obj.Z = str2double( eval(['Star' , num2str(i) , '.Z' ]) ); obj.U = str2double( eval(['Star' , num2str(i) , '.U' ]) ); obj.U = obj.U * 10^3; obj.V = str2double(eval(['Star' , num2str(i) , '.V' ]) ); obj.V = obj.V * 10^3; obj.W = str2double (eval(['Star' , num2str(i) , '.W' ]) ); obj.W = obj.W * 10^3; obj.XYZ = [obj.X , obj.Y , obj.Z]; switch i case 1 obj.H_Star_Color = [1 0 0]; obj.R = 30; case 2 obj.H_Star_Color = [0 1 0]; obj.R = 20; case 3 obj.H_Star_Color = [0 0 1]; obj.R = 10; end end function Display(obj) if obj.M ~= 0 [x , y , z] = sphere(20); x = (obj.X) + 20 .* obj.R .* x ; y = (obj.Y) + 12 .* obj.R .* y ; z = (obj.Z) + 10 .* obj.R .* z ; line('XData' , obj.X , 'YData' , obj.Y , 'ZData' , obj.Z , 'Color', obj.H_Star_Color ,'Marker','.','MarkerSize',4 ); obj.H_Star = mesh(x , y , z , 'FaceColor' , obj.H_Star_Color , 'FaceAlpha' , 0.1 , 'EdgeColor' , obj.H_Star_Color, 'EdgeAlpha' , 0.9); end end function MotionParameter(obj , f12 , f23 , f31 , i , t) %obj除外，剩下的是外部的输入 %前三个参数分别是星体之间的万有引力 %i 表示第几个星体 %t 表示运动的时间步长 if obj.M ~= 0 switch i case 1 A=(f31 - f12 ) / obj.M; obj.Ax = A(1); obj.Ay = A(2); obj.Az = A(3); case 2 A=(f12 - f23 ) / obj.M; obj.Ax = A(1); obj.Ay = A(2); obj.Az = A(3); case 3 A=(f23 - f31 ) / obj.M; obj.Ax = A(1); obj.Ay = A(2); obj.Az = A(3); end %必须！！！先算坐标，后算速度(因为位移方程中需要上一时刻的速度和位置) obj.X = obj.X + obj.U * t + 1/2 * obj.Ax * t^2; obj.Y = obj.Y + obj.V * t + 1/2 * obj.Ay * t^2; obj.Z = obj.Z + obj.W * t + 1/2 * obj.Az * t^2; obj.U = obj.U + obj.Ax*t; obj.V = obj.V + obj.Ay*t; obj.W = obj.W + obj.Az*t; %这个第一开始丢了，好久才发现，粗心了 obj.XYZ = [obj.X , obj.Y , obj.Z]; end end endend(eval 真是好用)
类中的Display方法就是画一个球心为(X,Y,Z)的Sphere；MotionParameter方法就是按照上面介绍的理论按步骤进行计算，所以这里也没有难度。
四. 下面该讲主程序了
显示运动的界面：
在这里我放了两个坐标轴1和2，坐标轴1用于显示背景，坐标轴二用于显示星体的空间运动。因为背景图片是二维的，运动是三维的运动，所以你想要背景只能用两个坐标系，而且放背景的坐标系要放在下方。
点击3中的设置参数就弹出了刚刚上面讲的设置参数的那个界面。点击运行就开始显示运动了。另外三个控件的回调不是重点，这里就不介绍了。所以下面只介绍 ‘运行Toggle’ 的回调
1. 在一开始我们需要把类实例化出来，生成三个星体对象；
2. 如果星体的质量不是零即该星体存在，就把它们显示出来；
3. 紧接着就开始迭代计算了，在每一次迭代的开始部分都要检查一下是否发生了碰撞，然后按照上面理论介绍最后那块的5步一步一步来。
注意时间步长的设定：如果是那种混沌状态的运动，时间步长大了就看不出来，因为在这短短的时间内星体可能已经运动过了很多很多位置。还有就是时间步长太大，一运行星体就飞的无影无踪了，所以这个步长要自己调。可以在界面上添加一个Slider控制步长，在下一次的版本中对这个进行优化。
Star1 = Star(1);%创建星体1Star2 = Star(2);%创建星体2Star3 = Star(3);%创建星体3t = 0.01; %时间步长% 如果星体的质量不是零就把它们显示出来if(Star1.M ~= 0) Star1.Display() endif(Star2.M~=0) Star2.Display()endif(Star3.M~=0) Star3.Display()end%迭代计算while true %% 计算星体间的相对距离 R12 = RelativeDistance(Star1.XYZ , Star2.XYZ) . R ; R13 = RelativeDistance(Star1.XYZ , Star3.XYZ) . R ; R23 = RelativeDistance(Star2.XYZ , Star3.XYZ) . R ; %% 先判断星球会不会碰撞爆炸；如果爆炸则对质量归零，把星体从坐标系中删除 %if ( (R12 <= max(Star1.R , Star2.R)) && (Star1.M~=0) && (Star2.M~=0) ) %碰撞距离大了一会就结束... if ( (R12 <= 10) && (Star1.M~=0) && (Star2.M~=0) ) Star1.M = 0; Star2.M = 0; Star1.H_Star = []; Star2.H_Star = []; end if ( (R13 <= 10) && (Star1.M~=0) && (Star3.M~=0) ) Star1.M = 0; Star3.M = 0; Star1.H_Star = []; Star3.H_Star = []; end if ( (R23 <= 10) && (Star2.M~=0) && (Star3.M~=0) ) Star2.M = 0; Star3.M = 0; Star2.H_Star = []; Star3.H_Star = []; end %% 计算万有引力 f12 = Gravitation(Star1.M , Star2.M , R12 , Star1.XYZ , Star2.XYZ) . f; f23 = Gravitation(Star2.M , Star3.M , R23 , Star2.XYZ , Star3.XYZ) . f; f31 = Gravitation(Star3.M , Star1.M , R13 , Star3.XYZ , Star1.XYZ) . f; %% 更新，计算星体的运动参数 Star1.MotionParameter(f12 , f23 , f31 , 1 , t) Star2.MotionParameter(f12 , f23 , f31 , 2 , t) Star3.MotionParameter(f12 , f23 , f31 , 3 , t) %% 先把上一次的星体清除 delete(Star1.H_Star) delete(Star2.H_Star) delete(Star3.H_Star) %% 更新星体位置 if Star1.M ~= 0 Star1.Display() end if Star2.M ~= 0 Star2.Display() end if Star3.M ~= 0 Star3.Display() end drawnow %如果停止运动这个toggle被按下了，就结束运动 if handles.StopMotion.Value == 1 return endend
可以十分清楚地看到，使用面向对象后，程序变得十分清楚。
至于OpeningFunction中的初始化的那些代码和其他控件的代码这里就不放了。
五. 不同参数的运行结果
需要说明一下，下面所有输入的参数都是自己看着给的，没算过。
1.默认的初始参数，运动差不多在水平面内
2.不知道怎么描述，自己看把
Emm 我就放这两个把，其他的，你们自己写程序去试一试。由于这个我还有不少的地方没有进行优化，所以今天就不把这个程序传上去了。下次我优化好了把这个程序放到百度云里，然后发个链接。
做完这个这个模拟后，我想起了流浪地球里开始说要利用木星来给地球加速。真的，设置了适当的参数后，你会发现小行星绕大行星半圈后直接甩出去了...。真的很直观。
六. 后记
这篇文章我其实主要是想说说MATLAB的面向对象编程。也建议大家以后编程都最好使用面向对象，只是面向对象写起来可能没有面向过程方便，但是使用和维护起来面向对象编程还是挺好的。我一个多月前学了点MATLAB中面向对象的语法，然后一直也没用过，这次正好拿来练练手，感觉和Python中的差不多。
这个程序我当时写的挺快的，只是大部分时间都在调试上了。我当时是一股脑地从头写到尾，写完虽然没有语法错误，也能运行，但是运行后的轨迹有问题。主要时间都花在这上面了。
以后应该每写一段程序就运行一下，看看结果是否是预期的。而不应该从头写到尾，然后再调试。
距离复试的日子也就十几天了，这是这个月的最后一篇推文，如果录取的话，下个月继续更，不然就要开始找工地搬砖喽也就不会有太多时间更新了，希望能过把，加油。
最后，我在这段时间经常感到遗憾，当时考数学时要是再仔细点那道积分的大题就可能不会从一开始就错，然后就有140多了，那样我跟某高打的赌就可能不会输，就差4分这个赌我就赢了....还有那英语阅读也不知道是为啥和平时练习的差距简直差了一半...哎，以后再加油吧。
