三体人生活在三个恒星统治下的极端世界,为了生活的正常进行,他们需要得到一份精准的万年历,根据这份万年历指导人们的生产生活,而三颗太阳的运行微分方程无法得到解析解,因此在文明发展到牛顿时代时,三体人利用人列计算机求解,得到了三颗太阳的近似运动轨迹。
假设三颗太阳的质量为:
在半人马座α星附近找到一个空间中的固定点建立坐标系,利用位移矢量表示这三个恒星的位置:
由于他们生活的行星质量比起三颗太阳的质量而言极小,因此可以略去。由牛顿万有引力定律,在空间中两个天体之间的万有引力为:
其中r1和r2分别是两个天体的位移矢量,m1和m2分别是两个天体的质量(单位:kg),G为万有引力常量,其取值为:G=6.67×10-11N·m²/kg²。由牛顿第二定律:
对于其中任意一个太阳,其受到的外力均来自其他两个太阳对它的引力。于是对于这三个太阳可以列出如下微分方程组:
并给出初始条件:
该微分方程组中有三个向量微分方程方程,可以在笛卡尔坐标系中拆分为9个标量微分方程。该微分方程的解应当为连续的空间点集,可以表示为:
然而上述方程已经被证明无法得到上述解析解,或者说无法写出解的表达式。为了求解上述微分方程,三体人只能使用数值解,即将上述方程的解离散化:
其中的Δt为两个解空间点之间的时间步长,为了达到更加精确的结果,这一数值应该取的足够小,因为当求解的时间步长Δt趋于0时,离散方程的解就是原方程的解。三体人利用这一形式的解,跳过了复杂的微分运算,将上述方程转化为离散的代数方程。为了用离散解表示方程中的微分运算,这里使用最简单的欧拉向前差分形式,即可以将位移矢量的微分近似为如下格式:
利用上述公式将微分方程离散化,例如对于第3个等式,微分方程方程可以化为:
上述方程中,唯有i+2项即:
为未知,i+2项可以直接由第i+1和第i项计算得出,因此该种差分方式也被叫做显式差分:
将第一、第二个等式用同样的方式离散化,可以写出该微分方程的代数方程组,利用计算机迭代求解,为了表述方便,将上述计算简写为以下形式:
例如其中函数F3表达式如下:
其他两个函数表达式可以此类推。于是三体人进行迭代,从i=1开始,首先利用初始速度向量得到i=2时三个天体的位置,而后的迭代都可以得到i+2次的解,依次进行下去,最终求出了该方程的数值解。利用人列计算机求解上述代数方程,首先通过天文学家的观测和计算,三颗太阳运行的初始条件为:
卡文迪许测得三颗太阳的质量为:
而后秦始皇宣布人列计算机启动,运行如下代码(MATLAB),人列计算机迭代求解,得到了三颗太阳最近的时间段内的近似运行规律:
clc
clear
%参数:
r10 = [5*10^14;0.1;0.5]; r20 = [0.3;15*10^14;50000]; r30 = [0.2;0.4;11*10^14];%位失
m1 = 1.8*10^27; m2 = 1.6*10^27 ; m3 = 1.2*10^27;%质量
v10 = [-90000;-600;-66]; v20 = [900;-2000;0]; v30 = [-600;-600;0.5];%速度矢量
G = 6.67*10^-11;
max_time = 15;
dt = 0.0001;
r1(1,:) = transpose(r10); r2(1,:) = transpose(r20); r3(1,:) = transpose(r30);
r1(2,:) = transpose(r10) + dt * transpose(v10);
r2(2,:) = transpose(r20) + dt * transpose(v20);
r3(2,:) = transpose(r30) + dt * transpose(v30);
for i=1:max_time/dt
r1_i2 = get_i_2(transpose(r1(i,:)),transpose(r1(i+1,:)),transpose(r2(i,:)),transpose(r3(i,:)),G * m1 * m2,G * m1 * m3,dt);
r2_i2 = get_i_2(transpose(r2(i,:)),transpose(r2(i+1,:)),transpose(r1(i,:)),transpose(r3(i,:)),G * m1 * m2,G * m2 * m3,dt);
r3_i2 = get_i_2(transpose(r3(i,:)),transpose(r3(i+1,:)),transpose(r1(i,:)),transpose(r2(i,:)),G * m1 * m3,G * m2 * m3,dt);
r1(i+2,:) = transpose(r1_i2);
r2(i+2,:) = transpose(r2_i2);
r3(i+2,:) = transpose(r3_i2);
end
t=1:max_time/dt+2;
r1_x = r1(:,1);r1_y = r1(:,2);r1_z = r1(:,3);
r2_x = r2(:,1);r2_y = r2(:,2);r2_z = r2(:,3);
r3_x = r3(:,1);r3_y = r3(:,2);r3_z = r1(:,3);
% for i = 1:1000:max_time/dt -1
% plot_x = [r1(i,1) r2(i,1) r3(i,1)];
% plot_y = [r1(i,2) r2(i,2) r3(i,2)];
% plot_z = [r1(i,3) r2(i,3) r3(i,3)];
% comet3(plot_x,plot_y,plot_z,0.1)
% pause(0.00001);
% end
plot(t,r3_x,'-b',LineWidth = 2)
hold on
plot(t,r3_y,'-r',LineWidth = 2)
hold on
plot(t,r3_z,'-k',LineWidth = 2)
legend({'r3x','r3y','r3z'},'Location','northwest');
function output_i_2 = get_i_2(this_i,this_i_1,other1_i,other2_i,K1,K2,dt)
%GET_I_2 此处显示有关此函数的摘要
% 此处显示详细说明
dec_other1_this = other1_i - this_i;
dec_other2_this = other2_i - this_i;
length_other1_this = vector_length(dec_other1_this);
length_other2_this = vector_length(dec_other2_this);
output_i_2 = (dt^2)*K1*dec_other1_this./(length_other1_this^3)+...
(dt^2)*K2*dec_other2_this./(length_other2_this^3) - this_i + 2 * this_i_1;
end
function output_ = vector_length(inputArg1)
%VECTOR_LENGTH 此处显示有关此函数的摘要
% 此处显示详细说明
output_ = sqrt(inputArg1(1,1)^2 + inputArg1(2,1)^2 + inputArg1(3,1)^2);
end
运行后的结果如下:
太阳1的运行规律:
太阳2的运行规律:
太阳3的运行规律: