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 \((\frac{du}{d\varphi})^2=R_su^3-u^2+\frac{1}{b^2}\), l’équation différentielle peut se traiter en évitant la singularité due à \(\pm\) la racine.
La dérivation de \((\frac{du}{d\varphi})^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}(u(\varphi)+\frac{\delta \varphi}{2}\frac{du}{d\varphi})\), \(k_3=\frac{d^2u}{d\varphi^2}(u(\varphi)+\frac{\delta \varphi}{2}\frac{du}{d\varphi}+(\frac{\delta \varphi}{2})^2\ k_1)\) et \(k_4=\frac{d^2u}{d\varphi^2}(u(\varphi)+\delta \varphi\ \frac{du}{d\varphi}+\frac{(\delta \varphi)^2}{2}\ k_2)\)
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}(u(\varphi)+\frac{\delta \varphi}{2}l_1)\), \(l_3=\frac{du}{d\varphi}(u(\varphi)+\frac{\delta \varphi}{2}l_2)\) 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\) (\(={1\over u}\)).
Les conditions initiales de \(\frac{du}{d\varphi}\) et de \(u\) se déterminent comme suit, à partir des coordonnées d’émission (\(r_{em}\) , \(\varphi_{em}\)) du photon (\(u_{em}={1\over r_{em}}\)), et deux cas possibles pour \(\frac{du}{d\varphi}\) :
1) photon entrant (\(\frac{du}{dt}>0)\) et \(\frac{d\varphi}{dt}>0\) ou photon sortant (\(\frac{du}{dt}<0)\) 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 (\(\frac{du}{dt}<0)\) et \(\frac{d\varphi}{dt}>0\) ou photon entrant (\(\frac{du}{dt}>0)\) 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 (\(\frac{du}{dt}>0\)) et les deux cas possibles deviennent :
1) \(\frac{du}{d\varphi}=+{1\over b}\),
2) \(\frac{du}{d\varphi}=-{1\over b}\).