程序模拟三体问题

傷城~ 2023-10-14 23:29 7阅读 0赞

#### 文章目录 ####

*  1. 拉格朗日函数
 *  2. 拉格朗日方程的带入
 *  3. 数 学 工 具！！
 *  4. 结果的呈现
 *  5. 尾声

注：此篇文章考虑的是二维平面上的三题运动。使用字母上的点表示关于时间的导数，例如：  
 x ˙ = ∂ x ∂ t x ¨ = ∂ 2 x ∂ t 2 \\dot\{x\}=\\frac\{\\partial x\}\{\\partial t\} \\\\ \\ddot\{x\}=\\frac\{\\partial ^2x\}\{\\partial t^2\} x˙=∂t∂xx¨=∂t2∂2x

## 1. 拉格朗日函数 ##

采用笛卡尔坐标系，三个质点的质量分别为  m 1 , m 2 , m 3 m\_1,m\_2,m\_3 m1,m2,m3，坐标分别为  ( x 1 , y 1 ) ， ( x 2 , y 2 ) ， ( x 3 , y 3 ) (x1,y1) ， (x2,y2) ， (x3,y3) (x1,y1)，(x2,y2)，(x3,y3)，质点的速度可以表述为  x ˙ 2 + y ˙ 2 \\sqrt\{\\dot\{x\}^2+\\dot\{y\}^2\} x˙2\+y˙2，便可以写出三个质点的动能：  
 T 1 = 1 2 m 1 ( x 1 ˙ 2 + y 1 ˙ 2 ) T 2 = 1 2 m 2 ( x 2 ˙ 2 + y 2 ˙ 2 ) T 3 = 1 2 m 3 ( x 3 ˙ 2 + y 3 ˙ 2 ) T\_1=\\frac\{1\}\{2\}m\_1(\\dot\{x\_1\}^2+\\dot\{y\_1\}^2) \\\\ T\_2=\\frac\{1\}\{2\}m\_2(\\dot\{x\_2\}^2+\\dot\{y\_2\}^2) \\\\ T\_3=\\frac\{1\}\{2\}m\_3(\\dot\{x\_3\}^2+\\dot\{y\_3\}^2) T1=21m1(x1˙2\+y1˙2)T2=21m2(x2˙2\+y2˙2)T3=21m3(x3˙2\+y3˙2)  
考虑系统的势能，此处的势能为引力势能，可写为：  
 E p = − G m 1 ∗ m 2 r E\_p=-G\\frac\{m\_1\*m\_2\}\{r\}\\\\ Ep=−Grm1∗m2  
其中： G G G为引力常数，  m 1 m\_1 m1 ，  m 2 m\_2 m2为相互作用的两个质点的质量， r 为两质点间的距离。  G , m 1 , m 2 G,m\_1,m\_2 G,m1,m2 作为常数带入，棘手的是两质点间的距离  r r r  
 r a , b = ( x a − x b ) 2 + ( y a − y b ) 2 r\_\{a,b\}=\\sqrt\{(x\_a-x\_b)^2+(y\_a-y\_b)^2\}\\\\ ra,b=(xa−xb)2\+(ya−yb)2  
 a , b a,b a,b为两个质点，在三体问题中可取值  1 , 2 , 3 1,2,3 1,2,3 且  a ≠ b a\\ne b a=b。势能为：  
 U = − G m 1 ∗ m 2 r 1 , 2 − G m 1 ∗ m 3 r 1 , 3 − G m 2 ∗ m 3 r 2 , 3 U=-G\\frac\{m\_1\*m\_2\}\{r\_\{1,2\}\}-G\\frac\{m\_1\*m\_3\}\{r\_\{1,3\}\}-G\\frac\{m\_2\*m\_3\}\{r\_\{2,3\}\}\\\\ U=−Gr1,2m1∗m2−Gr1,3m1∗m3−Gr2,3m2∗m3  
整合以上内容，带入  L = T − U L=T-U L=T−U ，即得：  
 L = T 1 + T 2 + T 3 − U L=T\_1+T\_2+T\_3-U\\\\ L=T1\+T2\+T3−U  
