Störmer-Verlet格式研究三维多体问题

clear;

clc;

close all

h = 0.5; %步长

t = 0:100000; %时间范围

N = length(t)/h; %迭代次数

G = 2.95912208286e-4; %万有引力常数

%初始位置

q1 = [0, 0, 0]’;

q2 = [-3.5023653, -3.8169847, -1.5507963]’;

q3= [9.0155314, -3.0458353, -1.6483708]’;

q4 = [8.3101420, -16.2901086, -7.2521278]’;

q5 = [11.4707666, -25.7294829, -10.8169456]’;

q6 = [-15.5387357, -25.2225594, -3.1902382]’;

q = [q1 q2 q3 q4 q5 q6];

%初始质量

m1 = 1; %太阳的质量

m2 = 0.000954786104043; %Jupiter木星的质量

m3 = 0.000285583733151; %Saturn土星的质量

m4 = 0.0000437273164546; %Uranus天王星的质量

m5 = 0.0000517759138449; %Neptune海王星的质量

m6 = 1/(1.3*10^8); %Pluto冥王星的质量

m = [m1 m2 m3 m4 m5 m6];

%初始速度

v1 = [0, 0, 0]’;

v2 = [0.00565429,-0.00412490,-0.00190589]’;

v3 = [0.00168318,0.00483525,0.00192462]’;

v4 = [0.00354178,0.00137102,0.00055209]’;

v5 = [0.0028930,0.00114527,0.00039677]’;

v6 = [0.00276725,-0.00170702,-0.00136504]’;

v = [v1 v2 v3 v4 v5 v6];

%动量

p1 = m1 * v1;

p2 = m2 * v2;

p3 = m3 * v3;

p4 = m4 * v4;

p5 = m5 * v5;

p6 = m6 * v6;

p = [p1 p2 p3 p4 p5 p6];

%万有引力

f1 = -G*(m1*m2*(q1-q2)/norm(q1-q2)^3+m1*m3*(q1-q3)/norm(q1-q3)^3+…

m1*m4*(q1-q4)/norm(q1-q4)^3+m1*m5*(q1-q5)/norm(q1-q5)^3+m1*m6*(q1-q6)/norm(q1-q6)^3);

f2 = -G*(m2*m1*(q2-q1)/norm(q2-q1)^3+m2*m3*(q2-q3)/norm(q2-q3)^3+…

m2*m4*(q2-q4)/norm(q2-q4)^3+m2*m5*(q2-q5)/norm(q2-q5)^3+m2*m6*(q2-q6)/norm(q2-q6)^3);

f3 = -G*(m3*m1*(q3-q1)/norm(q3-q1)^3+m3*m2*(q3-q2)/norm(q3-q2)^3+…

m3*m4*(q3-q4)/norm(q3-q4)^3+m3*m5*(q3-q5)/norm(q3-q5)^3+m3*m6*(q3-q6)/norm(q3-q6)^3);

f4 = -G*(m4*m1*(q4-q1)/norm(q4-q1)^3+m4*m2*(q4-q2)/norm(q4-q2)^3+…

m4*m3*(q4-q3)/norm(q4-q3)^3+m4*m5*(q4-q5)/norm(q4-q5)^3+m4*m6*(q4-q6)/norm(q4-q6)^3);

f5 = -G*(m5*m1*(q5-q1)/norm(q5-q1)^3+m5*m2*(q5-q2)/norm(q5-q2)^3+…

m5*m3*(q5-q3)/norm(q5-q3)^3+m5*m4*(q5-q4)/norm(q5-q4)^3+m5*m6*(q5-q6)/norm(q5-q6)^3);

f6 = -G*(m6*m1*(q6-q1)/norm(q6-q1)^3+m6*m2*(q6-q2)/norm(q6-q2)^3+…

m6*m3*(q6-q3)/norm(q6-q3)^3+m6*m4*(q6-q4)/norm(q6-q4)^3+m6*m5*(q6-q5)/norm(q6-q5)^3);

f = [f1 f2 f3 f4 f5 f6];

Time=[];H=[];

x1=[]; y1=[]; z1=[];

x2=[]; y2=[]; z2=[];

x3=[]; y3=[]; z3=[];

x4=[]; y4=[]; z4=[];

x5=[]; y5=[]; z5=[];

x6=[]; y6=[]; z6=[];

for k= 1:N-1 %Störmer-Verlet格式

x1=[x1 q1(1)]; y1=[y1 q1(2)]; z1=[z1 q1(3)];

x2=[x2 q2(1)]; y2=[y2 q2(2)]; z2=[z2 q2(3)];

x3=[x3 q3(1)]; y3=[y3 q3(2)]; z3=[z3 q3(3)];

x4=[x4 q4(1)]; y4=[y4 q4(2)]; z4=[z4 q4(3)];

x5=[x5 q5(1)]; y5=[y5 q5(2)]; z5=[z5 q5(3)];

x6=[x6 q6(1)]; y6=[y6 q6(2)]; z6=[z6 q6(3)];

p1 = p1 + h*f1/2;

p2 = p2 + h*f2/2;

p3 = p3 + h*f3/2;

p4 = p4 + h*f4/2;

p5 = p5 + h*f5/2;

p6 = p6 + h*f6/2;

q1 = q1 + h*p1/m1;

q2 = q2 + h*p2/m2;

q3 = q3 + h*p3/m3;

q4 = q4 + h*p4/m4;

q5 = q5 + h*p5/m5;

q6 = q6 + h*p6/m6;