弱重力場中の光の経路の近似解：別解法

シュバルツシルト時空中の光の経路を決める式を2階微分方程式の形にして解く。

光の経路を決める式

$$\begin{array}{rcl}{\left(\frac{du}{d\varphi }\right)}^2& =& \frac{1}{b^2}-u^2+r_g(u^3-\frac{1}{b^3})\end{array}$$

両辺を $\varphi$ で微分して2階微分方程式の形にする。

$$\begin{array}{rcl}2\frac{du}{d\varphi }\frac{d^{2}u}{d{\varphi }^2}& =& -2u\frac{du}{d\varphi }+3r_g{u}^2\frac{du}{d\varphi }\\ \therefore \frac{d^{2}u}{d{\varphi }^2}& =& -u+\frac{3}{2}{r}_g{u}^2\end{array}$$

初期条件

1階微分方程式のときには1個の積分定数を決めるための初期条件はひとつでよくて，${\displaystyle \varphi =\frac{\pi }{2}}$ のとき，${\displaystyle u=\frac{1}{b}}$ とした。2階微分方程式の場合は，2個の積分定数を決めるために初期条件もふたつ必要となるから，もとの1階微分方程式の形から以下を採用する。

• ${\displaystyle \varphi =\frac{\pi }{2}}$ のとき，${\displaystyle u=\frac{1}{b}}$ および，
• ${\displaystyle \varphi =\frac{\pi }{2}}$ のとき，${\displaystyle \frac{du}{d\varphi }=0}$

$r_g$ のゼロ次解

$r_g$ がかかった項を無視したときの解を $u_0$ とおくと

$$\begin{array}{rcl}\frac{d^2{u}_0}{d{\varphi }^2}& =& -u_0\end{array}$$

初期条件 ${\displaystyle \varphi =\frac{\pi }{2}}$ のとき，${\displaystyle u_0=\frac{1}{b}}$，${\displaystyle \frac{d{u}_0}{d\varphi }=0}$ より

$$u_0=\frac{\sin \varphi }{b}$$

$r_g$ の1次解

${\displaystyle u=u_0+\frac{r_g}{b}{u}_1}$ とおいて微分方程式に代入し，$r_g$ の1次の項を取り出すと

$$\frac{d^2{u}_1}{d{\varphi }^2}+u_1=\frac{3}{2b}{\sin }^2\varphi$$

この非同次方程式の解（同次方程式の基本解と非同次方程式の特殊解の和）に初期条件を課す。

$u$ 全体に対する初期条件は ${\displaystyle \varphi =\frac{\pi }{2}}$ のとき，${\displaystyle u=\frac{1}{b}}$，${\displaystyle \frac{du}{d\varphi }=0}$ であり，これは $u_1$ に対しては，${\displaystyle u_1=0,\frac{d{u}_1}{d\varphi }=0}$ となるから，解は

$$u_1=\frac{1}{2b}(2-\sin \varphi -{\sin }^2\varphi )$$

したがって $r_g$ の1次までの解は

$$u=u_0+\frac{r_g}{b}{u}_1=\frac{\sin \varphi }{b}+\frac{r_g}{2b^2}(2-\sin \varphi -{\sin }^2\varphi )$$

eq: 'diff(u1, phi, 2) + u1 = 3/(2*b) * sin(phi)**2;

sol: ode2(eq, u1, phi);

sol1: ic2(sol, phi = %pi/2, u1 = 0, 'diff(u1, phi) = 0);

ev(sol1, cos(2*phi)=1 - 2*sin(phi)**2), ratsimp;

from sympy.abc import *
from sympy import *

u1 = Function('u1')

eq = Eq(Derivative(u1(phi), phi, 2) + u1(phi), 
        Rational(3,2) * sin(phi)**2/b
)
eq

dsolve(eq, u1(phi), 
       ics = {u1(phi).subs(phi, pi/2):0,          # u1(pi/2) = 0
              diff(u1(phi), phi).subs(phi, pi/2):0 # d u1/d phi(pi/2) = 0
             })

factor(_)

_.subs(cos(phi)**2, 1 - sin(phi)**2)