展开式过于复杂，这里不做展开。

## 2. 拉格朗日方程的带入 ##

我们求得了拉格朗日函数 L 现在我们将其带入拉格朗日方程即：  
 d d t ( ∂ L ∂ q α ˙ ) − ∂ L ∂ q α = 0 ( α = 0 , 1 , 2 , . . . , s ) \\frac\{d\}\{dt\}(\\frac\{\\partial L\}\{\\partial \\dot\{q\_\\alpha\}\})-\\frac\{\\partial L\}\{\\partial q\_\\alpha\}=0\\ \\ (\\alpha=0,1,2,...,s)\\\\ dtd(∂qα˙∂L)−∂qα∂L=0(α=0,1,2,...,s)  
其中  q α q\_\\alpha qα为系统的广义坐标，  s s s 为系统的自由度，此处  s s s为  6 6 6 ，即可得到  6 6 6 个微分方程，确定初始条件后便可以解除数值解。三题问题目前没有解析解，三题问题的解析解也是数学上的一个大课题。

好在的是，我们的有强大的数学工具  M a t h e m a t i c a Mathematica Mathematica

## 3. 数 学 工 具！！ ##

注：接下来的阅读需要一定的  M a t h e m a t i c a Mathematica Mathematica基础，这儿的 x 1 , y 1 x1,y1 x1,y1对应  ( x 1 , y 1 ) (x\_1,y\_1) (x1,y1)  
 d i s 12 = S q r t \[ ( x 1 \[ t \] − x 2 \[ t \] ) 2 + ( y 1 \[ t \] − y 2 \[ t \] ) 2 \] ; d i s 13 = S q r t \[ ( x 1 \[ t \] − x 3 \[ t \] ) 2 + ( y 1 \[ t \] − y 3 \[ t \] ) 2 \] ; d i s 23 = S q r t \[ ( x 2 \[ t \] − x 3 \[ t \] ) 2 + ( y 2 \[ t \] − y 3 \[ t \] ) 2 \] ; v 1 = S q r t \[ x 1 ′ \[ t \] 2 + y 1 ′ \[ t \] 2 \] ; v 2 = S q r t \[ x 2 ′ \[ t \] 2 + y 2 ′ \[ t \] 2 \] ; v 3 = S q r t \[ x 3 ′ \[ t \] 2 + y 3 ′ \[ t \] 2 \] ; dis12 = Sqrt\[(x1\[t\] - x2\[t\])^2 + (y1\[t\] - y2\[t\])^2\]; dis13 = Sqrt\[(x1\[t\] - x3\[t\])^2 + (y1\[t\] - y3\[t\])^2\]; dis23 = Sqrt\[(x2\[t\] - x3\[t\])^2 + (y2\[t\] - y3\[t\])^2\]; v1 = Sqrt\[x1'\[t\]^2 + y1'\[t\]^2\]; v2 = Sqrt\[x2'\[t\]^2 + y2'\[t\]^2\]; v3 = Sqrt\[x3'\[t\]^2 + y3'\[t\]^2\]; dis12=Sqrt\[(x1\[t\]−x2\[t\])2\+(y1\[t\]−y2\[t\])2\];dis13=Sqrt\[(x1\[t\]−x3\[t\])2\+(y1\[t\]−y3\[t\])2\];dis23=Sqrt\[(x2\[t\]−x3\[t\])2\+(y2\[t\]−y3\[t\])2\];v1=Sqrt\[x1′\[t\]2\+y1′\[t\]2\];v2=Sqrt\[x2′\[t\]2\+y2′\[t\]2\];v3=Sqrt\[x3′\[t\]2\+y3′\[t\]2\];  
为了接下来的书写方便，我们先定义 r 和 v ，例如dis12对应 r\_\{1,2\}

我们便写出 Mathematica 下的拉格朗日函数 L :  
 L = 1 / 2 ∗ m 1 ∗ v 1 2 + 1 / 2 ∗ m 2 ∗ v 2 2 + 1 / 2 ∗ m 3 ∗ v 3 2 + G ∗ ( m 1 ∗ m 2 ) / d i s 12 + G ∗ ( m 1 ∗ m 3 ) / d i s 13 + G ∗ ( m 2 ∗ m 3 ) / d i s 23 ; L = 1/2\*m1\*v1^2 + 1/2\*m2\*v2^2 + 1/2\*m3\*v3^2 + G\*(m1\*m2)/dis12 + G\*(m1\*m3)/dis13 + G\*(m2\*m3)/dis23; L=1/2∗m1∗v12\+1/2∗m2∗v22\+1/2∗m3∗v32\+G∗(m1∗m2)/dis12\+G∗(m1∗m3)/dis13\+G∗(m2∗m3)/dis23;  
同样的，我们写出拉格朗日方程：

LEquation = Function[D[D[#1, #3], t] - D[#1, #2]];

我们写作函数形式，方便下文调用，接下来我们使用这个函数写出六个方程

Equations = {
        
       LEquation[L, x1[t], x1'[t]] == 0,
       LEquation[L, y1[t], y1'[t]] == 0,
       LEquation[L, x2[t], x2'[t]] == 0,
       LEquation[L, y2[t], y2'[t]] == 0,
       LEquation[L, x3[t], x3'[t]] == 0,
       LEquation[L, y3[t], y3'[t]] == 0
    };

之后我们写出参数和基础条件

parameters = {
        G -> 6.672*10^-11, m1 -> 10^12, m2 -> 2*10^12, 
       m3 -> 3*10^12};
    init = {
        
       x1[0] == 0, y1[0] == 0, (*以下三行（包括本行）为三个质点的初始位置*)
       
       x2[0] == 0, y2[0] == 0,
       x3[0] == 0, y3[0] == 0,
       
       (*以下三行为三个质点的初速度*)
       x1'[0] == 0, y1'[0] == 0,
       x2'[0] == 0, y2'[0] == 0,
       x3'[0] == 0, y3'[0] == 0
       };

有了以上这些玩意儿，我们便可以开始快乐的求解，善良的 Mathematica 给我们提供了NDSolve方法，此方法可以求出给定微分方程（组）的数值解，具体如下：

solution = 
     NDSolve[{
        Equations, init} /. parameters, {
        x1, y1, x2, y2, x3, 
        y3}, {
        t, 0, 50}, 
       Method -> 
        "IndexReduction" -> {
        "Pantelides", 
          "ConstraintMethod" -> "Projection"}][[1]]

（吐槽一下CSDN的代码排版）t为求解的时间范围（可改动），跟上的行参Method，那么大一坨，可以参考官方附带的Reference，毕竟，我也是从别人那学来的操作awa。解出来啥样呢？

我们看到，三个坐标的（对应每一个时间的）值已经求解出来了，接下来便是结果的呈现

## 4. 结果的呈现 ##

POWERFUL的 Mathematica 提供了Graphics方法以供我们作图（更多详细请参考官方的Reference，这里展开不久偏题了吗qwq），配合Manipulate，便可以作出完整的运动轨迹。

Manipulate[
     Graphics[
      {
        
        PointSize[0.015],
        {
        Black, Point[{
        x1[t], y1[t]}]},
        {
        Blue, Point[{
        x2[t], y2[t]}]},
        {
        Red, Point[{
        x3[t], y3[t]}]},
        {
        Black, 
         Line[Map[Function[Evaluate[{
        x1[#], y1[#]} /. solution]], 
           Range[0, t, 0.025]]]},
        {
        Blue, 
         Line[Map[Function[Evaluate[{
        x2[#], y2[#]} /. solution]], 
           Range[0, t, 0.025]]]},
        {
        Red, 
         Line[Map[Function[Evaluate[{
        x3[#], y3[#]} /. solution]], 
           Range[0, t, 0.025]]]}
        } /. solution,
      Axes -> {
        True, True},
      PlotRange -> {
        {
        -100, 100}, {
        -100, 100}},
      Ticks -> False
      ],
     {
        t, 0, 50}
     ]

此处的t可改动，最好与上文NDSolve中的t吻合

我们按动：Shift+Enter，啪唧！！

之后，跟着我的步伐，按动进度条条旁边的加号，点击播放键…

是不是很美妙呐！是不是呐。这里再附上两个案例：

## 5. 尾声 ##

通读全文，便可发现本文作者完全没有技术含量，拉格朗日力学的解题也十分的规范

希望本文能够给到你帮助，有错误也欢迎指出

Appendix 牛顿力学的求解  
我们采用笛卡尔坐标系，三个质点的坐标为 ( x 1 , y 1 ) (x1,y1) (x1,y1)，  ( x 2 , y 2 ) (x2,y2) (x2,y2)，  ( x 3 , y 3 ) (x3,y3) (x3,y3)

万有引力的大小可以写为：  
 F = G m 1 ∗ m 2 r 2 F=G\\frac\{m\_1\*m\_2\}\{r^2\}\\\\ F=Gr2m1∗m2  
其方向指向相互作用的质点。我们使用  F ⃗ a , b \\vec\{F\}\_\{a,b\} Fa,b表示质点 a,b 之间的引力，此处 a,b 可取值 1,2,3 且 a\\ne b ，我们便可以得出作用在每一个质点上的合力，这里以质点 1 举例：  
 F ⃗ 1 = F ⃗ 1 , 2 + F ⃗ 1 , 3 ∣ F ⃗ ∣ = F 1 , 2 2 + 2 F 1 , 2 ∗ F 1 , 3 ∗ cos ⁡ θ + F 1 , 3 2 \\vec\{F\}\_1=\\vec\{F\}\_\{1,2\}+\\vec\{F\}\_\{1,3\}\\\\ |\\vec\{F\}|=\\sqrt\{F\_\{1,2\}^2+2F\_\{1,2\}\*F\_\{1,3\}\*\\cos\\theta+F\_\{1,3\}^2\}\\\\ F1=F1,2\+F1,3∣F∣=F1,22\+2F1,2∗F1,3∗cosθ\+F1,32  
其中： θ \\theta θ为  F ⃗ 1 , 2 , F ⃗ 1 , 3 \\vec\{F\}\_\{1,2\},\\vec\{F\}\_\{1,3\} F1,2,F1,3之间的夹角，可表述为：  
 θ = π 2 − ( arctan ⁡ y 2 − y 1 x 2 − x 1 + arctan ⁡ y 3 − y 1 x 3 − x 1 ) ∣ F ⃗ 1 ∣ = F 1 , 2 2 + 2 F 1 , 2 ∗ F 1 , 3 ∗ sin ⁡ ( arctan ⁡ y 2 − y 1 x 2 − x 1 + arctan ⁡ y 3 − y 1 x 3 − x 1 ) + F 1 , 3 2 \\theta=\\frac\{\\pi\}\{2\}-(\\arctan\{\\frac\{y\_2-y\_1\}\{x\_2-x\_1\}\}+\\arctan\{\\frac\{y\_3-y\_1\}\{x\_3-x\_1\}\})\\\\ |\\vec\{F\}\_1|=\\sqrt\{F\_\{1,2\}^2+2F\_\{1,2\}\*F\_\{1,3\}\*\\sin(\\arctan\{\\frac\{y\_2-y\_1\}\{x\_2-x\_1\}\}+\\arctan\{\\frac\{y\_3-y\_1\}\{x\_3-x\_1\}\})+F\_\{1,3\}^2\}\\\\ θ=2π−(arctanx2−x1y2−y1\+arctanx3−x1y3−y1)∣F1∣=F1,22\+2F1,2∗F1,3∗sin(arctanx2−x1y2−y1\+arctanx3−x1y3−y1)\+F1,32  
我们把合力  F ⃗ 1 \\vec\{F\}\_\{1\} F1带入牛顿方程：  
 a = ∣ F ⃗ 1 ∣ m 1 a=\\frac\{|\\vec\{F\}\_1|\}\{m\_1\}\\\\ a=m1∣F1∣  
 a ⃗ \\vec\{a\} a 的方向与  F ⃗ 1 \\vec\{F\}\_1 F1的方向一致。后积分求出  x ⃗ ( t ) \\vec\{x\}(t) x(t) 即可。