matlab练习程序（三体模型）

最近在看三体电视剧，正好看到了计算三体数值解那一部分，就想起了上学时看三体，也用matlab实现了三体的运动模拟。

不过当时是通过时域外推的方式实现的，不是很严谨。

下面通过微分方程求解三体问题，三体模型的微分方程：

公式中x1是质点1的位置向量(px,py,pz)，x1''是质点1的加速度向量(ax,ay,az)。

微分方程列出来后，参考之前的文章，按套路实现即可。

matlab代码如下：

clear all;close all;clc;

G = 6.67259e-11;
m1 = rand() / G;
m2 = rand() / G;
m3 = rand() / G;

%x = [px1 py1 pz1 px2 py2 pz2 px3 py3 pz3 ...
%     vx1 vy1 vz1 vx2 vy2 vz2 vx3 vy3 vz3]
[t,x]=ode45(@(t,x) testfun(t,x,G,m1,m2,m3),0:0.001:10,rand(1,18)-0.5);

plot3(x(:,1),x(:,2),x(:,3),'r');
hold on;
plot3(x(:,4),x(:,5),x(:,6),'b');
plot3(x(:,7),x(:,8),x(:,9),'g');
grid on;

function dx=testfun(t,x,G,m1,m2,m3)

x1 = x(1:3);
x2 = x(4:6);
x3 = x(7:9);

v1 = x(10:12);
v2 = x(13:15);
v3 = x(16:18);

a1 = G*(m2*(x2-x1)/(norm(x2-x1))^3+m3*(x3-x1)/(norm(x3-x1))^3);
a2 = G*(m3*(x3-x2)/(norm(x3-x2))^3+m1*(x1-x2)/(norm(x1-x2))^3);
a3 = G*(m1*(x1-x3)/(norm(x1-x3))^3+m2*(x2-x3)/(norm(x2-x3))^3);

dx = [v1;v2;v3;a1;a2;a3];
end

结果如下：