blog
圆形限制性三体问题中的共线平动点位置
可能很多人对三体问题感兴趣是源于小说《三体》。其实三体问题是轨道力学中特别经典的问题， 近年随着深空探测的兴起，有不少航天任务利用了三体问题里的一些特殊轨道。如嫦娥四号任务中，中继卫星鹊桥就是利用了在L2的Halo轨道形成对月球背面的信号中继。 平动点问题也成为普通大众所了解的话题。
在经典意义下，三体是不可解的。数学上就会把三体问题做一些特殊处理，其中较为常见的一种方式就是所谓的圆形限制性三体问题。 可惜圆形限制性三体问题也没有一般解析结果，但是找到了一些特殊解。这里就简单介绍共线平动点的计算。
无量纲单位下，在会合坐标中圆形限制性三体问题运动方程为 \[\left\{ \begin{gathered} \ddot x - 2\dot y = - \frac{{\partial \overline U }}{{\partial x}} \hfill \\ \ddot y + 2\dot x = - \frac{{\partial \overline U }}{{\partial y}} \hfill \\ \ddot z = - \frac{{\partial \overline U }}{{\partial z}} \hfill \\ \end{gathered} \right.\] 其三个共线平动点满足以下方程 \[\left\{ \begin{gathered} x - \frac{{1 - \mu }}{{{{\left( {x + \mu } \right)}^2}}} + \frac{\mu }{{{{\left( {x - 1 + \mu } \right)}^2}}} = 0 \hfill \\ x - \frac{{1 - \mu }}{{{{\left( {x + \mu } \right)}^2}}} - \frac{\mu }{{{{\left( {x - 1 + \mu } \right)}^2}}} = 0 \hfill \\ x + \frac{{1 - \mu }}{{{{\left( {x + \mu } \right)}^2}}} + \frac{\mu }{{{{\left( {x - 1 + \mu } \right)}^2}}} = 0 \hfill \\ \end{gathered} \right.\] 以上方程有近似解析结果分别为 \[{\xi ^{\left( 1 \right)}}{\text{ = }}{\left( {\frac{\mu }{3}} \right)^{\frac{1}{3}}}\left[ {1{\text{ - }}\frac{1}{3}{{\left( {\frac{\mu }{3}} \right)}^{\frac{1}{3}}}{\text{ - }}\frac{1}{9}{{\left( {\frac{\mu }{3}} \right)}^{\frac{2}{3}}}{\text{ + }} \cdots } \right]\] \[{\xi ^{\left( 2 \right)}}{\text{ = }}{\left( {\frac{\mu }{3}} \right)^{\frac{1}{3}}}\left[ {1{\text{ + }}\frac{1}{3}{{\left( {\frac{\mu }{3}} \right)}^{\frac{1}{3}}}{\text{ - }}\frac{1}{9}{{\left( {\frac{\mu }{3}} \right)}^{\frac{2}{3}}}{\text{ + }} \cdots } \right]\] \[\left\{ \begin{gathered} {\xi ^{\left( 3 \right)}}{\text{ = }}1{\text{ - }}\nu \left[ {1{\text{ + }}\frac{{23}}{{84}}{\nu ^2}{\text{ + }}\frac{{23}}{{84}}{\nu ^3}{\text{ + }}\frac{{761}}{{2352}}{\nu ^4}{\text{ + }}\frac{{3163}}{{7056}}{\nu ^5}{\text{ + }} \cdots } \right] \hfill \\ \nu {\text{ = }}\frac{7}{{12}}\mu \hfill \\ \end{gathered} \right.\] 解析法精度较低，实际计算应该用数值法直接求解。下面给出代码：
package main
import (
"fmt"
"math"
)
func main() {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-18 10:29
! -----------------------------------------------------
! Input Parameters :
!
! Output Parameters :
!
------------------------------------------------------*/
u := 0.01215057
_, _, _ = LagrangianPoint(u)
}
func LagrangianPoint(u float64) (x1 float64, x2 float64, x3 float64) {
/*------------------------------------------------------
Created : Song Yezhi 2022-5-5 21:35
计算圆形限制性三体问题三个共线平动点位置
--------------------------------------------------------
Input Parameters :
Output Parameters :
-------------------------------------------------------*/
var xa, xb float64
xa, xb = L1X0(u)
fmt.Printf("L1 \n")
x1, _ = secant(xa, xb, u, L1fx)
fmt.Printf("L2 \n")
xa, xb = L2X0(u)
x2, _ = secant(xa, xb, u, L2fx)
fmt.Printf("L3 \n")
xa, xb = L3X0(u)
x3, _ = secant(xa, xb, u, L3fx)
//fmt.Printf("%15.8f %15.9f %15.9f %15.9f \n", u , x1, x2,x3)
return
}
func L1X0(u float64) (x0 float64, x1 float64) {
/*------------------------------------------------------
Created : Song Yezhi 2022-5-5 21:24
L1初值级数解
--------------------------------------------------------
Input Parameters :
u -----
Output Parameters :
两个初值
-------------------------------------------------------*/
x0 = math.Pow(u/3.0, 1.0/3.0)
x1 = x0 - 1.0/3.0*x0 - 1.0/9*x0*x0
return
}
func L2X0(u float64) (x0 float64, x1 float64) {
/*------------------------------------------------------
Created : Song Yezhi 2022-5-5 21:26
L2初值级数解
--------------------------------------------------------
Input Parameters :
Output Parameters :
-------------------------------------------------------*/
x0 = math.Pow(u/3.0, 1.0/3.0)
x1 = x0 + 1.0/3.0*x0 - 1.0/9*x0*x0
return
}
func L3X0(u float64) (x0 float64, x1 float64) {
/*------------------------------------------------------
Created : Song Yezhi 2022-5-5 21:26
L2初值级数解
--------------------------------------------------------
Input Parameters :
Output Parameters :
-------------------------------------------------------*/
niu := 7.0 / 12.0 * u
x0 = 1.0 - niu
x1 = 1.0 - niu - 23.0/84.0*math.Pow(niu, 3) - 23.0/84.0*math.Pow(niu, 4)
return
}
func L1fx(x float64, u float64) float64 {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-18 10:28
! L1点函数
------------------------------------------------------*/
tmp1 := (x - u) * (x - u)
tmp2 := (x + 1.0 - u) * (x + 1.0 - u)
fx := x + (1.0-u)/tmp1 - u/tmp2
return fx
}
func L2fx(x float64, u float64) float64 {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-18 10:31
! L1点函数
------------------------------------------------------*/
tmp1 := (x - u) * (x - u)
tmp2 := (x + 1.0 - u) * (x + 1.0 - u)
fx := x + (1.0-u)/tmp1 + u/tmp2
return fx
}
func L3fx(x float64, u float64) float64 {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-18 10:31
! L1点函数
------------------------------------------------------*/
tmp1 := (x - u) * (x - u)
tmp2 := (x + 1.0 - u) * (x + 1.0 - u)
fx := x - (1.0-u)/tmp1 - u/tmp2
return fx
}
func secant(x0 float64, x1 float64, u float64, funcX func(float64, float64) float64) (x, fx float64) {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-18 10:30
! secant方法求解方程
! 注意：只计算一次方程函数
! -----------------------------------------------------
! Input Parameters :
! x0-----初值
! funcX------目标函数
! Output Parameters :
!
-------------------------------------------------------*/
var imax int = 30
var tol float64 = 1e-8
var x2 float64
var fx0, fx1 float64
fx0 = funcX(x0, u)
fx1 = funcX(x1, u)
for i := 0; i < imax; i++ {
x2 = x1 - (x1-x0)/(fx1-fx0)*fx1
dx := math.Abs(x1 - x0)
if dx < tol {
break
}
fx0 = fx1
fx1 = funcX(x2, u)
x0 = x1
x1 = x2
fmt.Printf("i= %4d x= %12.7f f(x)=%12.7f \n", i, x1, fx)
}
x = x1
fx = fx1
return x, fx
}
对于地月系而言，其中 \(\mu {\text{ = 0}}{\text{.01215057}}\)，输出结果为
L1 i= 0 x= 0.1944830 f(x)= 0.0000000 ... i= 22 x= -0.8369574 f(x)= 0.0000000 i= 23 x= -0.8369151 f(x)= 0.0000000 i= 24 x= -0.8369152 f(x)= 0.0000000 i= 25 x= -0.8369152 f(x)= 0.0000000 L2 i= 0 x= 0.2732924 f(x)= 0.0000000 ... i= 23 x= -1.1556889 f(x)= 0.0000000 i= 24 x= -1.1556821 f(x)= 0.0000000 i= 25 x= -1.1556821 f(x)= 0.0000000 i= 26 x= -1.1556821 f(x)= 0.0000000 L3 i= 0 x= 1.0049123 f(x)= 0.0000000 i= 1 x= 1.0050608 f(x)= 0.0000000 i= 2 x= 1.0050626 f(x)= 0.0000000 i= 3 x= 1.0050626 f(x)= 0.0000000其中x就表示三个共线平动点在无量纲单位下的位置。