Back Geodesics Simulation The metric and geodesic equations are in the form:
d s 2 = − ( 1 − 2 r ) d t 2 + ( 1 − 2 r ) − 1 d r 2 + r 2 d ϕ 2 , d 2 t d λ 2 + 2 r ( r − 2 ) d t d λ d r d λ = 0 , d 2 r d λ 2 + r − 2 r 3 ( d t d λ ) 2 − 1 r ( r − 2 ) ( d r d λ ) 2 − ( r − 2 ) ( d ϕ d λ ) 2 = 0 , d 2 ϕ d λ 2 + 2 r d r d λ d ϕ d λ = 0 , \begin{align*} ds^2 = -\left(1 - \frac{2}{r}\right) dt^2 + \left(1 - \frac{2}{r}\right)^{-1} dr^2 + r^2 d\phi^2, \\ \frac{d^2 t}{d\lambda^2} + \frac{2}{r(r - 2)} \frac{dt}{d\lambda} \frac{dr}{d\lambda} &= 0, \\ \frac{d^2 r}{d\lambda^2} + \frac{r - 2}{r^3} \left(\frac{dt}{d\lambda}\right)^2 - \frac{1}{r(r - 2)} \left(\frac{dr}{d\lambda}\right)^2 - (r - 2) \left(\frac{d\phi}{d\lambda}\right)^2 &= 0, \\ \frac{d^2 \phi}{d\lambda^2} + \frac{2}{r} \frac{dr}{d\lambda} \frac{d\phi}{d\lambda} &= 0, \end{align*} d s 2 = − ( 1 − r 2 ) d t 2 + ( 1 − r 2 ) − 1 d r 2 + r 2 d ϕ 2 , d λ 2 d 2 t + r ( r − 2 ) 2 d λ d t d λ d r d λ 2 d 2 r + r 3 r − 2 ( d λ d t ) 2 − r ( r − 2 ) 1 ( d λ d r ) 2 − ( r − 2 ) ( d λ d ϕ ) 2 d λ 2 d 2 ϕ + r 2 d λ d r d λ d ϕ = 0 , = 0 , = 0 , implying:
d 2 t d λ 2 = − 2 r ( r − 2 ) d t d λ d r d λ , d 2 r d λ 2 = − r − 2 r 3 ( d t d λ ) 2 + 1 r ( r − 2 ) ( d r d λ ) 2 + ( r − 2 ) ( d ϕ d λ ) 2 , d 2 ϕ d λ 2 = − 2 r d r d λ d ϕ d λ , \begin{align*} \frac{d^2 t}{d\lambda^2} &= - \frac{2}{r(r - 2)} \frac{dt}{d\lambda} \frac{dr}{d\lambda}, \\ \frac{d^2 r}{d\lambda^2} &= -\frac{r - 2}{r^3} \left(\frac{dt}{d\lambda}\right)^2 + \frac{1}{r(r - 2)} \left(\frac{dr}{d\lambda}\right)^2 + (r - 2) \left(\frac{d\phi}{d\lambda}\right)^2, \\ \frac{d^2 \phi}{d\lambda^2} &= -\frac{2}{r} \frac{dr}{d\lambda} \frac{d\phi}{d\lambda}, \end{align*} d λ 2 d 2 t d λ 2 d 2 r d λ 2 d 2 ϕ = − r ( r − 2 ) 2 d λ d t d λ d r , = − r 3 r − 2 ( d λ d t ) 2 + r ( r − 2 ) 1 ( d λ d r ) 2 + ( r − 2 ) ( d λ d ϕ ) 2 , = − r 2 d λ d r d λ d ϕ , which can be solved using the Euler method:
U n t = U n − 1 t − h 2 r n − 1 ( r n − 1 − 2 ) U n − 1 t U n − 1 r , U n r = U n − 1 r + h [ − r n − 1 − 2 ( r n − 1 ) 3 ( U n − 1 t ) 2 + 1 r n − 1 ( r n − 1 − 2 ) ( U n − 1 r ) 2 + ( r n − 1 − 2 ) ( U n − 1 ϕ ) 2 ] , U n ϕ = U n − 1 ϕ − h 2 r n − 1 U n − 1 r U n − 1 ϕ , t n = t n − 1 + h U n − 1 t , r n = r n − 1 + h U n − 1 r , ϕ n = ϕ n − 1 + h U n − 1 ϕ , \begin{align*} U^t_n &= U^t_{n - 1} - h \frac{2}{r_{n - 1}(r_{n - 1} - 2)} U^t_{n - 1} U^r_{n - 1}, \\ U^r_n &= U^r_{n - 1} + h \left[- \frac{r_{n - 1} - 2}{(r_{n - 1})^3} \left(U^t_{n - 1}\right)^2 + \frac{1}{r_{n - 1}(r_{n - 1} - 2)} \left(U^r_{n - 1}\right)^2 + (r_{n - 1} - 2) \left(U^{\phi}_{n - 1}\right)^2\right], \\ U^{\phi}_n &= U^{\phi}_{n - 1} - h \frac{2}{r_{n - 1}} U^r_{n - 1} U^{\phi}_{n - 1}, \\[4ex] t_n &= t_{n - 1} + h U^t_{n - 1}, \\ r_n &= r_{n - 1} + h U^r_{n - 1}, \\ \phi_n &= \phi_{n - 1} + h U^{\phi}_{n - 1}, \end{align*} U n t U n r U n ϕ t n r n ϕ n = U n − 1 t − h r n − 1 ( r n − 1 − 2 ) 2 U n − 1 t U n − 1 r , = U n − 1 r + h [ − ( r n − 1 ) 3 r n − 1 − 2 ( U n − 1 t ) 2 + r n − 1 ( r n − 1 − 2 ) 1 ( U n − 1 r ) 2 + ( r n − 1 − 2 ) ( U n − 1 ϕ ) 2 ] , = U n − 1 ϕ − h r n − 1 2 U n − 1 r U n − 1 ϕ , = t n − 1 + h U n − 1 t , = r n − 1 + h U n − 1 r , = ϕ n − 1 + h U n − 1 ϕ , where h h h U μ = d x μ d λ U^{\mu} = \frac{dx^{\mu}}{d\lambda} U μ = d λ d x μ
From the geodesics chapter , the length of a tangent vector is given by:
ϵ = − ( 1 − 2 r ) ( d t d λ ) 2 + ( 1 − 2 r ) − 1 ( d r d λ ) 2 + r 2 ( d ϕ d λ ) 2 , \epsilon = -\left(1 - \frac{2}{r}\right) \left(\frac{dt}{d\lambda}\right)^2 + \left(1 - \frac{2}{r}\right)^{-1} \left(\frac{dr}{d\lambda}\right)^2 + r^2 \left(\frac{d\phi}{d\lambda}\right)^2, ϵ = − ( 1 − r 2 ) ( d λ d t ) 2 + ( 1 − r 2 ) − 1 ( d λ d r ) 2 + r 2 ( d λ d ϕ ) 2 , where ϵ \epsilon ϵ
ϵ < 0 , for timelike geodesics , ϵ = − 1 , for timelike geodesics where λ = τ , ϵ = 0 , for lightlike geodesics , ϵ > 0 , for spacelike geodesics , ϵ = 1 , for spacelike geodesics where λ = L 0 . \begin{align*} \epsilon &< 0, & &\textrm{for timelike geodesics}, \\ \epsilon &= -1, & &\textrm{for timelike geodesics where \(\lambda = \tau\)}, \\ \epsilon &= 0, & &\textrm{for lightlike geodesics}, \\ \epsilon &> 0, & &\textrm{for spacelike geodesics}, \\ \epsilon &= 1, & &\textrm{for spacelike geodesics where \(\lambda = L_0\)}. \end{align*} ϵ ϵ ϵ ϵ ϵ < 0 , = − 1 , = 0 , > 0 , = 1 , for timelike geodesics , for timelike geodesics where λ = τ , for lightlike geodesics , for spacelike geodesics , for spacelike geodesics where λ = L 0 . For timelike paths, I will parametrize path by the proper time τ \tau τ ϵ = − 1 \epsilon = -1 ϵ = − 1
d 2 t d τ 2 = − 2 r ( r − 2 ) d t d τ d r d τ , d 2 r d τ 2 = − r − 2 r 3 ( d t d τ ) 2 + 1 r ( r − 2 ) ( d r d τ ) 2 + ( r − 2 ) ( d ϕ d τ ) 2 , d 2 ϕ d τ 2 = − 2 r d r d τ d ϕ d τ . \begin{align*} \frac{d^2 t}{d\tau^2} &= - \frac{2}{r(r - 2)} \frac{dt}{d\tau} \frac{dr}{d\tau}, \\ \frac{d^2 r}{d\tau^2} &= -\frac{r - 2}{r^3} \left(\frac{dt}{d\tau}\right)^2 + \frac{1}{r(r - 2)} \left(\frac{dr}{d\tau}\right)^2 + (r - 2) \left(\frac{d\phi}{d\tau}\right)^2, \\ \frac{d^2 \phi}{d\tau^2} &= -\frac{2}{r} \frac{dr}{d\tau} \frac{d\phi}{d\tau}. \end{align*} d τ 2 d 2 t d τ 2 d 2 r d τ 2 d 2 ϕ = − r ( r − 2 ) 2 d τ d t d τ d r , = − r 3 r − 2 ( d τ d t ) 2 + r ( r − 2 ) 1 ( d τ d r ) 2 + ( r − 2 ) ( d τ d ϕ ) 2 , = − r 2 d τ d r d τ d ϕ . The initial time dilation can be derived from the tangent vector length:
1 = ( 1 − 2 r ) ( U 0 t ) 2 − ( 1 − 2 r ) − 1 ( U 0 r ) 2 − r 2 ( U 0 ϕ ) 2 , r − 2 r ( U 0 t ) 2 = 1 + r r − 2 ( U 0 r ) 2 + r 2 ( U 0 ϕ ) 2 , U 0 t = r r − 2 [ 1 + r r − 2 ( U 0 r ) 2 + r 2 ( U 0 ϕ ) 2 ] , \begin{align*} 1 &= \left(1 - \frac{2}{r}\right) \left(U^t_0\right)^2 - \left(1 - \frac{2}{r}\right)^{-1} \left(U^r_0\right)^2 - r^2 \left(U^{\phi}_0\right)^2, \\ \frac{r - 2}{r} \left(U^t_0\right)^2 &= 1 + \frac{r}{r - 2} \left(U^r_0\right)^2 + r^2 \left(U^{\phi}_0\right)^2, \\ U^t_0 &= \sqrt{\frac{r}{r - 2}\left[1 + \frac{r}{r - 2} \left(U^r_0\right)^2 + r^2 \left(U^{\phi}_0\right)^2\right]}, \\ \end{align*} 1 r r − 2 ( U 0 t ) 2 U 0 t = ( 1 − r 2 ) ( U 0 t ) 2 − ( 1 − r 2 ) − 1 ( U 0 r ) 2 − r 2 ( U 0 ϕ ) 2 , = 1 + r − 2 r ( U 0 r ) 2 + r 2 ( U 0 ϕ ) 2 , = r − 2 r [ 1 + r − 2 r ( U 0 r ) 2 + r 2 ( U 0 ϕ ) 2 ] , then the Euler method can be used to solve for the motion.
For lightlike paths, ϵ = 0 \epsilon = 0 ϵ = 0
0 = − ( 1 − 2 r ) ( d t d λ ) 2 + ( 1 − 2 r ) − 1 ( d r d λ ) 2 + r 2 ( d ϕ d λ ) 2 . 0 = -\left(1 - \frac{2}{r}\right) \left(\frac{dt}{d\lambda}\right)^2 + \left(1 - \frac{2}{r}\right)^{-1} \left(\frac{dr}{d\lambda}\right)^2 + r^2 \left(\frac{d\phi}{d\lambda}\right)^2. 0 = − ( 1 − r 2 ) ( d λ d t ) 2 + ( 1 − r 2 ) − 1 ( d λ d r ) 2 + r 2 ( d λ d ϕ ) 2 . And the geodesic equations don't change:
d 2 t d λ 2 = − 2 r ( r − 2 ) d t d λ d r d λ , d 2 r d λ 2 = − r − 2 r 3 ( d t d λ ) 2 + 1 r ( r − 2 ) ( d r d λ ) 2 + ( r − 2 ) ( d ϕ d λ ) 2 , d 2 ϕ d λ 2 = − 2 r d r d λ d ϕ d λ . \begin{align*} \frac{d^2 t}{d\lambda^2} &= - \frac{2}{r(r - 2)} \frac{dt}{d\lambda} \frac{dr}{d\lambda}, \\ \frac{d^2 r}{d\lambda^2} &= -\frac{r - 2}{r^3} \left(\frac{dt}{d\lambda}\right)^2 + \frac{1}{r(r - 2)} \left(\frac{dr}{d\lambda}\right)^2 + (r - 2) \left(\frac{d\phi}{d\lambda}\right)^2, \\ \frac{d^2 \phi}{d\lambda^2} &= -\frac{2}{r} \frac{dr}{d\lambda} \frac{d\phi}{d\lambda}. \end{align*} d λ 2 d 2 t d λ 2 d 2 r d λ 2 d 2 ϕ = − r ( r − 2 ) 2 d λ d t d λ d r , = − r 3 r − 2 ( d λ d t ) 2 + r ( r − 2 ) 1 ( d λ d r ) 2 + ( r − 2 ) ( d λ d ϕ ) 2 , = − r 2 d λ d r d λ d ϕ . Consider a light source S S S r 0 r_0 r 0 M M M L L L ( r ~ , ϕ ~ ) (\tilde{r}, \tilde{\phi}) ( r ~ , ϕ ~ )
All right rays travel from S S S d r ~ d λ = c = 1 \frac{d\tilde{r}}{d\lambda} = c = 1 d λ d r ~ = c = 1 d ϕ ~ d λ = 0 \frac{d\tilde{\phi}}{d\lambda} = 0 d λ d ϕ ~ = 0 λ = 0 \lambda = 0 λ = 0
d r d λ = − cos ϕ ~ , v ⊥ = r 0 d ϕ d λ = sin ϕ ~ , d ϕ d λ = sin ϕ ~ r 0 . \begin{align*} \frac{dr}{d\lambda} = -\cos \tilde{\phi}, \\ v_{\perp} = r_0 \frac{d\phi}{d\lambda} = \sin \tilde{\phi}, \\ \frac{d\phi}{d\lambda} = \frac{\sin \tilde{\phi}}{r_0}. \end{align*} d λ d r = − cos ϕ ~ , v ⊥ = r 0 d λ d ϕ = sin ϕ ~ , d λ d ϕ = r 0 sin ϕ ~ . Note: ϕ ~ \tilde{\phi} ϕ ~ d r ~ d λ \frac{d\tilde{r}}{d\lambda} d λ d r ~ r 0 r_0 r 0
Substituting into the length of tangent vector formula (at λ = 0 , r = r 0 \lambda = 0, r = r_0 λ = 0 , r = r 0
0 = − ( 1 − 2 r 0 ) ( d t d λ ) 2 + ( 1 − 2 r 0 ) − 1 cos 2 ϕ ~ + r 0 2 sin 2 ϕ ~ r 0 2 , ( 1 − 2 r 0 ) ( d t d λ ) 2 = ( 1 − 2 r 0 ) − 1 cos 2 ϕ ~ + sin 2 ϕ ~ , ( d t d λ ) 2 = ( 1 − 2 r 0 ) − 1 [ ( 1 − 2 r 0 ) − 1 cos 2 ϕ ~ + sin 2 ϕ ~ ] , d t d λ = ( 1 − 2 r 0 ) − 1 [ ( 1 − 2 r 0 ) − 1 cos 2 ϕ ~ + sin 2 ϕ ~ ] . \begin{align*} 0 &= -\left(1 - \frac{2}{r_0}\right) \left(\frac{dt}{d\lambda}\right)^2 + \left(1 - \frac{2}{r_0}\right)^{-1} \cos^2 \tilde{\phi} + r_0^2 \frac{\sin^2 \tilde{\phi}}{r_0^2}, \\ \left(1 - \frac{2}{r_0}\right) \left(\frac{dt}{d\lambda}\right)^2 &= \left(1 - \frac{2}{r_0}\right)^{-1} \cos^2 \tilde{\phi} + \sin^2 \tilde{\phi}, \\ \left(\frac{dt}{d\lambda}\right)^2 &= \left(1 - \frac{2}{r_0}\right)^{-1} \left[\left(1 - \frac{2}{r_0}\right)^{-1} \cos^2 \tilde{\phi} + \sin^2 \tilde{\phi}\right], \\ \frac{dt}{d\lambda} &= \sqrt{\left(1 - \frac{2}{r_0}\right)^{-1} \left[\left(1 - \frac{2}{r_0}\right)^{-1} \cos^2 \tilde{\phi} + \sin^2 \tilde{\phi}\right]}. \end{align*} 0 ( 1 − r 0 2 ) ( d λ d t ) 2 ( d λ d t ) 2 d λ d t = − ( 1 − r 0 2 ) ( d λ d t ) 2 + ( 1 − r 0 2 ) − 1 cos 2 ϕ ~ + r 0 2 r 0 2 sin 2 ϕ ~ , = ( 1 − r 0 2 ) − 1 cos 2 ϕ ~ + sin 2 ϕ ~ , = ( 1 − r 0 2 ) − 1 [ ( 1 − r 0 2 ) − 1 cos 2 ϕ ~ + sin 2 ϕ ~ ] , = ( 1 − r 0 2 ) − 1 [ ( 1 − r 0 2 ) − 1 cos 2 ϕ ~ + sin 2 ϕ ~ ] . Now, the Euler method can be employed. But since we're working with multiple rays at once, I will use the coordinate time t t t λ \lambda λ
0 = − ( 1 − 2 r ) ( d t d λ ) 2 + ( 1 − 2 r ) − 1 ( d r d λ ) 2 + r 2 ( d ϕ d λ ) 2 , ( 1 − 2 r ) ( d t d λ ) 2 = ( 1 − 2 r ) − 1 ( d r d λ ) 2 + r 2 ( d ϕ d λ ) 2 , ( d t d λ ) 2 = ( 1 − 2 r ) − 2 ( d r d λ ) 2 + r 2 ( 1 − 2 r ) − 1 ( d ϕ d λ ) 2 , d t d λ = ( 1 − 2 r ) − 2 ( d r d λ ) 2 + r 2 ( 1 − 2 r ) − 1 ( d ϕ d λ ) 2 , d λ = [ ( 2 − 2 r ) − 2 ( d r d λ ) 2 + r 2 ( 1 − 2 r ) − 1 ( d ϕ d λ ) 2 ] − 1 / 2 d t , Δ λ = [ ( 2 − 2 r ) − 2 ( d r d λ ) 2 + r 2 ( 1 − 2 r ) − 1 ( d ϕ d λ ) 2 ] − 1 / 2 Δ t . \begin{align*} 0 &= -\left(1 - \frac{2}{r}\right) \left(\frac{dt}{d\lambda}\right)^2 + \left(1 - \frac{2}{r}\right)^{-1} \left(\frac{dr}{d\lambda}\right)^2 + r^2 \left(\frac{d\phi}{d\lambda}\right)^2, \\ \left(1 - \frac{2}{r}\right) \left(\frac{dt}{d\lambda}\right)^2 &= \left(1 - \frac{2}{r}\right)^{-1} \left(\frac{dr}{d\lambda}\right)^2 + r^2 \left(\frac{d\phi}{d\lambda}\right)^2, \\ \left(\frac{dt}{d\lambda}\right)^2 &= \left(1 - \frac{2}{r}\right)^{-2} \left(\frac{dr}{d\lambda}\right)^2 + r^2 \left(1 - \frac{2}{r}\right)^{-1} \left(\frac{d\phi}{d\lambda}\right)^2, \\ \frac{dt}{d\lambda} &= \sqrt{\left(1 - \frac{2}{r}\right)^{-2} \left(\frac{dr}{d\lambda}\right)^2 + r^2 \left(1 - \frac{2}{r}\right)^{-1} \left(\frac{d\phi}{d\lambda}\right)^2}, \\ d\lambda &= \left[\left(2 - \frac{2}{r}\right)^{-2} \left(\frac{dr}{d\lambda}\right)^2 + r^2 \left(1 - \frac{2}{r}\right)^{-1} \left(\frac{d\phi}{d\lambda}\right)^2\right]^{-1/2} dt, \\ \Delta \lambda &= \left[\left(2 - \frac{2}{r}\right)^{-2} \left(\frac{dr}{d\lambda}\right)^2 + r^2 \left(1 - \frac{2}{r}\right)^{-1} \left(\frac{d\phi}{d\lambda}\right)^2\right]^{-1/2} \Delta t. \end{align*} 0 ( 1 − r 2 ) ( d λ d t ) 2 ( d λ d t ) 2 d λ d t d λ Δ λ = − ( 1 − r 2 ) ( d λ d t ) 2 + ( 1 − r 2 ) − 1 ( d λ d r ) 2 + r 2 ( d λ d ϕ ) 2 , = ( 1 − r 2 ) − 1 ( d λ d r ) 2 + r 2 ( d λ d ϕ ) 2 , = ( 1 − r 2 ) − 2 ( d λ d r ) 2 + r 2 ( 1 − r 2 ) − 1 ( d λ d ϕ ) 2 , = ( 1 − r 2 ) − 2 ( d λ d r ) 2 + r 2 ( 1 − r 2 ) − 1 ( d λ d ϕ ) 2 , = [ ( 2 − r 2 ) − 2 ( d λ d r ) 2 + r 2 ( 1 − r 2 ) − 1 ( d λ d ϕ ) 2 ] − 1/2 d t , = [ ( 2 − r 2 ) − 2 ( d λ d r ) 2 + r 2 ( 1 − r 2 ) − 1 ( d λ d ϕ ) 2 ] − 1/2 Δ t . We can substitute h = Δ λ h = \Delta \lambda h = Δ λ
U n t = U n − 1 t − Δ λ 2 r n − 1 ( r n − 1 − 2 ) U n − 1 t U n − 1 r , U n r = U n − 1 r + Δ λ [ − r n − 1 − 2 ( r n − 1 ) 3 ( U n − 1 t ) 2 + 1 r n − 1 ( r n − 1 − 2 ) ( U n − 1 r ) 2 + ( r n − 1 − 2 ) ( U n − 1 ϕ ) 2 ] , U n ϕ = U n − 1 ϕ − Δ λ 2 r n − 1 U n − 1 r U n − 1 ϕ , t n = t n − 1 + Δ λ U n − 1 t , r n = r n − 1 + Δ λ U n − 1 r , ϕ n = ϕ n − 1 + Δ λ U n − 1 ϕ , Δ λ = [ ( 2 − 2 r n − 1 ) − 2 ( U n − 1 r ) 2 + r n − 1 2 ( 1 − 2 r n − 1 ) − 1 ( U n − 1 ϕ ) 2 ] − 1 / 2 Δ t . \begin{align*} U^t_n &= U^t_{n - 1} - \Delta \lambda \frac{2}{r_{n - 1}(r_{n - 1} - 2)} U^t_{n - 1} U^r_{n - 1}, \\ U^r_n &= U^r_{n - 1} + \Delta \lambda \left[- \frac{r_{n - 1} - 2}{(r_{n - 1})^3} \left(U^t_{n - 1}\right)^2 + \frac{1}{r_{n - 1}(r_{n - 1} - 2)} \left(U^r_{n - 1}\right)^2 + (r_{n - 1} - 2) \left(U^{\phi}_{n - 1}\right)^2\right], \\ U^{\phi}_n &= U^{\phi}_{n - 1} - \Delta \lambda \frac{2}{r_{n - 1}} U^r_{n - 1} U^{\phi}_{n - 1}, \\[4ex] t_n &= t_{n - 1} + \Delta \lambda U^t_{n - 1}, \\ r_n &= r_{n - 1} + \Delta \lambda U^r_{n - 1}, \\ \phi_n &= \phi_{n - 1} + \Delta \lambda U^{\phi}_{n - 1}, \\[4ex] \Delta \lambda &= \left[\left(2 - \frac{2}{r_{n - 1}}\right)^{-2} \left(U^r_{n - 1}\right)^2 + r_{n - 1}^2 \left(1 - \frac{2}{r_{n - 1}}\right)^{-1} \left(U^{\phi}_{n - 1}\right)^2\right]^{-1/2} \Delta t. \end{align*} U n t U n r U n ϕ t n r n ϕ n Δ λ = U n − 1 t − Δ λ r n − 1 ( r n − 1 − 2 ) 2 U n − 1 t U n − 1 r , = U n − 1 r + Δ λ [ − ( r n − 1 ) 3 r n − 1 − 2 ( U n − 1 t ) 2 + r n − 1 ( r n − 1 − 2 ) 1 ( U n − 1 r ) 2 + ( r n − 1 − 2 ) ( U n − 1 ϕ ) 2 ] , = U n − 1 ϕ − Δ λ r n − 1 2 U n − 1 r U n − 1 ϕ , = t n − 1 + Δ λ U n − 1 t , = r n − 1 + Δ λ U n − 1 r , = ϕ n − 1 + Δ λ U n − 1 ϕ , = [ ( 2 − r n − 1 2 ) − 2 ( U n − 1 r ) 2 + r n − 1 2 ( 1 − r n − 1 2 ) − 1 ( U n − 1 ϕ ) 2 ] − 1/2 Δ t .
