Méthode d'Euler Richardson appliqué au problème de Kepler

Feuille python inspirée du travail de Damien Decout

Les méthodes classiques utilisées pour approcher les solutions des équations différentielles (Euler, Runge-Kutta ou Verlet) présentent l'inconvénient d'avoir un pas fixe. Or, dans de nombreuses situations, une grandeur physique peut varier rapidement un moment, puis plus lentement ensuite. Il est alors judicieux d'adapter le pas en fonction des besoins, c'est-à-dire de le diminuer lorsque les variations sont rapides et de l'augmenter dans l'autre cas, de façon à économiser du temps de calcul. C'est la philosophie des algorithmes à pas adaptatif dont la méthode décrite ici fait partie.

Pour en savoir plus, voir La Méthode d'Euler Richardson

Chargement des bibliothèques

import numpy as np
import matplotlib.pyplot as plt

Problème de gravité

L'équation différentielle du problème de Képler peut se mettre sous la forme : $$\dot{\mathbf{y}} = f(t,\mathbf{y}(t)),\quad 0\leq t \leq T \quad\text{avec}\quad \mathbf{y}(0)=\mathbf{y}_{0}$$ avec: $$\mathbf{y}=\left(\begin{array}{c}x\\y\\v_{x}\\v_{y}\end{array}\right) \quad\text{et}\quad f(\mathbf{y})=\left(\begin{array}{c}v_{x}\\v_{y}\\[1mm]-\dfrac{4\pi^{2}x}{(x^{2}+y^{2})^{3/2}}\\[3mm] -\dfrac{4\pi^{2}y}{(x^{2}+y^{2})^{3/2}} \end{array}\right)$$

def f(t,u):
    G = 4*(np.pi)**2
    x = u[0]
    y = u[1]
    ax = -G*x/(x**2+y**2)**1.5
    ay = -G*y/(x**2+y**2)**1.5
    return np.array([u[2],u[3],ax,ay])

Algorithme d'Euler Richardson :

• Initialisation du pas $h$, de la durée $T$ et du seuil de précision $\varepsilon_{\text{seuil}}$.
• Initialisation des conditions initiales : $t=0$ et $\mathbf{y}=\mathbf{y}(0)$.
• Tant que $t\leq T$ faire :
  • Calcul de $\mathbf{k}=f(t,\mathbf{y})$.
  • Calcul de $\mathbf{k'}=f(t+h/2,y+\mathbf{k}h/2)$.
  • Calcul de $\varepsilon=\frac{h}{2}\left|\mathbf{k}'-\mathbf{k}\right|$ et $\alpha=\varepsilon/\varepsilon_{\rm seuil}$.
  • Si $\alpha > 1$ faire $h=0,9\,h/\sqrt{\alpha}$.
  • Si $\alpha <1$ faire $\mathbf{y}=\mathbf{y}+\mathbf{k'}h$ ; $h=0,9\,h/\sqrt{\alpha}$ ; $t=t+h$ ; enregistrement des données.

def eulerRichardson(f,t0,y0,h0,tmax,precision):
    X = []
    Y = []
    T = []
    h = h0
    seuil = precision
    t = t0
    y = np.array(y0)
    while t<=tmax:
        epsilon=0
        k = f(t,y)
        kp = f(t+h/2,y+k*h/2)
        for i in range(0,len(y0)):
            epsilon += (h/2)*abs(kp[i]-k[i])    
        alpha = epsilon/seuil
        if alpha > 1:
            h = 0.9*h/np.sqrt(alpha)
        if alpha < 1:
            y = y+kp*h
            h = 0.9*h/np.sqrt(alpha)
            t = t+h
        X.append(y[0])
        Y.append(y[1])
        T.append(t)
    return T,X,Y

On adopte les conditions initiales : $x=0,5\quad y=0\quad v_x=0\quad v_y=12$

On choisit un pas initial de $h=0,04$ (environ 15 jours) et un seuil de 0,018.

ER = eulerRichardson(f,0,[0.5,0,0,12],0.04,4.768,0.018)

plt.plot(ER[1],ER[2],marker='.',label='Euler-Richardson')
plt.xlabel(r'$x(u.a)$')
plt.ylabel(r'$y(u.a)$')
plt.scatter(0,0,200,c='y')
plt.legend();

Comparaison avec RK4

Algorithme de Runge-Kutta :

• Initialisation du pas $h$, de la durée $T$.
• Initialisation des conditions initiales : $t=0$ et $\mathbf{y}=\mathbf{y}(0)$.
• Définition de la fonction $f(t,y).$
• Tant que $t\leq T$ faire :
  • Calcul de $k_1=f(t,y)$.
  • Calcul de $k_2=f(t+h/2,y+hk_1/2)$.
  • Calcul de $k_3= f(t+h/2,y+hk_2/2)$.
  • Calcul de $k_4=f(t+h,y+hk_3)$.
  • $y=y+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$ ; $t=t+h$.
  • Enregistrement des données

def RK4_2D(f,t0,y0,h,tmax):
    t = t0
    y = np.array(y0)
    X = [y[0]]
    Y = [y[1]]
    T = [0]
    while t<=tmax:
        k1 = f(t,y)
        k2 = f(t+h/2,y+h*k1/2)
        k3 = f(t+h/2,y+h*k2/2)
        k4 = f(t+h,y+h*k3)
        y = y + (h/6)*(k1+2*k2+2*k3+k4)
        t = t + h
        X.append(y[0])
        Y.append(y[1])
        T.append(t)
    return T,X,Y

RK4 = RK4_2D(f,0,[0.5,0,0,12],0.04,4.768)

plt.plot(RK4[1],RK4[2],marker='.',label='RK4')
plt.xlabel(r'$x(u.a)$')
plt.ylabel(r'$y(u.a)$')
plt.scatter(0,0,200,c='y')
plt.legend();

Comparaison

plt.plot([0,0], [-2,2],color='.8',lw=1)
plt.plot([-5.5,1], [0,0],color='.8',lw=1)
plt.plot(RK4[1],RK4[2],label='Runge-Kutta 4',c='orange')
plt.plot(ER[1],ER[2],label='Euler-Richardson',c='b')
plt.xlabel('$x$ (u.a)')
plt.ylabel('$y$ (u.a)')
plt.scatter(0,0,200,c='y')
plt.scatter(RK4[1][-1],RK4[2][-1],25,c='orange')
plt.scatter(ER[1][-1],ER[2][-1],25,c='b')
plt.legend(loc='upper right', bbox_to_anchor=(1,1));
plt.title('Comparaison pour un pas initial $h=0,04\simeq$ 15 jours',loc='center')
plt.savefig("test.svg", format="svg")