from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './customNEW.css'
HTML(open(css_file, "r").read())

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

Python version 3.10.6 (main, Nov 14 2022, 16:10:14) [GCC 11.3.0]

M62 TP 2 - Implémentation des schémas "classiques" et étude de la convergence¶

1 Rappel de CM : implémentation des schémas d'Euler explicite et implicite
2 Rappel de CM : étude empirique de la convergence des schémas d'Euler explicite et implicite
3 Exercice : implémentation et étude empirique de la convergence des méthodes d'Euler modifié, RK1_m, de Crank-Nicolson et de Heun
4 Exercice : Fonction intégrale
5 Exercice : Système de deux EDO avec invariant

Rappel de CM : 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:

%reset -f
%matplotlib inline

import sympy as sym
sym.init_printing()

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)

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

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

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:

from matplotlib.pylab import *
rcParams.update({'font.size': 16}) # rcParams.update(rcParamsDefault)
from scipy.optimize import fsolve # pour les schémas implicites

On initialise le problème de Cauchy

# 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$.

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} $$

def EE(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = [y0]
    for i in range(len(tt) - 1):
        uu.append(uu[i] + h * phi(tt[i], uu[i]))
    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)$$

def EI(phi, tt, y0):
    h = tt[1] - tt[0]
    uu = [y0]
    for i in range(len(tt) - 1):
        temp = fsolve(lambda x: -x + uu[i] + h * phi(tt[i + 1], x), uu[i])
        uu.append(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:

N  = 8 
tt = linspace(t0,tfinal,N+1)
# print(tt)

yy    = [sol_exacte(t) for t in tt]
# 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| $$

figure(figsize=(18, 7))

subplot(1, 2, 1)
plot(tt, yy, 'b-', tt, uu_ep, 'r-D')
err_ep = norm(array(uu_ep)-array(yy),inf)
title(f'Euler explicite - max(|erreur|)= {err_ep:1.10f}')
grid()

subplot(1, 2, 2)
plot(tt, yy, 'b-', tt, uu_er, 'r-D')
err_er = norm(array(uu_er)-array(yy),inf)
title(f'Euler implicite - max(|erreur|)= {err_er:1.10f}')
grid();

Rappel de CM : é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}|$.

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

N = 10
for k in range(7):
    N+=20
    tt = linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = [sol_exacte(t) for t in tt]
    uu_ep = EE(phi, tt, y0)
    uu_er = EI(phi, tt, y0)
    H.append(h)
    err_ep.append(norm(array(uu_ep)-array(yy),inf))
    err_er.append(norm(array(uu_er)-array(yy),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})$.

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

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).

order_ep = polyfit(log(H),log(err_ep), 1)[0]
print (f'Euler progressif : ordre = {order_ep:1.2f}')
order_er = polyfit(log(H),log(err_er), 1)[0]
print (f'Euler regressif  : ordre = {order_er:1.2f}')

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

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 = h\varphi(t_n,u_n),\\ u_{n+1}=u_n+h\varphi\left(t_n+\frac{h}{2}, u_n+\frac{k_1}{2}\right)& n=0,1,2,\dots N-1 \end{cases} $$

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

Schéma RK1_M: $$ \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi\left(\frac{t_n+t_{n+1}}{2},\frac{u_n+u_{n+1}}{2}\right)& n=0,1,2,\dots N-1 \end{cases} $$

