TP 2 - Implémentation des schémas “classiques” et étude de la convergence

Auteur·rice

Gloria FACCANONI

Date de publication

22 mars 2024

Code

import sys #only needed to determine Python version number
print('Python version ' + sys.version)

Python version 3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0]

1 Exercice (resolu) : implémentation des schémas d’Euler explicite et implicite

Considérons le problème de Cauchy

trouver la fonction \(y \colon I\subset \mathbb{R} \to \mathbb{R}\) définie sur l’intervalle \(I=[0,1]\) telle que
\( \begin{cases} y'(t) = -4y(t)+t^2, &\forall t \in I=[0,1],\\ y(0) = 1 \end{cases} \)

Calculer la solution exacte en utilisant le module sympy.
Calculer la solution approchée obtenue avec la méthode d’Euler explicite avec \(h=1/N\) et \(N=8\) (pour bien visualiser les erreurs);
Même exercice pour la méthode d’Euler implicite.

Correction 1
Calculons la solution exacte en utilisant le module sympy:

Code

%reset -f
%matplotlib inline
%autosave 300

import sympy as sym
sym.init_printing()

t0 = 0
y0 = 1

t = sym.Symbol('t')
y = sym.Function('y')
edo = sym.Eq( sym.diff(y(t),t) , -4*y(t)+t**2 )
display(edo)

# solgen = sym.dsolve(edo)
# display(solgen)
# cond_init = sym.Eq( y0, solgen.rhs.subs(t,t0))
# display(cond_init)
# consts = sym.solve( cond_init , dict=True)[0]
# display(consts)
# solpar = solgen.subs(consts).simplify()

solpar = sym.dsolve(edo, ics={y(t0): y0})
display(solpar)

Autosaving every 300 seconds

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = t^{2} - 4 y{\left(t \right)}\)

\(\displaystyle y{\left(t \right)} = \frac{t^{2}}{4} - \frac{t}{8} + \frac{1}{32} + \frac{31 e^{- 4 t}}{32}\)

On définit la solution exacte à utiliser pour estimer les erreurs:

Code

sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')

Correction 2 et 3
On commence par importer le module matplotlib et la fonction fsolve du module scipy.optimize pour résoudre l’équation implicite présente dans le schéma implicite:

Code

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 12}) # plt.rcParams.update(rcParamsDefault)
from scipy.optimize import fsolve # pour les schémas implicites

On initialise le problème de Cauchy

Code

# variables globales
t0     = 0
tfinal = 1
y0     = 1

On définit l’équation différentielle : phi est une lambda function qui contient la fonction mathématique \(\varphi(t, y)=-4y+t^2\) dépendant des variables \(t\) et \(y\).

Code

phi = lambda t,y : -4*y+t**2

On écrit les schémas numériques : les valeurs \([u_0,u_1,\dots,u_{N}]\) pour chaque méthode sont contenues dans le vecteur uu.

Schéma d’Euler progressif :

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_n,u_n)& n=0,1,2,\dots N-1 \end{cases} \)

Code

def EE(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        uu[n+1] = uu[n] + h * phi(tt[n], uu[n])
    return uu

Schéma d’Euler régressif :

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1})& n=0,1,2,\dots N-1 \end{cases} \)

Attention : \(u_{n+1}\) est solution de l’équation \(x=u_n+h\varphi(t_{n+1},x)\), c’est-à-dire un zéro de la fonction (en générale non linéaire)

\( x\mapsto -x+u_n+h\varphi(t_{n+1},x) \)

Code

def EI(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        temp = fsolve(lambda x: -x + uu[n] + h * phi(tt[n + 1], x), uu[n])
        uu[n+1] = temp[0]
    return uu

On introduit la discrétisation: les nœuds d’intégration \([t_0,t_1,\dots,t_{N}]\) sont contenus dans le vecteur tt.
On a \(N+1\) points espacés de \(h=\frac{t_N-t_0}{N}\).

On calcule les solutions exacte et approchées en ces points:

Code

N  = 8 
tt = np.linspace(t0,tfinal,N+1)

# yy    = np.array([sol_exacte(t) for t in tt], dtype=object)
yy    = sol_exacte(tt)  # si vectorisée

uu_ep = EE(phi, tt, y0)
uu_er = EI(phi, tt, y0)

On compare les graphes des solutions exacte et approchées et on affiche le maximum de l’erreur: \( \max_n \left| y(t_n)-u_n \right| \)

Code

plt.figure(figsize=(18, 7))

plt.subplot(1, 2, 1)
plt.plot(tt, yy,    'b-')
plt.plot(tt, uu_ep, 'r-D')
err_ep = np.linalg.norm(uu_ep-yy,np.inf)
plt.title(f'Euler explicite \nmax(|erreur|)= {err_ep:1.10f}')
plt.grid()

plt.subplot(1, 2, 2)
plt.plot(tt, yy,    'b-')
plt.plot(tt, uu_er, 'r-D')
err_er = np.linalg.norm(uu_er-yy,np.inf)
plt.title(f'Euler implicite\nmax(|erreur|)= {err_er:1.10f}')
plt.grid();

2 Exercice (resolu) : étude empirique de la convergence des schémas d’Euler explicite et implicite

Considérons le même problème de Cauchy.

On se propose d’estimer empiriquement l’ordre de convergence de la méthode d’Euler explicite.
- On calcule d’abord la solution approchée avec différentes valeurs de \(h_k=1/N_k\)
- Pour chaque valeur de \(h_k\), on calcule le maximum de la valeur absolue de l’erreur et on la sauvegarde dans le vecteur err_ep de sort que err_ep[k] contient \(e_k=\max_{i=0,\dots,N_k}|y(t_i)-u_{i}|\).
- Pour estimer l’ordre de convergence on affiche les points (h[k],err_ep[k]) en echèlle logarithmique. On trouve ainsi une droite qui relie l’erreur au pas \(k\) à l’erreur au pas \(k+1\). Pour estimer la pente globale de cette droite (par des moindres carrés) on utilisers la fonction polyfit avec une régression linéaire.
Même exercice pour estimer l’ordre de convergence de la méthode d’Euler implicite.

Correction

Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de \(h_k=1/N_k\) et on sauvegarde les valeurs de \(h_k\) dans le vecteur H.

Pour chaque valeur de \(h_k\), on calcule le maximum de la valeur absolue de l’erreur et on sauvegarde toutes ces erreurs dans le vecteur err_schema de sort que err_schema[k] contient \(e_k=\max_{i=0,\dots,N_k}|y(t_i)-u_{i}|\).

Code

H = []
err_ep = []
err_er = []

N = 10
for k in range(7):
    N += 20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)    
    uu_ep = EE(phi, tt, y0)
    uu_er = EI(phi, tt, y0)
    H.append(h)
    err_ep.append(np.linalg.norm(uu_ep-yy,np.inf))
    err_er.append(np.linalg.norm(uu_er-yy,np.inf))

