ラグランジュ点（２）

（ラグランジュ点（１）の続き）

地球−月系のラグランジュ点の位置を実際に計算してみる。$x \equiv \frac{r}{R}$，$a \equiv \frac{M_2}{M_1}$ とすると，次の式の解がそれぞれ $\rm{L}_1, \rm{L}_2, \rm{L}_3$ を与える。
\begin{equation*}
\begin{aligned}
\frac{1}{(1-x)^2}=\frac{a}{x^2}+ 1 -x*(1+a)\\
\frac{1}{(1-x)^2}+\frac{a}{x^2}= 1 +x*(1+a)\\
\frac{1}{(1-x)^2}+\frac{a}{(2-x)^2}= a +(1-x)*(1+a)
\end{aligned}
\end{equation*}
これは$x$の５次方程式であるが，Mathematicaならば，NSolve[]でもよいけれど，複素数解は不要なので，FindRoot[]を使うところだ。

Juliaではどうするのか調べてみたところ，NLsolveという非線形方程式の求解パッケージがあったが，GitHubのオリジナルサイトをみてもJacobian がなんたらかんたらでと五月蝿くて分かりにくい。阪大サイバーメディアセンターの降旗大介さんのページがたいへん親切な解説をつけてくれていたので参考にする。

using NLsolve

function nls(func, params...; ini = [0.0])
    if typeof(ini) <: Number
        return nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini]).zero[1]
    else
        return nlsolve((vout,vin)->vout .= func(vin,params...), ini).zero
    end
end

function lagrange(x,a,i)
    if i==1
        1.0/(1.0-x)^2-a/(x^2)-1.0+x*(1+a)
    elseif i==2
        1.0/(1.0+x)^2+a/(x^2)-1.0-x*(1+a)
    elseif i==3
        1.0/(1.0-x)^2+a/(2-x)^2-a-(1-x)*(1+a)
    end
end

# earth - moon system  M2/M1=0.0123
L1=nls(lagrange, 0.0123, 1, ini=0.1)
L2=nls(lagrange, 0.0123, 2, ini=0.1)
L3=nls(lagrange, 0.0123, 3, ini=0.1)
println([L1,L2,L3]," ",[L1,L2,L3]*38.4)

# sun - earth system  M2/M1=0.000003
S1=nls(lagrange, 0.000003, 1, ini=0.1)
S2=nls(lagrange, 0.000003, 2, ini=0.1)
S3=nls(lagrange, 0.000003, 3, ini=0.1)
println([S1,S2,S3]," ",[S1,S2,S3]*14960)


[0.150934, 0.167833, 0.00708792] [5.79587, 6.44477, 0.272176]
[0.00996655, 0.0100332, 1.74999e-6] [149.1, 150.097, 0.0261799]

ということで，無事に鵲橋のハロー軌道の元になる地球−月系の$\rm{L}_2$点の6.4万kmが得られた。ここでは，ラグランジュ点の安定性の議論はすっとばしているが，京大の三宅雄紀さんの「ラグランジュ点の安定性を調べてみた」もたいへん参考になる。

「わが十指われにかしずく寒の入り」（岡本眸　1928-2018）

（ラグランジュ点（３）に続く）