def RK1_M(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        uu.append( fsolve(lambda x: -x+uu[i]+h*phi( (tt[i]+tt[i+1])/2,(uu[i]+x)/2 ), uu[i])[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))$.

def CN(phi,tt, y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        temp = fsolve(lambda x: -x+uu[i]+0.5*h*( phi(tt[i+1],x)+phi(tt[i],uu[i]) ), uu[i])
        uu.append(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 = h\varphi(t_n,u_n),\\ k_2 = h\varphi(t_{n+1},u_n+k_1),\\ u_{n+1}=u_n+\frac{k_1+k_2}{2}& n=0,1,2,\dots N-1 \end{cases} $$

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

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

yy      = [sol_exacte(t) for t in tt]
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)

figure(1, figsize=(18, 10))

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

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


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

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

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

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

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

for k in range(7):
    N = 2**(k + 3)
    tt = linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = [sol_exacte(t) for t in tt]
    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(max([abs(uu_ee[i] - yy[i]) for i in range(len(yy))]))
    err_ei.append(max([abs(uu_ei[i] - yy[i]) for i in range(len(yy))]))
    err_em.append(max([abs(uu_em[i] - yy[i]) for i in range(len(yy))]))
    err_rk1m.append(max([abs(uu_rk1m[i] - yy[i]) for i in range(len(yy))]))
    err_cn.append(max([abs(uu_cn[i] - yy[i]) for i in range(len(yy))]))
    err_heun.append(max([abs(uu_heun[i] - yy[i]) for i in range(len(yy))]))

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

order_ee = polyfit(log(H),log(err_ee), 1)[0]
print (f'Euler progressif : ordre = {order_ee:1.2f}')
order_ei = polyfit(log(H),log(err_ei), 1)[0]
print (f'Euler regressif  : ordre = {order_ei:1.2f}')
order_em = polyfit(log(H),log(err_em), 1)[0]
print (f'Euler modifié  : ordre = {order_em:1.2f}')
order_rk1m = polyfit(log(H),log(err_rk1m), 1)[0]
print (f'RK1_M  : ordre = {order_rk1m:1.2f}')
order_cn = polyfit(log(H),log(err_cn), 1)[0]
print (f'Cranck Nicolson  : ordre = {order_cn:1.2f}')
order_heun = polyfit(log(H),log(err_heun), 1)[0]
print (f'Heun  : ordre = {order_heun:1.2f}')

Euler progressif : ordre = 1.05
Euler regressif  : ordre = 0.96
Euler modifié  : ordre = 2.08
RK1_M  : ordre = 2.00
Cranck Nicolson  : ordre = 2.01
Heun  : ordre = 2.08

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$.

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$.

Correction

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.

%reset -f
%matplotlib inline

from matplotlib.pylab import *
rcParams.update({'font.size': 16}) # rcParams.update(rcParamsDefault)

# schéma
def EE(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = [y0]
    for i in range(len(tt) - 1):
        uu.append(uu[i] + h * phi(tt[i], uu[i]))
    return uu

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

# CI
t0, y0 = 0, 0

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

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

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

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

f(pi/6.0)=0.5179648116355586
f(pi/2.0)=1.4678129540121247

Et sympy ?

%reset -f
%matplotlib inline

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()

Exercice : Système de deux EDO avec invariant¶

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.

Tracer $(t_n,u_n)$ et $(t_n,w_n)$ sur un même graphique;

sur un autre graphique tracer $(t_n,J_n)$;

tracer enfin $(u_n,w_n)$ sur un autre graphique et vérifier que la solution numérique tourne vers l’extérieur

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.

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)$.

%reset -f
%matplotlib inline

from matplotlib.pylab import *
rcParams.update({'font.size': 16}) # rcParams.update(rcParamsDefault)
# from scipy.optimize import fsolve # non necessaire car on resout l'equation implicite a la main

t0 = 0
tfinal = 4*pi

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

x0 = 1
y0 = 0

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

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

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

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$.

[uu_ep, ww_ep] = EE(phi1,phi2,tt,x0,y0)
JJ_ep = [(uu_ep[i])**2+(ww_ep[i])**2 for i in range(len(tt))]


figure(figsize=(24,7))

subplot(1,3,1)
plot(tt,uu_ep,tt,ww_ep)
grid()
xlabel('t')
legend(['x(t)','y(t)'])
title('Euler progressif - x(t) et y(t)') 

subplot(1,3,2)
plot(tt,JJ_ep)
grid()
xlabel('t')
title('Euler progressif - Invariant')

subplot(1,3,3)
Y1,Y2 = meshgrid(linspace(min(uu_ep),max(uu_ep),21),linspace(min(ww_ep),max(ww_ep),21))
V1,V2 = phi1(tt,Y1,Y2),phi2(tt,Y1,Y2)
r=sqrt(V1**2+V2**2)
quiver(Y1, Y2, V1/r, V2/r) 
plot(uu_ep,ww_ep,'r')
xlabel('x(t)')
ylabel('y(t)')
title('Euler progressif - y(x)');

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} $$

# Notons que ici il est inutile de passer phi1,phi2 car non utilisées explicitement dans le code
def EI(phi1,phi2,tt,x0,y0):
	uu = [x0]
	ww = [y0]
	h  = tt[1]-tt[0]
	for i in range(len(tt)-1):
		uu.append((uu[i]-h*ww[i])/(1+h**2))
		ww.append((ww[i]+h*uu[i])/(1+h**2))
	return [uu,ww]

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

def EI(phi1,phi2,tt,x0,y0):
	uu = [x0]
	ww = [y0]
	h  = tt[1]-tt[0]
	for i in range(len(tt)-1):
		sys = lambda z : [ -z[0]+uu[i]+h*phi1(tt[i+1],z[0],z[1]) , -z[1]+ww[i]+h*phi2(tt[i+1],z[0],z[1]) ]
		utemp,wtemp = fsolve( sys , (uu[i],ww[i]) ) 
		uu.append(utemp)
		ww.append(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$.

[uu_er, ww_er] = EI(phi1,phi2,tt,x0,y0)
JJ_er = [(uu_er[i])**2+(ww_er[i])**2 for i in range(len(tt))]

figure(figsize=(24,7))

subplot(1,3,1)
plot(tt,uu_er,tt,ww_er)
grid()
xlabel('t')
legend(['x(t)','y(t)'])
title('Euler regressif - x(t) et y(t)') 

subplot(1,3,2)
plot(tt,JJ_er)
grid()
xlabel('t')
title('Euler regressif - Invariant')

subplot(1,3,3)
Y1,Y2 = meshgrid(linspace(min(uu_er),max(uu_er),21),linspace(min(ww_er),max(ww_er),21))
V1,V2 = phi1(tt,Y1,Y2),phi2(tt,Y1,Y2)
r = sqrt(V1**2+V2**2)
quiver(Y1, Y2, V1/r, V2/r) 
# streamplot(Y1, Y2, V1/r, V2/r) 
plot(uu_er,ww_er)
xlabel('x(t)')
ylabel('y(t)')
title('Euler regressif - y(x)');