In [20]:
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
plt.rc('figure',figsize=(14,9))

Résolution d'équations différentielles

Premier cas : une équation avec un phénomène d'explosion

$$ \dot x(t) = [x(t)]^2 \quad \hbox{pour }t>0,\; x(0)=1\; .$$

La solution est :

$$ x(t) = \frac{1}{1-t} \quad \hbox{pour }0 \leq t < 1\; .$$

Mise en place du schéma d'Euler :

In [21]:
h=0.067
t=np.arange(0.0,1.32,h)
n=len(t)
x=np.zeros(n)
x[0]=1.0
for i in range(n-1):
    x[i+1]=x[i]+h*x[i]*x[i]
plt.plot(t,x)
Out[21]:
[<matplotlib.lines.Line2D at 0x109bc7fd0>]

Remarque : Le schéma d'Euler prédit mal le temps d'explosion!

Mise en place de la méthode de Runge-Kutta

In [35]:
def F2(y):
    z=y*y
    return z
In [36]:
def RKF(y,FK,h):
    k1 = h * FK(y)
    k2 = h * FK(y + k1/2.0)
    k3 = h * FK(y + k2/2.0)
    k4 = h * FK(y + k3)
    
    return y + ( k1 + 2.0 * k2 + 2.0 * k3 + k4 ) / 6.0
In [37]:
s=[0.0]
y=[1.0]
for i in range(n-1):
    if y[i]<100.0:
        s.append(s[i]+h)
        y.append(RKF(y[i],F2,h))
    else: break

plt.plot(s,y)
Out[37]:
[<matplotlib.lines.Line2D at 0x10a6c5210>]

Remarque : avec Runge-Kutta, la prédiction est bien meilleure!

Équations de Képler en dimension 2 :

In [38]:
def FK(y):
    N=(y[0]*y[0]+y[1]*y[1])**1.5
    g=np.zeros(4)
    g[0]=y[2]
    g[1]=y[3]
    g[2]=-3.0*y[0]/N
    g[3]=-3.0*y[1]/N 

    return g

Vérification :

In [39]:
y=np.array([0,1.0,0,0])

g=FK(y)

print g

[ 0.  0. -0. -3.]

Le solveur avec la méthode d'Euler :

In [40]:
def solLV(h,T=80.0,x0=4.0,y0=4.2,u0=-0.40,v0=0.40):
    t=np.arange(0.0,T,h)
    n=len(t)
    x=np.zeros((n,4))
    x[0,0],x[0,1],x[0,2],x[0,3]=x0,y0,u0,v0
    for i in range(n-1):
        x[i+1,:]=x[i,:]+h*FK(x[i,:])
          
    return (t,x)
In [41]:
t,x=solLV(0.01)
plt.plot(x[:,0],x[:,1])
plt.xlim([-10.0,10.0])
plt.ylim([-10.0,10.0])
Out[41]:
(-10.0, 10.0)

Remarque : la trajectoire n'est pas vraiment périodique (comme on s'y attendrait pourtant...).

In [43]:
from IPython.html.widgets import interactive
from IPython.display import display
In [44]:
def visu(k):
    t,x=solLV(0.01)
    plt.plot(x[:k,0],x[:k,1])
    plt.xlim([-8.0,8.0])
    plt.ylim([-8.0,8.0])
In [45]:
wdgt=interactive(visu,k=(2000,4000,200))
display(wdgt)

Résolution avec Runge-Kutta :

In [46]:
def solRK(h,T=40.0,x0=4.0,y0=4.2,u0=-0.40,v0=0.40):
    t=np.arange(0.0,T,h)
    n=len(t)
    x=np.zeros((n,4))
    x[0,0],x[0,1],x[0,2],x[0,3]=x0,y0,u0,v0
    for i in range(n-1):
        x[i+1,:]=RKF(x[i,:],FK,h)
          
    return (t,x)
In [47]:
t,x=solRK(0.01)
plt.plot(x[:,0],x[:,1])
plt.xlim([-8.0,8.0])
plt.ylim([-8.0,8.0])
Out[47]:
(-8.0, 8.0)

Remarque : du point de vue de la périodicité, Runge-Kutta fait bien mieux...

In [48]:
def visuRK(k):
    t,x=solRK(0.01)
    plt.plot(x[:k,0],x[:k,1])
    plt.xlim([-8.0,8.0])
    plt.ylim([-8.0,8.0])
In [49]:
wdgt=interactive(visuRK,k=(2000,4000,200))
display(wdgt)
In [ ]:
In [ ]: