Matlab：秦始皇算法-三体微分方程数值解

伪程序猿l S x

已于 2023-02-01 21:18:19 修改

阅读量680

点赞数 5

分类专栏：数学算法文章标签：算法 matlab 线性代数矩阵

于 2023-01-29 00:34:01 首次发布

本文链接：https://blog.csdn.net/lsxhhh/article/details/128782716

版权

数学算法专栏收录该内容

4 篇文章 2 订阅

订阅专栏

三体人生活在三个恒星统治下的极端世界，为了生活的正常进行，他们需要得到一份精准的万年历，根据这份万年历指导人们的生产生活，而三颗太阳的运行微分方程无法得到解析解，因此在文明发展到牛顿时代时，三体人利用人列计算机求解，得到了三颗太阳的近似运动轨迹。

假设三颗太阳的质量为：

$M_1,M_2,M_3$

在半人马座α星附近找到一个空间中的固定点建立坐标系，利用位移矢量表示这三个恒星的位置：

$\left\{\begin{matrix} \overrightarrow{r_1}=(x_1,y_1,z_1)^T \\ \overrightarrow{r_2}=(x_2,y_2,z_2)^T \\ \overrightarrow{r_3}=(x_3,y_3,z_3)^T \end{matrix}\right.$

由于他们生活的行星质量比起三颗太阳的质量而言极小，因此可以略去。由牛顿万有引力定律，在空间中两个天体之间的万有引力为：

$F_1=\frac{Gm_1 m_2}{\begin{vmatrix} \overrightarrow{r_{2}}-\overrightarrow{r_{1}} \end{vmatrix}^3}(\overrightarrow{r_{2}}-\overrightarrow{r_{1}})$

$F_2=\frac{Gm_1 m_2}{\begin{vmatrix} \overrightarrow{r_{1}}-\overrightarrow{r_{2}} \end{vmatrix}^3}(\overrightarrow{r_{1}}-\overrightarrow{r_{2}})$

其中r1和r2分别是两个天体的位移矢量，m1和m2分别是两个天体的质量(单位：kg)，G为万有引力常量，其取值为：G=6.67×10-11N·m²/kg²。由牛顿第二定律：

$\frac{d^2\overrightarrow{r}}{dt^2}=\sum_{i}^{}\overrightarrow{F_i^e}$

对于其中任意一个太阳，其受到的外力均来自其他两个太阳对它的引力。于是对于这三个太阳可以列出如下微分方程组：

$\left\{\begin{matrix} \frac{d^2\overrightarrow{r_1}}{dt^2}=\frac{GM_1 M_2}{\begin{vmatrix} \overrightarrow{r_{2}}-\overrightarrow{r_{1}} \end{vmatrix}^3}(\overrightarrow{r_{2}}-\overrightarrow{r_{1}})+ \frac{GM_1 M_3}{\begin{vmatrix} \overrightarrow{r_{3}}-\overrightarrow{r_{1}} \end{vmatrix}^3}(\overrightarrow{r_{3}}-\overrightarrow{r_{1}}) \\ \frac{d^2\overrightarrow{r_2}}{dt^2}=\frac{GM_2 M_1}{\begin{vmatrix} \overrightarrow{r_{1}}-\overrightarrow{r_{2}} \end{vmatrix}^3}(\overrightarrow{r_{1}}-\overrightarrow{r_{2}})+ \frac{GM_2 M_3}{\begin{vmatrix} \overrightarrow{r_{3}}-\overrightarrow{r_{2}} \end{vmatrix}^3}(\overrightarrow{r_{3}}-\overrightarrow{r_{2}}) \\ \frac{d^2\overrightarrow{r_3}}{dt^2}=\frac{GM_3 M_1}{\begin{vmatrix} \overrightarrow{r_{1}}-\overrightarrow{r_{3}} \end{vmatrix}^3}(\overrightarrow{r_{1}}-\overrightarrow{r_{3}})+ \frac{GM_3 M_2}{\begin{vmatrix} \overrightarrow{r_{2}}-\overrightarrow{r_{3}} \end{vmatrix}^3}(\overrightarrow{r_{2}}-\overrightarrow{r_{3}}) \end{matrix}\right.$

并给出初始条件：

$\left\{\begin{matrix} \frac{\mathrm{d} \overrightarrow{r_i}}{\mathrm{d} t}=\overrightarrow{v_{i0}} \\ \overrightarrow{r_{i}}=\overrightarrow{r_{i0}} \end{matrix}\right.\ \,i=1,2,3$

该微分方程组中有三个向量微分方程方程，可以在笛卡尔坐标系中拆分为9个标量微分方程。该微分方程的解应当为连续的空间点集，可以表示为：

$\begin{pmatrix} \overrightarrow{r_1(t)} \\ \overrightarrow{r_2(t)} \\ \overrightarrow{r_3(t)} \end{pmatrix}\ \,t\in (0,+\infty )$

然而上述方程已经被证明无法得到上述解析解，或者说无法写出解的表达式。为了求解上述微分方程，三体人只能使用数值解，即将上述方程的解离散化：

$\begin{pmatrix} \overrightarrow{r_1^i} \\ \overrightarrow{r_2^i} \\ \overrightarrow{r_3^i} \end{pmatrix}\ \,i\in \begin{Bmatrix} 0,\Delta t,2\Delta t,...,N\Lambda t \end{Bmatrix}$

其中的Δt为两个解空间点之间的时间步长，为了达到更加精确的结果，这一数值应该取的足够小，因为当求解的时间步长Δt趋于0时，离散方程的解就是原方程的解。三体人利用这一形式的解，跳过了复杂的微分运算，将上述方程转化为离散的代数方程。为了用离散解表示方程中的微分运算，这里使用最简单的欧拉向前差分形式，即可以将位移矢量的微分近似为如下格式：

$\frac{d\overrightarrow{r}}{dt}(t=i\Delta t) \approx \frac{\overrightarrow{r_{i+1}}-\overrightarrow{r_{i}}}{\Delta t}$

$\frac{d^2\overrightarrow{r}}{dt^2}(t=i\Delta t)\approx \frac{\overrightarrow{r_{i+2}}-2\overrightarrow{r_{i+1}}+\overrightarrow{r_i}}{\Delta t^2}$

利用上述公式将微分方程离散化，例如对于第3个等式，微分方程方程可以化为：

$\frac{\overrightarrow{r_{3}^{i+2}}-2\overrightarrow{r_3^{i+1}}+\overrightarrow{r_3^i}}{\Delta t^2}=\frac{GM_3 M_1}{\begin{vmatrix} \overrightarrow{r_{1}^i}-\overrightarrow{r_{3}^i} \end{vmatrix}^3}(\overrightarrow{r_{1}^i}-\overrightarrow{r_{3}^i})+ \frac{GM_3 M_2}{\begin{vmatrix} \overrightarrow{r_{2}^i}-\overrightarrow{r_{3}^i} \end{vmatrix}^3}(\overrightarrow{r_{2}^i}-\overrightarrow{r_{3}^i})$

上述方程中，唯有i+2项即：

$\overrightarrow{r_{3}^{i+2}}$

为未知，i+2项可以直接由第i+1和第i项计算得出，因此该种差分方式也被叫做显式差分:

$\overrightarrow{r_{3}^{i+2}}=\Delta t^2(\frac{GM_3 M_1}{\begin{vmatrix} \overrightarrow{r_{1}^i}-\overrightarrow{r_{3}^i} \end{vmatrix}^3}(\overrightarrow{r_{1}^i}-\overrightarrow{r_{3}^i})+ \frac{GM_3 M_2}{\begin{vmatrix} \overrightarrow{r_{2}^i}-\overrightarrow{r_{3}^i} \end{vmatrix}^3}(\overrightarrow{r_{2}^i}-\overrightarrow{r_{3}^i}))+2\overrightarrow{r_3^{i+1}}-\overrightarrow{r_3^{i}}$

将第一、第二个等式用同样的方式离散化，可以写出该微分方程的代数方程组，利用计算机迭代求解，为了表述方便，将上述计算简写为以下形式：

$\begin{pmatrix} \overrightarrow{r_1^{i+2}} \\ \overrightarrow{r_2^{i+2}} \\ \overrightarrow{r_3^{i+2}} \end{pmatrix}=\begin{pmatrix} F_1(\overrightarrow{r_1^i},\overrightarrow{r_1^{i+1}}) \\ F_2(\overrightarrow{r_2^i},\overrightarrow{r_2^{i+1}}) \\ F_3(\overrightarrow{r_3^i},\overrightarrow{r_3^{i+1}}) \end{pmatrix}\ \ i=1,2,...$

例如其中函数F3表达式如下：

$F_3(\overrightarrow{r_3^i},\overrightarrow{r_3^{i+1}})=\Delta t^2(\frac{GM_3 M_1}{\begin{vmatrix} \overrightarrow{r_{1}^i}-\overrightarrow{r_{3}^i} \end{vmatrix}^3}(\overrightarrow{r_{1}^i}-\overrightarrow{r_{3}^i})+ \frac{GM_3 M_2}{\begin{vmatrix} \overrightarrow{r_{2}^i}-\overrightarrow{r_{3}^i} \end{vmatrix}^3}(\overrightarrow{r_{2}^i}-\overrightarrow{r_{3}^i}))+2\overrightarrow{r_3^{i+1}}-\overrightarrow{r_3^{i}}$

其他两个函数表达式可以此类推。于是三体人进行迭代，从i=1开始，首先利用初始速度向量得到i=2时三个天体的位置，而后的迭代都可以得到i+2次的解，依次进行下去，最终求出了该方程的数值解。利用人列计算机求解上述代数方程，首先通过天文学家的观测和计算，三颗太阳运行的初始条件为：

$\overrightarrow{r_{10}}=\begin{pmatrix} 5\times 10^{14} \\ 0.1 \\ 0.5 \end{pmatrix} \ \ ,\ \ \overrightarrow{v_{10}}=\begin{pmatrix} -90000 \\ -600 \\ -66 \end{pmatrix}$

$\overrightarrow{r_{20}}=\begin{pmatrix} 0.3 \\ 15\times\10^{14} \\ 50000 \end{pmatrix} \ \ ,\ \ \overrightarrow{v_{20}}=\begin{pmatrix} 900 \\ -2000 \\ 0 \end{pmatrix}$

$\overrightarrow{r_{10}}=\begin{pmatrix} 0.2 \\ 0.4 \\ 11\times\10^{14} \end{pmatrix} \ \ ,\ \ \overrightarrow{v_{10}}=\begin{pmatrix} -600 \\ -600 \\ 0.5 \end{pmatrix}$

卡文迪许测得三颗太阳的质量为：

$M_1=1.8\times\10^{27}\ \ ,\ \ M_2 = 1.6\times\10^{27}\ \ ,\ \ M_3=1.2\times\10^{27}$

而后秦始皇宣布人列计算机启动，运行如下代码(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的运行规律：