Pour afficher l’ordre de convergence on utilise une échelle logarithmique, i.e. on représente \(\ln(h)\) sur l’axe des abscisses et \(\ln(\text{err})\) sur l’axe des ordonnées.
En effet, si \(\text{err}=Ch^p\) alors \(\ln(\text{err})=\ln(C)+p\ln(h)\).
En échelle logarithmique, \(p\) représente donc la pente de la ligne droite \(\ln(\text{err})\).

Code

plt.figure(figsize=(10,7))
plt.plot(np.log(H), np.log(err_ep), 'r-o', label='Euler Explicite')
plt.plot(np.log(H), np.log(err_er), 'g-+', label='Euler Implicite')
# plt.loglog(H, err_ep, 'r-o', label='Euler Explicite')
# plt.loglog(H, err_er, 'g-+', label='Euler Implicite')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True)

Pour estimer l’ordre de convergence on doit estimer la pente de la droite qui relie l’erreur au pas \(k\) à l’erreur au pas \(k+1\) en echelle logarithmique. Pour estimer la pente globale de cette droite (par des moindres carrés) on peut utiliser la fonction polyfit (du module numpy que nous avons déjà importé avec matplotlib.pylab).

Code

order_ep = np.polyfit(np.log(H), np.log(err_ep), 1)[0]
order_er = np.polyfit(np.log(H), np.log(err_er), 1)[0]

print (f'Euler progressif : ordre = {order_ep:1.2f}')
print (f'Euler regressif  : ordre = {order_er:1.2f}')

Euler progressif : ordre = 1.03
Euler regressif  : ordre = 0.98

3 Exercice : implémentation et étude empirique de la convergence des méthodes d’Euler modifié, RK1_m, de Crank-Nicolson et de Heun

Ajouter l’implémentation et l’étude empirique de la convergence des méthodes d’Euler modifié, RK1_M, de Crank-Nicolson et de Heun .

Pour cela modifier directement le code ci-dessus : les endroits où il faut écrire du code sont indiqués par les commentaire commençant par # !!!!!!.

Correction

On écrit les quatre schémas:

Schéma d’Euler modifié:

\( \begin{cases} u_0=y_0,\\ \tilde u = u_n+\frac{h}{2}\varphi(t_n,u_n),\\ u_{n+1}=u_n+h\varphi\left(t_n+\frac{h}{2},\tilde u\right)& n=0,1,2,\dots N-1 \end{cases} \)

qu’on peut réécrire sous la forme

\( \begin{cases} u_0=y_0,\\ k_1 = \varphi(t_n,u_n),\\ u_{n+1}=u_n+h\varphi\left(t_n+\frac{h}{2}, u_n+h\frac{k_1}{2}\right)& n=0,1,2,\dots N-1 \end{cases} \)

Code

def EM(phi,tt, y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        uu[n+1] = uu[n]+h*phi( tt[n]+h/2 , uu[n]+h*k1/2 ) 
    return uu

Schéma RK1_M:

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi\left(\dfrac{t_n+t_{n+1}}{2},\dfrac{u_n+u_{n+1}}{2}\right)& n=0,1,2,\dots N-1 \end{cases} \)

avec \(u_{n+1}\) zéro de la fonction \(x\mapsto -x+u_n+h\varphi(t_n+\frac{h}{2},\frac{u_n+x}{2})\).

Code

def RK1_M(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+h*phi( (tt[n]+tt[n+1])/2,(uu[n]+x)/2 ), uu[n])
        uu[n+1] = temp[0]
    return uu

Schéma de Crank-Nicolson :

\( \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(\varphi(t_n,u_n)+\varphi(t_{n+1},u_{n+1})\Bigr)& n=0,1,2,\dots N-1 \end{cases} \)

