Avant-propos : l’équation de la trajectoire des photons avec la métrique de Schwarzschild \(\frac{du}{d\varphi}=\pm\sqrt{R_su^3-u^2+\frac{1}{b^2}}\) permet de calculer \(\frac{d\varphi}{du}\), et \(\varphi(u)\) peut alors se calculer analytiquement par une intégrale elliptique de \(\frac{d\varphi}{du}\) dans le cas où il existe un périastre, ce qui n’est pas toujours le cas suivant la valeur de \(b\). De plus, le cas d’un photon émis à une coordonnée radiale très proche de \(R_s\) ne peut être traité par cette méthode.
En écrivant \(\left(\frac{du}{d\varphi}\right)^2=R_su^3-u^2+\frac{1}{b^2}\), l’équation différentielle peut se traiter en évitant la gestion du \(\pm\) de la racine.
La dérivation de \(\left(\frac{du}{d\varphi}\right)^2\) donne en effet :
\(2\frac{du}{d\varphi}\frac{d^2u}{d\varphi^2}=3R_su^2\frac{du}{d\varphi}-2u\frac{du}{d\varphi}\) soit \(\frac{d^2u}{d\varphi^2}={3\over 2}R_su^2-u\).
Cette équation différentielle peut alors se résoudre numériquement par double intégration avec les conditions initiales et lever les limites de l’intégrale elliptique de \(\frac{d\varphi}{du}\) vues précédemment.
Double intégration – Méthode Runge-Kutta d’ordre 4
L’équation avec dérivée seconde vue ci-dessus peut s’intégrer numériquement comme suit par une méthode de Runge-Kutta d’ordre 4 :
1) avec \(k_1=\frac{d^2u}{d\varphi^2}(u(\varphi))\), \(k_2=\frac{d^2u}{d\varphi ^2}\left(u\left(\varphi\right)+\frac{\delta \varphi}{2}\frac{du}{d\varphi}\right)\), \(k_3=\frac{d^2u}{d\varphi^2}\left(u\left(\varphi\right)+\frac{\delta \varphi}{2}\frac{du}{d\varphi}+\left(\frac{\delta \varphi}{2}\right)^2\ k_1\right)\) et \(k_4=\frac{d^2u}{d\varphi^2}\left(u\left(\varphi\right)+\delta \varphi\ \frac{du}{d\varphi}+\frac{\left(\delta \varphi\right)^2}{2}\ k_2\right)\)
la valeur de \(\frac{du}{d\varphi}\) au pas suivant est :
\(\frac{du}{d\varphi}(\varphi+\delta \varphi)=\frac{du}{d\varphi}(\varphi)+\frac{\delta \varphi}{6}(k_1+2k_2+2k_3+k_4)\), et
2) avec \(l_1=\frac{du}{d\varphi}(u(\varphi))\), \(l_2=\frac{du}{d\varphi}\left(u\left(\varphi\right)+\frac{\delta \varphi}{2}l_1\right)\), \(l_3=\frac{du}{d\varphi}\left(u\left(\varphi\right)+\frac{\delta \varphi}{2}l_2\right)\) et \(l_4=\frac{du}{d\varphi}(u(\varphi)+\delta \varphi \ l_3)\)
la valeur de \(u\) au pas suivant est :
\(u(\varphi+\delta \varphi)=u(\varphi)+\frac{\delta \varphi}{6}(l_1+2l_2+2l_3+l_4)\)
ce qui s’écrit avec \(l_2=\frac{du}{d\varphi}(\varphi)+\frac{\delta \varphi}{2}k_1\), \(l_3=\frac{du}{d\varphi}(\varphi)+\frac{\delta \varphi}{2}k_2\) et \(l_4=\frac{du}{d\varphi}(\varphi)+\delta \varphi\ k_3\)
et après développement :
\(u(\varphi+\delta \varphi)=u(\varphi)+\delta \varphi \frac{du}{d\varphi}(\varphi)+\frac{(\delta \varphi)^2}{6}(k_1+k_2+k_3)\).
Le pas de calcul \(\delta \varphi \) peut être choisi constant ou bien adaptatif suivant la précision souhaitée sur \(r\left(={1\over u}\right)\).
Les conditions initiales de \(\frac{du}{d\varphi}\) et de \(u\) se déterminent comme suit, à partir des coordonnées d’émission \(\left(r_{em},\varphi_{em}\right)\) du photon avec \(u_{em}={1\over r_{em}}\), et deux cas possibles pour \(\frac{du}{d\varphi}\) :
1) photon entrant \(\left(\frac{du}{dt}>0\right)\) et \(\frac{d\varphi}{dt}>0\) ou photon sortant \(\left(\frac{du}{dt}<0\right)\) et \(\frac{d\varphi}{dt}<0\) : \(\frac{du}{d\varphi}=+\sqrt{R_su_{em}^3-u_{em}^2+\frac{1}{b^2}}\)
2) photon sortant \(\left(\frac{du}{dt}<0\right)\) et \(\frac{d\varphi}{dt}>0\) ou photon entrant \(\left(\frac{du}{dt}>0\right)\) et \(\frac{d\varphi}{dt}<0\) : \(\frac{du}{d\varphi}=-\sqrt{R_su_{em}^3-u_{em}^2+\frac{1}{b^2}}\)
Si le photon arrive depuis l’infini, \(u_{em}=0\), le photon est entrant \(\left(\frac{du}{dt}>0\right)\) et les deux cas possibles deviennent :
1) \(\frac{du}{d\varphi}=+{1\over b}\),
2) \(\frac{du}{d\varphi}=-{1\over b}\).