TP 3 Schéma des trapèzes implicites

import numpy as np
from scipy.integrate import solve_ivp  # solveur EDO de Python (initial value problem)
from matplotlib import pyplot as plt

def fun(t, y): 
    return 3* y - 4*np.exp(-t)

t0, tfinal = 0, 5
y0 = 1.0

Code du schéma des trapèzes implicites

pas=[];
erreur=[];
h=1;#pas de depart
nbiter=10; #on passe DIX fois dans la boucle
#le dernier pas de temps utilise est 1/2^9 
# la mise à jour du pas s'effectue en effet après le calcul de 
# la solution approchee

# boucle sur les pas de temps de plus en plus petits
for k in range(nbiter):
    #initialisation
    t=np.arange(t0,tfinal+h,h)
    N=len(t)
    yapp=np.zeros(N)
    y = y0
    yapp[0]=y0
    # code la methode des trapezes pour le pas h
    for n in range(N-1):
        y=yapp[n]
        #initialisation avec Euler
        z=y+h* fun(t[n],y)
        # iterations de point fixe
        for i in range(3):
            z = y + 0.5*h*(fun(t[n],y) + fun(t[n]+h,z))
        yapp[n+1]=z
    err=np.max(np.abs(yapp-np.exp(-t)));
    pas.append(h);
    erreur.append(err)
    h=h/2 #divise le pas par deux. c'est pour cela que nbiter=10 et pas 9

#trace la solutions exacte
xexact=np.arange(t0,tfinal,1e-2)
yexact=np.exp(-xexact);
plt.plot(xexact,yexact,'r');
plt.plot(t,yapp,'-k')#trace la solution approchee
plt.legend(['solution exacte', 'solution approchee trapèzes implicites'])
plt.xlabel('t')
plt.ylabel('y')
plt.show()

#trace l'erreur en fonction du pas en coordonnées log-log 
plt.loglog(pas, erreur,'-og')
plt.loglog(pas,np.power(pas,2),'--og')
plt.legend(['schéma trapèze implicite','pente 2'],loc='best')
plt.xlabel('log pas')
plt.ylabel('log erreur')
plt.title('ordre de convergence du schéma des trapèzes implicites')
plt.show()

#solveur basé sur schéma RK45 de Scipy
sol = solve_ivp(fun, [t0, tfinal], [1], 'RK45',dense_output=True)

#les temps choisis automatiquement par le solveur
sol.t

array([0.        , 0.10001999, 0.61293194, 1.20083981, 1.80458359,
       2.41190597, 3.02019605, 3.62771743, 4.1587852 , 4.67089781,
       5.        ])

plt.plot(sol.t, sol.y.T, '-ob')
# il faut transposer y pour qu'il ait le meme format que t 
plt.plot(t,yapp,'-k')
t_eval=np.arange(t0,tfinal,1e-2)
plt.plot(t_eval, np.exp(-t_eval),'-r')
plt.xlabel('t')
plt.ylabel('y')
plt.legend(['RK45','Trapèzes implicites', 'sol. exacte'], shadow=True)
plt.show()

La solution générale de l'équation différentielle est $ y(t) = C \, \exp(3t) + \exp(-t).$ Pour la solution exacte vérifiant la condition initiale, $y(t) = \exp(-t)$ i.e. $C=0.$ Les schémas numériques en s'écartant un peu de la solution exacte calculent des approximations avec $C \ll 1$ mais non rigoureusement nul. Au bout d'un certain temps le terme $C \, \exp(3t)$ devient dominant.

# tracé de solutions voisines de la solution qui s'en écartent progressivement.
plt.plot(t_eval, np.exp(-t_eval),'r')
for C in [-3e-7, -2e-7, -1e-7,-0.5e-7,0.5e-7,1e-7, 2e-7, 3e-7 ]:
    plt.plot(t_eval, C*np.exp(3*t_eval)+ np.exp(-t_eval))
plt.xlabel('t')
plt.ylabel('y')
plt.title('sol C exp(3t) + exp(-t) pour C entre -3e-7 et 3e-7')
plt.show()