avec \(u_{n+1}\) zéro de la fonction \(x\mapsto -x+u_n+\frac{h}{2}(\varphi(t_n,u_n)+\varphi(t_{n+1},x))\).

Code

def CN(phi,tt, y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        temp = fsolve(lambda x: -x+uu[n]+0.5*h*( phi(tt[n+1],x)+phi(tt[n],uu[n]) ), uu[n])
        uu[n+1] = temp[0]
    return uu

Schéma de Heun:

\( \begin{cases} u_0=y_0,\\ \tilde u = u_n+h\varphi(t_n,u_n)\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(\varphi(t_n,u_n)+\varphi(t_{n+1},\tilde u)\Bigr)& n=0,1,2,\dots N-1 \end{cases} \)

qu’on peut réécrire sous la forme

\( \begin{cases} u_0=y_0,\\ k_1 = \varphi(t_n,u_n),\\ k_2 = \varphi(t_{n+1},u_n+hk_1),\\ u_{n+1}=u_n+h\frac{k_1+k_2}{2}& n=0,1,2,\dots N-1 \end{cases} \)

Code

def heun(phi,tt, y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        k2 = phi( tt[n+1], uu[n] + h*k1  )
        uu[n+1] = uu[n] + h*(k1+k2)/2 
    return uu

Code

N = 8 
tt = np.linspace(t0,tfinal,N+1)

# yy      = np.array([sol_exacte(t) for t in tt], dtype=object)
yy      = sol_exacte(tt)  # si vectorisée
uu_ep   = EE(phi, tt, y0)
uu_er   = EI(phi, tt, y0)
uu_em   = EM(phi, tt, y0)
uu_rk1m = RK1_M(phi, tt, y0)
uu_cn   = CN(phi, tt, y0)
uu_heun = heun(phi, tt, y0)

plt.figure(figsize=(24, 14))

plt.subplot(2, 3, 1)
plt.plot(tt, yy, 'b-')
plt.plot(tt, uu_ep, 'r-D')
# err_ep = max( [abs(uu_ep[i] - yy[i]) for n in range(N)] )
err_ep = np.linalg.norm( uu_ep-yy, np.inf)
plt.title(f'Euler explicite\nmax(|erreur|)= {err_ep:1.10f}')
plt.grid()

plt.subplot(2, 3, 2)
plt.plot(tt, yy, 'b-', tt, uu_er, 'r-D')
# err_er = max( [abs(uu_er[i] - yy[i]) for n in range(N)] )
err_er = np.linalg.norm(uu_er-yy, np.inf)
plt.title(f'Euler implicite\nmax(|erreur|)={err_er:1.2f}')
plt.grid();


plt.subplot(2, 3, 3)
plt.plot(tt, yy, 'b-', tt, uu_em, 'r-D')
# err_em = max( [abs(uu_em[i] - yy[i]) for n in range(N)] )
err_em = np.linalg.norm(uu_em-yy, np.inf)
plt.title(f'Euler modifié\nmax(|erreur|)={err_em:1.5f}')
plt.grid()

plt.subplot(2, 3, 4)
plt.plot(tt, yy, 'b-', tt, uu_rk1m, 'r-D')
# err_rk1m = max( [abs(uu_rk1m[i] - yy[i]) for n in range(N)] )
err_rk1m = np.linalg.norm(uu_rk1m-yy, np.inf)
plt.title(f'RK1_m\nmax(|erreur|)={err_rk1m:1.5f}')
plt.grid()

plt.subplot(2, 3, 5)
plt.plot(tt, yy, 'b-', tt, uu_cn, 'r-D')
# err_cn = max( [abs(uu_cn[i] - yy[i]) for n in range(N)] )
err_cn = np.linalg.norm(uu_cn-yy,np.inf)
plt.title(f'Crank Nicolson\nmax(|erreur|)={err_cn:1.5f}')
plt.grid()

plt.subplot(2, 3, 6)
plt.plot(tt, yy, 'b-', tt, uu_heun, 'r-D')
# err_heun = max( [abs(uu_heun[i] - yy[i]) for n in range(N)] )
err_heun = np.linalg.norm(uu_heun-yy, np.inf)
plt.title(f'Heun\nmax(|erreur|)={err_heun:1.5f}')
plt.grid();

Code

H = []
err_ee = []
err_ei = []
err_em = []
err_rk1m = []
err_cn = []
err_heun = []

for k in range(7):
    N = 2**(k + 3)
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    # yy = np.array([sol_exacte(t) for t in tt], dtype=object)
    yy = sol_exacte(tt)  # si vectorisée
    uu_ee   = EE(phi, tt, y0)
    uu_ei   = EI(phi, tt, y0)
    uu_em   = EM(phi, tt, y0)
    uu_rk1m = RK1_M(phi, tt, y0)
    uu_cn   = CN(phi, tt, y0)
    uu_heun = heun(phi, tt, y0)
    H.append(h)
    err_ee.append(np.linalg.norm(uu_ee-yy,np.inf))
    err_ei.append(np.linalg.norm(uu_ei-yy,np.inf))
    err_em.append(np.linalg.norm(uu_em-yy,np.inf))
    err_rk1m.append(np.linalg.norm(uu_rk1m-yy,np.inf))
    err_cn.append(np.linalg.norm(uu_cn-yy,np.inf))
    err_heun.append(np.linalg.norm(uu_heun-yy,np.inf))


# Plot des droites des erreurs vs taille de grille
plt.figure(figsize=(10,7))
plt.loglog(H, err_ee,   'r-o', label='Euler Explicite')
plt.loglog(H, err_ei,   'g-+', label='Euler Implicite')
plt.loglog(H, err_em,   'm-x', label='Euler Modife')
plt.loglog(H, err_rk1m, 'b-^', label='RK1_m')
plt.loglog(H, err_cn,   'y-d', label='Crank Nicolson')
plt.loglog(H, err_heun, 'c-v', label='Heun')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True);


# Calcul des ordres de convergence = pente des droites de regression
order_ee   = np.polyfit(np.log(H), np.log(err_ee), 1)[0]
order_ei   = np.polyfit(np.log(H), np.log(err_ei), 1)[0]
order_em   = np.polyfit(np.log(H), np.log(err_em), 1)[0]
order_rk1m = np.polyfit(np.log(H), np.log(err_rk1m), 1)[0]
order_cn   = np.polyfit(np.log(H), np.log(err_cn), 1)[0]
order_heun = np.polyfit(np.log(H), np.log(err_heun), 1)[0]


# Print des ordres de convergence

# Version avec print
# print (f'Euler progressif : ordre = {order_ee:1.2f}')
# print (f'Euler regressif  : ordre = {order_ei:1.2f}')
# print (f'Euler modifié  : ordre = {order_em:1.2f}')
# print (f'RK1_M  : ordre = {order_rk1m:1.2f}')
# print (f'Cranck Nicolson  : ordre = {order_cn:1.2f}')
# print (f'Heun  : ordre = {order_heun:1.2f}')

# Version avec tabulate pour une affichage ameliore
from tabulate import tabulate
T = []
T.append( ["Schéma", "Ordre"] )
T.append( ["Euler progressif", order_ee] )
T.append( ["Euler regressif", order_ei] )
T.append( ["Euler modifié", order_em] )
T.append( ["RK1_M", order_rk1m] )
T.append( ["Cranck Nicolson", order_cn] )
T.append( ["Heun", order_heun] )
display(tabulate(T,headers="firstrow",tablefmt="html",floatfmt=".2f"))

Schéma	Ordre
Euler progressif	1.05
Euler regressif	0.96
Euler modifié	2.08
RK1_M	2.00
Cranck Nicolson	2.01
Heun	2.08

4 Exercice : Fonction intégrale

Il existe des fonctions définies par des intégrales qu’on ne peut pas calculer explicitement. Il est pourtant possible de calculer des valeurs approchées de ces fonctions.

Pour \(x\in\left[0;\frac{\pi}{2}\right]\) calculer et afficher la table des valeurs et tracer le graphe de la fonction

\( x \mapsto f(x)=\int_0^x \sqrt{1-\frac{1}{4}\sin^2(t)}\mathrm{d}t \)

en approchant numériquement un problème de Cauchy (lequel?) avec la méthode d’Euler explicite.

Vérifier que \(f\left(\frac{\pi}{6}\right)\simeq0.51788193\) et \(f\left(\frac{\pi}{2}\right)\simeq1.46746221\).

Correction

Rappel sur la dérivée d’une fonction définie par une intégrale:
si \(F(x)=\displaystyle\int_{u(x)}^{v(x)}g(x,t)\mathrm{d}t\) alors \(F'(x)=\displaystyle g\left(x,v(x)\right)v'(x)-g\left(x,u(x)\right)u'(x)+ \int_{u(x)}^{v(x)} \frac{\partial}{\partial x} g(x,t)\mathrm{d}t\).

La fonction \(f\) est solution du problème de Cauchy

\( \begin{cases} y'(t)= \sqrt{1-\frac{1}{4}\sin^2(t)},\\ y(0)=0. \end{cases} \)

En effet, \(f(0)=0\) et si on intègre l’EDO entre \(0\) et \(x\) on a

\( f(x)=\int_0^x \sqrt{1-\frac{1}{4}\sin^2(t)}\mathrm{d}t=\int_0^x y'(t)\mathrm{d}t=y(x)-y(0)=y(x). \)

Dans le code ci-dessous:

On définit l’équation différentielle : phi est une lambda function qui contient la fonction mathématique \(\varphi(t, y)\) dépendant des variables \(t\) et \(y\) (même si \(y\) n’apparaît pas explicitement).
On définit la condition initiale \(t_0=0\) et \(y_0=0\).
On introduit la discrétisation: les nœuds d’intégration \([t_0,t_1,\dots,t_{N}]\) sont contenus dans le vecteur tt.

Comme nous devons éstimer ensuite \(f\) en \(\frac{\pi}{6}\) et \(\frac{\pi}{2}\), faisons en sort davoir ces deux valeurs dans nos noeuds de discrétisation.
Étant donné que

\(t_0=0\),
\(t_N=\frac{\pi}{2}\) et
\(t_i=t_0+ih\) avec \(h=\frac{t_N-t_0}{N}=\frac{\pi}{2N}\),

on a \(\frac{\pi}{6}=i\frac{\pi}{2N}\) ssi \(i=\frac{2N}{6}=\frac{N}{3}\) et bien sûr \(\frac{\pi}{2}=t_{N}\). Pour que \(i\in\mathbb{N}\) on choisira \(N\) multiple de \(3\).

On calcule enfin la solution approchée et on vérifie que \(f(\pi/6)\simeq0.51788193\) et \(f(\pi/2)\simeq1.46746221\) et on affiche le graphe de la fonction approchée.

Code

%reset -f
%matplotlib inline

import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 12}) # rcParams.update(rcParamsDefault)
import numpy as np


def EE(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt) - 1):
        uu[n+1] = uu[n] + h * phi(tt[n], uu[n])
    return uu

# EDO
phi = lambda t,y : np.sqrt(1-0.25*(np.sin(t))**2)

# CI
t0, y0 = 0, 0

# DISCRETISATION
tfinal = np.pi/2
k = 100
N = 3*k # ainsi t_N=pi/2=3pi/6 et t_k=pi/6
tt = np.linspace(t0,tfinal,N+1)

# CALCUL DE LA SOLUTION APPROCHÉE
uu = EE(phi,tt, y0)

print(f"f(pi/{np.pi/tt[k]}) = {uu[k]:g}") 
print(f"f(pi/{np.pi/tt[-1]}) = {uu[-1]:g}") 

# AFFICHAGE, NON DEMANDÉ
plt.plot(tt,uu)
plt.plot([tt[0],tt[k],tt[k]],[uu[k],uu[k],0],'m--o');
plt.plot([tt[0],tt[-1],tt[-1]],[uu[-1],uu[-1],0],'c--o');

f(pi/6.0) = 0.517965
f(pi/2.0) = 1.46781

Et sympy ?

Code

import sympy as sym
sym.init_printing()

t   = sym.Symbol('t')
y   = sym.Function('y')
edo = sym.Eq( sym.diff(y(t),t) , sym.sqrt(1-sym.Rational(1,4)*(sym.sin(t))**2) )
display(edo)
# solgen = sym.dsolve(edo)
# display(solgen)

t  = sym.Symbol('t')
f = sym.Integral(sym.sqrt(1-sym.Rational(1,4)*(sym.sin(t))**2))
display(f)
f.doit()

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = \sqrt{1 - \frac{\sin^{2}{\left(t \right)}}{4}}\)

\(\displaystyle \int \sqrt{1 - \frac{\sin^{2}{\left(t \right)}}{4}}\, dt\)

\(\displaystyle E\left(t\middle| \frac{1}{4}\right)\)

5 Exercice : Système de deux EDO avec invariant - influence de la méthode d’intégration

Il faut rester toujours critique par rapport aux résultats numériques obtenus. Les outils numériques ne pouvant fournir qu’une approximation de la solution, ils ne donnent aucune garantie que le résultat est forcément bon. Il faut donc penser à vérifier que les résultats obtenus s’approchent de ce à quoi on s’attend. Si la solution exacte n’est pas connue, il faudra se donner des critères pour juger de la qualité de la solution approchée.

Considérons le problème de Cauchy

\( \begin{cases} x'(t)= -y(t) ,\\ y'(t)= x(t), & t\in[0;4\pi]\\ x(0)=1,\\ y(0)=0. \end{cases} \)

Invariant
Vérifier que \(I(t)=x^2(t) + y^2(t)\) est un invariant pour le système.

Approximations
Dans une simulation numérique on aimerait que l’invariant soit préservé aussi bien que possible. Nous allons regarder ce qui se passe avec certaines méthodes vues en cours.
On notera \(u_n\approx x_n=x(t_n)\) et \(w_n\approx y_n=y(t_n)\).
À chaque instant \(t_n\), on calculera \(J_n\approx I_n=I(t_n)\).
Choisissons \(N=193\) points de discrétisation ainsi \(h=t_{n+1}-t_n=\pi/48\).

