多体运动的matlab动画演示

多体运动的matlab动画演示

1 问题说明

当考虑的体系中的对象超过两个时（比如三个），由于其相互作用的复杂性，使得多体动力学问题变得极其复杂，要用解析的办法通过求解动力学微分方程来求得多体系统的每个对象的运动状态几乎是不可能的。

不过，如果只是为了获得多体系统的粗略而简要的认知，那么可以利用matlab等软件在数值上求解多体系统的动力学微分方程，从而给出多体系统的大致运动状况。

需要说明的是，matlab求解动力学微分方程所得到的结果毕竟是存在误差的，这些误差主要来自于其算法的迭代过程，即舍入误差和截断误差。并且随着迭代次数增加，所产生的误差会不断累积，使得求解的时间尺度越大，后面的误差也就越大。所以要想获得比较可靠的计算结果，则所求解的时间尺度不能太大。

另外，该文档中所采用的动画演示方法是先用ode45在初始条件下求解出所有粒子的运动状态，包括速度和位置，然后按照所求解出来的位置，在每个对应的时间下在图中画出该粒子，从而形成可以连续演示的动画。

由于不是在求解微分方程的同时给出粒子的运动，这种方法演示更加流畅；但是依然能看出在长时间演示后期，动画演示也会迟滞，这是需要改进的地方。

2 理论求解

在这里，所考虑的多体系统可以包括3个粒子，4个粒子甚至更多粒子，只要计算机足够强大，可以放入100个粒子也行。所考虑的粒子分为两类，一种是以质子为原型，一种以中子为原型。

所考虑的作用力包括了两种：强相互作用力，电磁相互作用力。

在程序中求解时，需要把两种作用力加起来（具体见求解微分方程的m文件）。

这里只处理了平面中的运动情况，也可以将其更改为立体空间中的运动，那样会更加复杂。

程序已经写为多体（N体）的通用形式，改变N的大小时，不需要再去更改求解微分方程的m文件的内容，而只需要在参数设置区更改相应的粒子数以及每个粒子的具体信息。

注意，Rx,Ry,m,ke这几个代表着每个粒子信息的量，这几个变量所包含的数值个数（向量维度）应该等于N（n+p）的值，否则将出现不匹配的情况而无法运行。

比如：

n=3;p=3;

Rx=[0,1,1,1,1/2,0];Ry=[0,0,1/2,1,1/2,1];

m=[1.0012,1.0012,1,1.0012,1.0012,1];

ke=[0,0,1,0,1,0];%Rx,Ry,m,ke都应该包含6个数值，注意理清6个粒子对应的参量。

n=3;p=2;

Rx=[0,1/2,1/2,1/2,1];Ry=[1/2,0,1/2,1,1/2];

m=[1,1.0012,1.0012,1.0012,1];

ke=[1,0,0,0,1];%Rx,Ry,m,ke都应该包含5个数值

3 matlab程序

注：如果要copy该段程序直接放入matlab，需要调整注释（绿色的字体部分）

文件1

figure('name','多体运动演示'); %设置标题名字

global N m ke r0 %定义全局变量，使得求解微分方程的m文件可以使用这些变量

%****************************参数设置区********************************%

n=3;p=3; %n,p分别为中子和质子数

Rx=[0,1,1,1,1/2,0];Ry=[0,0,1/2,1,1/2,1]; %Rx和Ry分别为起始位置的坐标

m=[1.0012,1.0012,1,1.0012,1.0012,1]; %m以质子的质量为单位1的相对质量值

ke=[0,0,1,0,1,0]; %ke为以e为单位的电荷值

%*********************************************************************%

N=n+p;r0=0.4; %N为所有的粒子总数，r0为坐标尺度相对值，可调整

pausetime=.01;%设置暂停时间

set(gca,'xlim',[-2 2],'ylim',[-2 2]);%设置图形窗口的坐标显示范围（可根据实际情况进行更改。）

set(gcf,'doublebuffer','on') %消除抖动

axis equal

hold on

x0=zeros(4*N,1); %x0为求解方程组的初始值

for i=1:N

x0(4*i-2)=Rx(i);x0(4*i)=Ry(i);

if ke(i)==1

pp(i)=plot(x0(4*i-2),x0(4*i),'color','r','marker','.','markersize',15); %p为红色点

else

pp(i)=plot(x0(4*i-2),x0(4*i),'color','k','marker','.','markersize',15); %n为黑色点

end

end

t0=0;tf=15; %求解的时间范围

[t,x]=ode45('dohezi',t0,tf,x0); %调用m文件求解微分方程

len=length(t);

for i=1:len %在图中作出运动状况for j=1:N

set(pp(j),'xdata',x(i,4*j-2),'ydata',x(i,4*j));

plot(x(i,4*j-2),x(i,4*j));

end

pause(pausetime); %暂停一会

drawnow

end

下面是求解微分方程时调用的m文件的内容：

注：如果要copy该段程序进matlab，其文件名需要命名为dohezi.m保存，同时注意调整注释（绿色部分）。

文件2

function solhe=dohezi(t,x)

global N m ke r0

solhe=zeros(4*N,1); %x(4*j-3),x(4*j-2),x(4*j-1),x(4*j)分别为第j个粒子的vx,x,vy,y；solhe(4*j-3),solhe(4*j-2),solhe(4*j-1),solhe(4*j)分别对应vx,x,vy,y四个值对时间t的导数；

for j=1:N %j对应每个粒子

solhe(4*j-3)=0;solhe(4*j-1)=0;

for k=1:N

if j~=k %第j个粒子受到除了第j个粒子之外的其他所有粒子的作用力，所以要把其他粒子施加的作用全部求和。这里的作用力又分为两种：强相互作用的力和电磁相互作用的力。

r=((x(4*j-2)-x(4*k-2))^2+(x(4*j)-x(4*k))^2)^0.5; %r为第j个粒子和第k个粒子的相对距离

solhe(4*j-3)=solhe(4*j-3)+15/m(j)/r^2*(1/r0-1/r)*exp(-r/r0)*(x(4*k-2)-x (4*j-2))+ke(j)*ke(k)/m(j)/137*(x(4*j-2)-x(4*k-2))/r^3;

solhe(4*j-1)=solhe(4*j-1)+15/m(j)/r^2*(1/r0-1/r)*exp(-r/r0)*(x(4*k)-x(4 *j))+ke(j)*ke(k)/m(j)/137*(x(4*j)-x(4*k))/r^3;

else

end

end

solhe(4*j-2)=x(4*j-3);

solhe(4*j)=x(4*j-1);

end

4 运行结果

（1）三个粒子的运动

其中有两个p粒子（红色），一个n粒子（黑色）。从运动轨迹可以看出，三体运动是极其复杂的。三个粒子的运动相互交错融合，毫无规律可循。

（2）四个粒子的运动

其中有两个p粒子（红色），两个n粒子（黑色）。由于初始位置的对称，在开始一段时间四个粒子的运动几乎各自表现为简谐振动，但后来由于求解微分方程不断累积误差，使得各个粒子的轨迹不断偏离理论值，最后呈现出杂乱无章的状态。

（3）六个粒子的运动

其中有两个p粒子（红色），四个n粒子（黑色）。六个粒子的运动状态最初很有规律，但随着时间增加，运动轨迹的复杂度急速增加，最后变成一团毛线。

5 结语

该文档中的程序只处理了平面中的运动情况，而且设计并非完美，只希望起到抛砖引玉的作用，希望读者能根据该内容，做出更加完美的多体运动的matlab动画演示，得到更加准确的数值解，或者能在更大的时间范围内求解数值解。

更多内容见文库内本人的其他matlab 动画演示。