song.yz@foxmail.com wechat: math-box

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就表示三个共线平动点在无量纲单位下的位置。