Dans un premier temps on se propose d’appliquer la méthode d’Euler explicite à la résolution du système.
1. Tracer \((t_n,u_n)\) et \((t_n,w_n)\) sur un même graphique;
2. sur un autre graphique tracer \((t_n,J_n)\);
3. tracer enfin \((u_n,w_n)\) sur un autre graphique et vérifier que la solution numérique tourne vers l’extérieur
4. Après avoir constaté que l’invariant \(J_n\) n’est pas conservé mais augment au cours du temps, montrer analytiquement que \(J_n=(1+h^2)^n\).
Même exercice pour le schéma d’Euler implicite (la linéarité du second membre permet de rendre ce schéma explicite). Que peut-on constater? Commenter les résultats et écrire l’expression analytique de \(J_n\).
Même exercice pour le schéma de Crank Nicolson (la linéarité du second membre permet de rendre ce schéma explicite). Que peut-on constater? Commenter les résultats et écrire l’expression analytique de \(J_n\).
Même exercice avec la fonction odeint du module scipy.optimize.
Bonus : même exercice avec la fonction solve_ivp du module scipy.optimize avec comparaison de différentes méthodes d’intégration.

Solution exacte
Dans ce cas il est même possible de calculer analytiquement la solution exacte. Calculer donc la solution exacte avec le module sympy.

Correction

Invariant

Pour tout \(t\) on a

\( I'(t)=\left(x^2(t) + y^2(t)\right)' = 2x'(t)x(t)+2y'(t)y(t)=-2x(t)y(t)+2x(t)y(t)=0 \)

donc \(I(t)=I(0)=0^2+1^2=1\).

Approximations

On notera \(u_n\approx x_n=x(t_n)\) et \(w_n\approx y_n=y(t_n)\).

À chaque instant \(t_n\), on calculera \(J_n=u_n^2+w_n^2\approx I_n=I(t_n)\).

Code

%reset -f
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 12}) 
# from scipy.optimize import fsolve # non necessaire car on resout l'equation implicite a la main

import sympy as sym
sym.init_printing()

Code

# variables globales

t0 = 0
tfinal = 4*np.pi

#h = np.pi/48
N = 192 #*100
tt = np.linspace(0, tfinal, N+1)  # np.arange(t0, tfinal+h/2, h)

x0 = 1
y0 = 0

phi1 = lambda t,x,y : -y
phi2 = lambda t,x,y : x

Code

def affichage(uu, ww, JJ, tt, x0, y0, phi1, phi2, schema):

    plt.figure(figsize=(24,7))

    plt.subplot(1,3,1)
    plt.plot(tt,uu,tt,ww)
    plt.grid()
    plt.xlabel('t')
    plt.legend(['x(t)','y(t)'])
    plt.title(f'{schema} - x(t) et y(t)') 

    plt.subplot(1,3,2)
    plt.plot(tt,JJ)
    plt.grid()
    plt.xlabel('t')
    plt.title(f'{schema} - Invariant')

    plt.subplot(1,3,3)
    u_min, u_max, w_min, w_max = min(uu), max(uu), min(ww), max(ww)
    Y1, Y2 = np.meshgrid( np.linspace(u_min,u_max,20), np.linspace(w_min,w_max,20) )
    V1, V2 = phi1(tt,Y1,Y2), phi2(tt,Y1,Y2)
    r = np.sqrt(V1**2+V2**2)
    plt.quiver(Y1, Y2, V1/r, V2/r) 
    plt.plot(uu, ww, 'r')
    plt.plot([x0], [y0], 'o', ms=10) # point initial
    plt.xlabel('x(t)')
    plt.ylabel('y(t)')
    plt.title(f'{schema} - y(x)')
    plt.axis('equal');

Euler explicite

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\varphi_1(t_n,u_n,w_n),\\ w_{n+1}=w_n+h\varphi_2(t_n,u_n,w_n) \end{cases} \)

Code

def EE(phi1,phi2,tt,x0,y0):
    uu = np.zeros_like(tt)
    ww = np.zeros_like(tt)
    uu[0] = x0
    ww[0] = y0
    h = tt[1]-tt[0]
    for n in range(len(tt)-1):
        uu[n+1] = uu[n]+h*phi1(tt[n],uu[n],ww[n])
        ww[n+1] = ww[n]+h*phi2(tt[n],uu[n],ww[n])
    return uu, ww

On remarque que

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n-hw_n,\\ w_{n+1}=w_n+hu_n \end{cases} \)

On a

\( J_{n+1} =u_{n+1}^2+w_{n+1}^2 =(u_n-hw_n)^2+(w_n+hu_n)^2 =(1+h^2)(u_n^2+w_n^2) =(1+h^2)J_n \)

soit encore

\( J_n=(1+h^2)^nJ_0=(1+h^2)^n. \)

On voit que

l’invariant n’est conservé que si \(h=0\), ce qui est impossible,
pour \(h\) fixé, \(\lim\limits_{n\to+\infty}J_n=+\infty\),
pour \(n\) fixé, l’erreur \(|J_n-I_n|\) diminue comme \(h^2\).

Code

uu_ep, ww_ep = EE(phi1,phi2,tt,x0,y0)
JJ_ep = uu_ep**2+ww_ep**2
affichage(uu_ep, ww_ep, JJ_ep, tt, x0, y0, phi1, phi2, 'Euler Explicite')

Euler implicite

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\varphi_1(t_{n+1},u_{n+1},w_{n+1})=u_n-hw_{n+1},\\ w_{n+1}=w_n+h\varphi_2(t_{n+1},u_{n+1},w_{n+1})=w_n+hu_{n+1} \end{cases} \)

En resolvant le système linéaire on obtient une écriture explicite

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=\dfrac{u_n-hw_{n}}{1+h^2},\\ w_{n+1}=\dfrac{w_n+hu_{n}}{1+h^2}. \end{cases} \)

Code

# Notons que même si ici il est inutile de passer phi1, phi2 car non utilisées explicitement dans le code,
# il est préférable de les passer pour une meilleure modularité du code (on ne change pas la signature de la fonction)

def EI(phi1,phi2,tt,x0,y0):
    uu = np.zeros_like(tt)
    ww = np.zeros_like(tt)
    uu[0] = x0
    ww[0] = y0
    h  = tt[1]-tt[0]
    for n in range(len(tt)-1):
        uu[n+1] = (uu[n]-h*ww[n])/(1+h**2)
        ww[n+1] = (ww[n]+h*uu[n])/(1+h**2)
    return  uu, ww

Code

#  VERSION AVEC fsolve (non demandé)
from scipy.optimize import fsolve

