Avant-propos : l’équation \(\frac{du}{d\phi}=\pm\sqrt{R_su^3-u^2+\frac{1}{b^2}}\) permet de calculer \(\frac{d\phi}{du}\), et \(\phi(u)\) peut alors se calculer analytiquement par une intégrale elliptique de \(\frac{d\phi}{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\phi})^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\phi})^2\) donne en effet :
\(2\frac{du}{d\phi}\frac{d^2u}{d\phi^2}=3R_su^2\frac{du}{d\phi}-2u\frac{du}{d\phi}\) soit \(\frac{d^2u}{d\phi^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\phi}{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\phi^2}(u(\phi))\), \(k_2=\frac{d^2u}{d\phi^2}(u(\phi)+\frac{\delta \phi}{2}\frac{du}{d\phi})\), \(k_3=\frac{d^2u}{d\phi^2}(u(\phi)+\frac{\delta \phi}{2}\frac{du}{d\phi}+(\frac{\delta \phi}{2})^2\ k_1)\) et \(k_4=\frac{d^2u}{d\phi^2}(u(\phi)+\delta \phi\ \frac{du}{d\phi}+\frac{(\delta \phi)^2}{2}\ k_2)\)
la valeur de \(\frac{du}{d\phi}\) au pas suivant est :
\(\frac{du}{d\phi}(\phi+\delta \phi)=\frac{du}{d\phi}(\phi)+\frac{\delta \phi}{6}(k_1+2k_2+2k_3+k_4)\), et
2) avec \(l_1=\frac{du}{d\phi}(u(\phi))\), \(l_2=\frac{du}{d\phi}(u(\phi)+\frac{\delta \phi}{2}l_1)\), \(l_3=\frac{du}{d\phi}(u(\phi)+\frac{\delta \phi}{2}l_2)\) et \(l_4=\frac{du}{d\phi}(u(\phi)+\delta \phi\ l_3)\)
la valeur de \(u\) au pas suivant est :
\(u(\phi+\delta \phi)=u(\phi)+\frac{\delta \phi}{6}(l_1+2l_2+2l_3+l_4)\)
ce qui s’écrit avec \(l_2=\frac{du}{d\phi}(\phi)+\frac{\delta \phi}{2}k_1\), \(l_3=\frac{du}{d\phi}(\phi)+\frac{\delta \phi}{2}k_2\) et \(l_4=\frac{du}{d\phi}(\phi)+\delta \phi\ k_3\)
et après développement :
\(u(\phi+\delta \phi)=u(\phi)+\delta \phi\frac{du}{d\phi}(\phi)+\frac{(\delta \phi)^2}{6}(k_1+k_2+k_3)\).
Le pas de calcul \(\delta \phi\) 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\phi}\) et de \(u\) se déterminent comme suit, à partir des coordonnées d’émission (\(r_{em}\) , \(\phi_{em}\)) du photon (\(u_{em}={1\over r_{em}}\)), et deux cas possibles pour \(\frac{du}{d\phi}\) :
1) photon entrant (\(\frac{du}{dt}>0)\) et \(\frac{d\phi}{dt}>0\) ou photon sortant (\(\frac{du}{dt}<0)\) et \(\frac{d\phi}{dt}<0\) : \(\frac{du}{d\phi}=+\sqrt{R_su_{em}^3-u_{em}^2+\frac{1}{b^2}}\)
2) photon sortant (\(\frac{du}{dt}<0)\) et \(\frac{d\phi}{dt}>0\) ou photon entrant (\(\frac{du}{dt}>0)\) et \(\frac{d\phi}{dt}<0\) : \(\frac{du}{d\phi}=-\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\phi}=+{1\over b}\)
2) \(\frac{du}{d\phi}=-{1\over b}\)