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

光の経路を決める式

\begin{eqnarray}
\left(\frac{du}{d\phi}\right)^2
&=& \frac{1}{b^2} -u^2 + r_g \left(\,u^3 -\frac{1}{b^3}\right)
\end{eqnarray}

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

\begin{eqnarray}
2\frac{du}{d\phi}\,\frac{d^2u}{d\phi^2}
&=& -2 u\, \frac{du}{d\phi}+ 3 r_g \,u^2 \, \frac{du}{d\phi} \\
\therefore\ \ \frac{d^2u}{d\phi^2} &=& -u + \frac{3}{2} r_g \, u^2
\end{eqnarray}

初期条件

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

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

$r_g$ のゼロ次解

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

\begin{eqnarray}
\frac{d^2 u_0}{d\phi^2} &=& -u_0
\end{eqnarray}

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

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

$r_g$ の1次解

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

$$\frac{d^2 u_1}{d\phi^2} + u_1 = \frac{3}{2 b} \sin^2 \phi$$

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

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

$$ u_1 = \frac{1}{2 b} \left( 2 -\sin \phi -\sin^2 \phi \right)$$

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

$$u = u_0 + \frac{r_g}{b} u_1 = \frac{\sin\phi}{b} + \frac{r_g}{2 b^2} \left( 2 -\sin\phi -\sin^2\phi\right)$$

Maxima で非同次方程式を解いてみる

In [1]:

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

Out[1]:

\[\tag{${\it \%o}_{1}$}\frac{d^2}{d\,\varphi^2}\,u_{1}+u_{1}=\frac{3\,\sin ^2\varphi}{2\,b}\]

In [2]:

sol: ode2(eq, u1, phi);

Out[2]:

\[\tag{${\it \%o}_{2}$}u_{1}=\frac{\cos \left(2\,\varphi\right)+3}{4\,b}+{\it \%k}_{1}\,\sin \varphi+{\it \%k}_{2}\,\cos \varphi\]

初期条件 $\displaystyle \phi = \frac{\pi}{2}$ で $\displaystyle u_1 = 0, \frac{d u_1}{d\phi} = 0$ を課すと…

In [3]:

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

Out[3]:

\[\tag{${\it \%o}_{3}$}u_{1}=\frac{\cos \left(2\,\varphi\right)+3}{4\,b}-\frac{\sin \varphi}{2\,b}\]

$\cos 2 \phi \Rightarrow 1 – 2 \sin^2 \phi$ と置き換えて…

In [4]:

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

Out[4]:

\[\tag{${\it \%o}_{4}$}u_{1}=-\frac{\sin ^2\varphi+\sin \varphi-2}{2\,b}\]

SymPy で非同次方程式を解いてみる

モジュールの import

In [1]:

from sympy.abc import *
from sympy import *

In [2]:

u1 = Function('u1')

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

Out[2]:

$\displaystyle u_{1}{\left(\phi \right)} + \frac{d^{2}}{d \phi^{2}} u_{1}{\left(\phi \right)} = \frac{3 \sin^{2}{\left(\phi \right)}}{2 b}$

In [3]:

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
             })

Out[3]:

$\displaystyle u_{1}{\left(\phi \right)} = – \frac{\sin{\left(\phi \right)}}{2 b} + \frac{\cos^{2}{\left(\phi \right)}}{2 b} + \frac{1}{2 b}$

In [4]:

factor(_)

Out[4]:

$\displaystyle u_{1}{\left(\phi \right)} = \frac{- \sin{\left(\phi \right)} + \cos^{2}{\left(\phi \right)} + 1}{2 b}$

In [5]:

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

Out[5]:

$\displaystyle u_{1}{\left(\phi \right)} = \frac{- \sin^{2}{\left(\phi \right)} – \sin{\left(\phi \right)} + 2}{2 b}$

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