def EI(phi1,phi2,tt,x0,y0):
    uu = np.zeros_like(tt)
    ww = np.zeros_like(tt)
    uu[0] = x0
    ww[0] = y0
    h  = tt[1]-tt[0]
    for n in range(len(tt)-1):
        sys = lambda z : [ -z[0]+uu[n]+h*phi1(tt[n+1],z[0],z[1]) , 
                           -z[1]+ww[n]+h*phi2(tt[n+1],z[0],z[1]) ]
        utemp,wtemp = fsolve( sys , (uu[n],ww[n]) ) 
        uu[n+1] = utemp
        ww[n+1] = wtemp
    return uu, ww

On a

\( J_{n+1} =u_{n+1}^2+w_{n+1}^2 =\frac{(u_n-hw_n)^2+(w_n+hu_n)^2}{(1+h^2)^2} =(u_n^2+w_n^2)\frac{1+h^2}{(1+h^2)^2} =\frac{1+h^2}{(1+h^2)^2}J_n =\frac{1}{1+h^2}J_n \)

soit encore

\( J_n=\left(\frac{1}{1+h^2}\right)^nJ_0=\left(\frac{1}{1+h^2}\right)^n. \)

On voit que

l’invariant n’est conservé que si \(h=0\), ce qui est impossible,
pour \(h\) fixé, \(\lim\limits_{n\to+\infty}J_n=0^+\),
pour \(n\) fixé, l’erreur \(|J_n-I_n|\) diminue comme \(h^2\).

Code

[uu_er, ww_er] = EI(phi1,phi2,tt,x0,y0)
JJ_er = uu_er**2+ww_er**2
affichage(uu_er, ww_er, JJ_er, tt, x0, y0, phi1, phi2, 'Euler Implicite')

Crank Nicolson

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+\frac{h}{2}\Big(\varphi_1(t_{n},u_{n},w_{n})+\varphi_1(t_{n+1},u_{n+1},w_{n+1})\big)=u_n-\frac{h}{2}w_{n}-\frac{h}{2}w_{n+1},\\ w_{n+1}=w_n+\frac{h}{2}\Big(\varphi_2(t_{n},u_{n},w_{n})+\varphi_2(t_{n+1},u_{n+1},w_{n+1})\big)=w_n+\frac{h}{2}u_{n}+\frac{h}{2}u_{n+1}. \end{cases} \)

En resolvant le système linéaire on obtient une écriture explicite:

\( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=\dfrac{(4-h^2)u_n-4hw_{n}}{4+h^2},\\ w_{n+1}=\dfrac{(4-h^2)w_n+4hu_{n}}{4+h^2}. \end{cases} \)

De plus, on a

\( \begin{aligned} J_{n+1} &=u_{n+1}^2+w_{n+1}^2 \\ &=\frac{((4-h^2)u_n-4hw_n)^2+((4-h^2)w_n+4hu_n)^2}{(4+h^2)^2} \\ &=\frac{(4-h^2)^2u_n^2+16h^2w_n^2-8(4-h^2)hu_nw_n+(4-h^2)^2w_n^2+16h^2u_n^2+8(4-h^2)hu_nw_n}{(4+h^2)^2} \\ &=\frac{\big((4-h^2)^2+16h^2\big)(u_n^2+w_n^2)}{(4+h^2)^2} \\ &=\frac{(4+h^2)^2(u_n^2+w_n^2)}{(4+h^2)^2} \\ &=(u_n^2+w_n^2)=J_n \end{aligned} \)

On voit que l’invariant est conservé pour tout \(h\)

Vérifions ces calculs avec sympy:

Code

dt = sym.Symbol('h') # pour ne pas redefinir h comme un symbole, on écrit dt mais on affiche h
u_n = sym.Symbol('u_n')
u_np1 = sym.Symbol('u_{n+1}')
w_n = sym.Symbol('w_n')
w_np1 = sym.Symbol('w_{n+1}')

eq1 = sym.Eq( u_np1 , u_n-dt/2*w_n-dt/2*w_np1 )
eq2 = sym.Eq( w_np1 , w_n+dt/2*u_n+dt/2*u_np1 )
sol = sym.solve([eq1,eq2],[u_np1,w_np1],dict=True)[0]
display(sol)

J_np1 = u_np1**2+w_np1**2
J_np1.subs(sol).simplify()

\(\displaystyle \left\{ u_{n+1} : \frac{- h^{2} u_{n} - 4 h w_{n} + 4 u_{n}}{h^{2} + 4}, \ w_{n+1} : \frac{- h^{2} w_{n} + 4 h u_{n} + 4 w_{n}}{h^{2} + 4}\right\}\)

\(\displaystyle u_{n}^{2} + w_{n}^{2}\)

Code

def cn(phi1,phi2,tt,x0,y0):
    uu = np.zeros_like(tt)
    ww = np.zeros_like(tt)
    uu[0] = x0
    ww[0] = y0
    h  = tt[1]-tt[0]
    for n in range(len(tt)-1):
        uu[n+1] = ((4-h**2)*uu[n]-4*h*ww[n])/(4+h**2)
        ww[n+1] = ((4-h**2)*ww[n]+4*h*uu[n])/(4+h**2)
    return uu, ww

Code

#  VERSION AVEC fsolve (non demandé)
from scipy.optimize import fsolve

