レーン・エムデン方程式(Lane-Emden equation)をExcel VBAでルンゲクッタ8次で計算。: Fallen Physicist, Rising Engineer

$\frac{1}{\xi^2} \frac{d}{d\xi} \left({\xi^2 \frac{d\theta}{d\xi}}\right) + \theta^n = 0$

というような形をしている。ここでθ(0)=1, dθ(0)/dξ=0。

これはξ=0で発散するような形の式になるので数値計算しにくい。

そこで、θをξのべきで展開して、ちょっとξ=0から離れた位置から計算する、ということはよくやられている。

たいてい、ξの4次の項くらいまで計算されているが、、、いつもこのブログではルンゲクッタ8次のDormand&Princeを使っていることから！

今回はMaximaさんの力も借りて、8次まで計算してみよう。

θ(ξ)= c₀+c₁ξ+c₂ξ²+c₃ξ³+c₄ξ⁴+c₅ξ⁵+c₆ξ⁶+c₇ξ⁷+c₈ξ⁸+...

で、c0=1, c1=0, c2=1/6, c3=0, c4 = n/120, c5=0,c6=-n(n/1890-1/3024), c7=0

c8=-(n ^ 3 - 3 * n ^ 2 + 2 * n) * c2 ^ 3 + 6 * c4 * n * (n - 1) * c2 + 3 * c3 ^ 2 * n * n + (-3 * c3 ^ 2 + 6 * c6 * nn) / 432

ふう。こりゃ手計算ではしんどいのでMaximaさんさまさま。
これで、ξ=Δξを出発点にしてルンゲクッタ8次を使ってn=0～6まで計算。

厳密解がある、n=0,1,5とも比べてぴったり（上には描いてないけど）。

*というか、これだけべき展開するならもうルンゲクッタいらないんじゃ、、、

土

Fallen Physicist, Rising Engineer