def cn(phi1,phi2,tt,x0,y0):
    uu = np.zeros_like(tt)
    ww = np.zeros_like(tt)
    uu[0] = x0
    ww[0] = y0
    h  = tt[1]-tt[0]
    for n in range(len(tt)-1):
        sys = lambda z : [ -z[0]+uu[n]+h/2*(phi1(tt[n],uu[n],ww[n])+phi1(tt[n+1],z[0],z[1])) , 
                           -z[1]+ww[n]+h/2*(phi2(tt[n],uu[n],ww[n])+phi2(tt[n+1],z[0],z[1])) ]
        utemp,wtemp = fsolve( sys , (uu[n],ww[n]) ) 
        uu[n+1] = utemp
        ww[n+1] = wtemp
    return uu, ww

Code

[uu_cn, ww_cn] = cn(phi1,phi2,tt,x0,y0)
JJ_cn = uu_cn**2+ww_cn**2
affichage(uu_cn, ww_cn, JJ_cn, tt, x0, y0, phi1, phi2, 'Crank Nicolson')

Avec odeint du module scipy.integrate:

Code

from scipy.integrate import odeint

pphi_odeint = lambda t,yy : [ phi1(t,yy[0],yy[1]) , phi2(t,yy[0],yy[1]) ]
uu_odeint, ww_odeint = odeint(pphi_odeint,[x0,y0],tt,tfirst=True).T

JJ_odeint = uu_odeint**2+ww_odeint**2

affichage(uu_odeint, ww_odeint, JJ_odeint, tt, x0, y0, phi1, phi2, 'ODEINT')

Bonus : avec solve_ivp du module scipy.integrate:

Code

from scipy.integrate import solve_ivp

pphi_solve_ivp = lambda t,yy : [ phi1(t,yy[0],yy[1]) , phi2(t,yy[0],yy[1]) ]

Code

sol_RK45 = solve_ivp(pphi_solve_ivp,[t0,tfinal],[x0,y0],t_eval=tt, method='RK45')
uu_RK45, ww_RK45 = sol_RK45.y

JJ_RK45 = uu_RK45**2+ww_RK45**2

affichage(uu_RK45, ww_RK45, JJ_RK45, tt, x0, y0, phi1, phi2, 'RK45')

Code

sol_Radau = solve_ivp(pphi_solve_ivp,[t0,tfinal],[x0,y0],t_eval=tt, method='Radau')
uu_Radau, ww_Radau = sol_Radau.y

JJ_Radau = uu_Radau**2+ww_Radau**2

affichage(uu_Radau, ww_Radau, JJ_Radau, tt, x0, y0, phi1, phi2, 'Radau')

Code

sol_BDF = solve_ivp(pphi_solve_ivp,[t0,tfinal],[x0,y0],t_eval=tt, method='BDF')
uu_BDF, ww_BDF = sol_BDF.y

JJ_BDF = uu_BDF**2+ww_BDF**2

affichage(uu_BDF, ww_BDF, JJ_BDF, tt, x0, y0, phi1, phi2, 'BDF')

Code

sol_LSODA = solve_ivp(pphi_solve_ivp,[t0,tfinal],[x0,y0],t_eval=tt, method='LSODA')
uu_LSODA, ww_LSODA = sol_LSODA.y

JJ_LSODA = uu_LSODA**2+ww_LSODA**2

affichage(uu_LSODA, ww_LSODA, JJ_LSODA, tt, x0, y0, phi1, phi2, 'LSODA')

Solution exacte

La solution exacte est

\(\begin{align*} x(t)&=\cos(t),\\ y(t)&=\sin(t). \end{align*}\)

Cette solution est périodique de période \(2\pi\) et on a bien

\( I(t)=x^2(t) + y^2(t) = 1 \quad \forall t. \)

Vérifions-le avec le module sympy:

Code

t = sym.Symbol('t')
x = sym.Function('x')
y = sym.Function('y')

t_0 = 0
x_0 = 1
y_0 = 0


edo1 = sym.Eq( sym.diff(x(t),t) , phi1(t,x(t),y(t)) )
edo2 = sym.Eq( sym.diff(y(t),t) , phi2(t,x(t),y(t)) )
display(edo1)
display(edo2)

# solgen = sym.dsolve([edo1,edo2],[x(t),y(t)])
# display(solgen)

# cond_init_1 = sym.Eq( x_0, solgen[0].rhs.subs(t,t_0))
# cond_init_2 = sym.Eq( y_0, solgen[1].rhs.subs(t,t_0))
# consts = sym.solve( [ cond_init_1 , cond_init_2 ] , dict=True)[0]
# display(consts)

# solpar_1 = solgen[0].subs(consts)
# solpar_2 = solgen[1].subs(consts)
# display(solpar_1)
# display(solpar_2)

solpar_1, solpar_2 = sym.dsolve([edo1,edo2],[x(t),y(t)], ics={x(t_0): x_0, y(t_0): y_0})


# AFFICHAGE
func_1 = sym.lambdify(t,solpar_1.rhs,'numpy')
func_2 = sym.lambdify(t,solpar_2.rhs,'numpy')

tt = np.linspace(0,4*np.pi,193)
xx = func_1(tt)
yy = func_2(tt)
JJ = xx**2+yy**2

affichage(xx, yy, JJ, tt, x0, y0, phi1, phi2, 'Exacte')

\(\displaystyle \frac{d}{d t} x{\left(t \right)} = - y{\left(t \right)}\)

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = x{\left(t \right)}\)

Retour au sommet