In [5]:
import math
from matplotlib.pylab import *
%matplotlib inline

R43 - Devoir maison (mai 2018)

Considérons le problème de Cauchy : trouver une fonction $y \colon I\subset \mathbb{R} \to \mathbb{R}$ définie sur un intervalle $I=[t_0,T]$ telle que $$ \begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in I=[t_0,T],\\ y(t_0) = y_0, \end{cases} $$ avec $y_0$ une valeur donnée et supposons que l'on ait montré l'existence et l'unicité d'une solution $y$ pour $t\in I$.

On subdivise l'intervalle $I=[t_0;T]$, avec $T<+\infty$, en $N$ intervalles $[t_{n};t_{n+1}]$ de largeur $h=(T-t_0)/N$ avec $t_n=t_0+nh$ pour $n=0,\dots,N$.

Pour chaque nœud $t_n$, on note $y_n=y(t_n)$ et on cherche la valeur inconnue $u_n$ qui approche la valeur exacte $y_n$

Exercice 1

Considérons le problème de Cauchy $$ \begin{cases} y'(t)= \begin{cases}y(t)+t&\text{si }-1\le t\le 0,\\ y(t)-t&\text{si }0\le t\le 1, \end{cases} \\ y(-1)=1. \end{cases} $$

  1. Tracer le graphe de la solution approchée obtenue avec le schéma Nyström-2 avec un pas de discrétisation $h=0.01$.
  2. Calculer la solution exacte et calculer l'erreur maximale commise entre la solution exacte et la solution approchée.

Solution exercice 1

On définit l'équation différentielle : phi est une fonction python qui contient la fonction mathématique $\varphi(t, y)$, définie par morceaux, dépendant des variables $t$ et $y$.

In [6]:
def phi(t,y):
	if t<=0.0:
		return y+t
	else:
		return y-t

et la condition initiale

In [7]:
t0,y0,tfinal = -1.0, 1.0, 1.0

On calcule la solution exacte $ y(t)=\begin{cases}y_L(t)&\text{si }-1\le t\le 0,\\ y_R(t)&\text{si }0\le t\le 1, \end{cases} $

  • Sur $[-1;0]$ on a $a(t)=1$, $b(t)=-1$, $g(t)=t$,

    • $A(t)=\int \frac{b(t)}{a(t)}\mathrm{d}t=\int -1\mathrm{d}t=-t$,
    • $B(t)=\int \frac{g(t)}{a(t)}e^{A(t)}\mathrm{d}t=-\int -te^{-t}\mathrm{d}t=-(1+t)e^{-t}$,

    donc pour $t\in [-1;0]$ on a $y_L(t)=ce^{-A(t)}+B(t)e^{-A(t)}=ce^{t}-(1+t)e^{-t}e^{t}=ce^{t}-1-t$.

    La condition initiale $1=y_L(-1)=c/e$ impose $c=e$ ainsi $y_L(t)=e^{t+1}-1-t$.

  • On recommence sur $[0;1]$: on a $a(t)=1$, $b(t)=-1$, $g(t)=-t$,

    • $A(t)=\int \frac{b(t)}{a(t)}\mathrm{d}t=\int -1\mathrm{d}t=-t$,
    • $B(t)=\int \frac{g(t)}{a(t)}e^{A(t)}\mathrm{d}t=\int -te^{-t}\mathrm{d}t=(1+t)e^{-t}$,

    donc pour $t\in [0;1]$ on a $y_R(t)=ce^{t}+(1+t)e^{-t}e^{t}=ce^{t}+1+t$.

    La condition initiale est obtenue en imposant la continuité avec $y_L$: $y_L(0)=e-1=y_R(0)=c+1$ ce qui impose $c=e-2$ ainsi $y_R(t)=e^{t+1}-2e^t+1+t$.

Conclusion: la solution est $$ y(t)= \begin{cases} e^{t+1}-1-t&\text{si }-1\le t\le 0,\\ e^{t+1}-2e^t+1+t&\text{si }0\le t\le 1. \end{cases} $$ La fonction ainsi calculée vérifie l'EDO, la condition initiale, et est bien de classe $\mathcal{C}^1([-1;1])$.

In [8]:
def exact(t):
	if t<=0.0:
		return math.exp(t+1.0)-t-1.0
	else:
		return math.exp(t+1.0)-2.0*math.exp(t)+t+1.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

In [9]:
h=0.01
tt=arange(t0,tfinal+0.1,h)

On écrit le schéma numérique : les valeurs $[u_0,u_1,\dots,u_{N}]$ pour la méthode sont contenues dans le vecteur uu

Schéma de Nylström d'ordre 2: $$ \begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi(t_0,u_0),\\ u_{n+1}=u_{n-1}+2h\varphi(t_{n},u_{n})& n=1,2,3,4,5,\dots N-1 \end{cases} $$

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

On calcule les solutions exacte et approchée:

In [11]:
uu=N2(phi,tt)
yy=[exact(t) for t in tt]

On compare les graphes des solutions exacte et approchée et on affiche le maximum de l'erreur:

In [12]:
err=[abs(yy[i]-uu[i]) for i in range(len(yy))]
print('Erreur max =',max(err))

plot(tt,uu,'b-o')
show()   
Erreur max = 0.00023819339299091524

Exercice 2

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[0;\pi/2]$ 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 $$ avec un pas de discrétisation $h=\dfrac{\pi}{48}$ en approchant numériquement un problème de Cauchy (lequel?) avec la méthode de Heun.

Vérifier que $f(\pi/6)\simeq0.51788193$ et $f(\pi/2)\simeq1.46746221$.

Solution exercice 2

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

On définit l'équation différentielle : phi est une fonction python qui contient la fonction mathématique $\varphi(t, y)$ dépendant des variables $t$ et $y$.

In [13]:
def phi(t,y):
	return sqrt(1.0-0.25*(sin(t))**2)

et la condition initiale

In [14]:
t0=0.
y0=0.
tfinal = math.pi/2

On introduit la discrétisation: les nœuds d'intégration $[t_0,t_1,\dots,t_{N}]$ sont contenus dans le vecteur tt

In [15]:
N=25
tt=linspace(t0,tfinal,N)

On écrit le schéma numérique : les valeurs $[u_0,u_1,\dots,u_{N}]$ pour la méthode sont contenues dans le vecteur 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} $$

In [16]:
def heun(phi,tt):
    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.0 )
    return uu

On calcule la solution approchée:

In [17]:
uu=heun(phi,tt)

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.

In [18]:
print("u(%1.2f pi)=%1.10f" %(tt[8]/math.pi, uu[8]))
print("u(%1.2f pi)=%1.10f" %(tt[-1]/math.pi, uu[-1]))

plot(tt,uu)
show() 
u(0.17 pi)=0.5178420139
u(0.50 pi)=1.4674622093

Exercice 3 - problème mal conditionné (= numériquement mal posé)

Considérons le problème de Cauchy $$ \begin{cases} y'(t)= 10y(t)+11t-5t^2-1, &t\in[0;3]\\ y(0)=\varepsilon. \end{cases} $$ Calculer la solution exacte en fonction de la donnée initiale.

Pour $\varepsilon=0$, tracer les graphes de la solution exacte et de la solution approchée obtenue avec le schéma RK4 avec un pas de discrétisation $h=\dfrac{1}{2^8}$. Expliquer les résultats obtenus.

Solution exercice 3

On définit l'équation différentielle : phi est une fonction python qui contient la fonction mathématique $\varphi(t, y)$ dépendant des variables $t$ et $y$.

In [32]:
def phi(t,y):
	return 10*y+11*t-5*t**2-1

et la condition initiale

In [33]:
t0=0.0
y0=0.0
tfinal = 3.0

On calcule la solution exacte: on a $a(t)=1$, $b(t)=-10$ et $g(t)=11t-5t^2-1$ donc

  • $A(t)=\int -10dt= -10t$,
  • $B(t)=\int (11t-5t^2-1) e^{-10t}dt=\frac{t(t-2)}{2}e^{-10t}$,

d'où $y(t)=ce^{10 t}+\frac{t(t-2)}{2}$.

En imposant la condition $\varepsilon=y(0)$ on trouve l'unique solution du problème de Cauchy donné: $$ y(t)=\varepsilon e^{10 t}+\frac{t(t-2)}{2}. $$ On constate qu'une petite perturbation de la condition initiale entraîne une grosse perturbation de la solution, en effet si $\varepsilon=0$ la solution est une parabole, si $\varepsilon\neq0$ le terme en $e^{10t}$ est dominant.

Par exemple $$ y(4) =\varepsilon e^{40}+4 =\begin{cases} 4 &\text{si }\varepsilon=0,\

10^7&\text{si }\varepsilon=10^{-10}. \end{cases} $$ Par conséquent, si $\varepsilon=0$ le problème de Cauchy est numériquement mal posé car toute (petite) erreur de calcul a le même effet qu'une perturbation de la condition initiale: on *réveil* le terme dormant $e^{10t}$.

In [34]:
def exact(t):
	return y0*math.exp(10.*t)+0.5*(t**2)-t

On introduit la discrétisation: les nœuds d'intégration $[t_0,t_1,\dots,t_{N}]$ sont contenus dans le vecteur tt

In [35]:
h=1./(2.0**8)
tt=arange(t0,tfinal+0.1,h)

On écrit le schéma numérique : les valeurs $[u_0,u_1,\dots,u_{N}]$ pour la méthode sont contenues dans le vecteur uu

Schéma de Runke-Kutta RK4: $$\begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n+\frac{h}{2} \varphi(t_{n},u_{n}),\\ \check u_{n+1/2}=u_n+\frac{h}{2} \varphi(t_{n}+\frac{h}{2},\tilde u_{n+1/2}),\\ \hat u_{n+1}=u_n+h\varphi(t_{n+1},\check u_{n+1/2}),\\ u_{n+1}=u_n+\frac{h}{6}\left(\varphi(t_{n},u_{n})+2\varphi(t_{n}+\frac{h}{2},\tilde u_{n+1/2} )+2\varphi(t_{n}+\frac{h}{2}, \check u_{n+1/2})+\varphi(t_{n+1},\hat u_{n+1} \right)& n=0,1,2,\dots N-1 \end{cases} $$

In [36]:
def RK4(phi,tt):
    uu = [y0]
    for i in range(len(tt)-1):
        k1 = h * phi( tt[i], uu[i] )
        k2 = h * phi( tt[i]+h/2. , uu[i]+k1/2. )
        k3 = h * phi( tt[i]+h/2. , uu[i]+k2/2. )
        k4 = h * phi( tt[i+1]    , uu[i]+k3 )
        uu.append( uu[i] + (k1+2.0*k2+2.0*k3+k4) /6.0 )
    return uu

On calcule les solutions exacte et approchée:

In [37]:
yy=[exact(t) for t in tt]
uu=RK4(phi,tt) 

On compare les graphes des solutions exacte et approchée et on affiche le maximum de l'erreur:

In [38]:
plot(tt,yy,'r-',tt,uu,'b-o')
show()  

Exercice 4 - Système de deux EDO avec invariant

Considérons le problème de Cauchy $$ \begin{cases} x'(t)= \varphi_1(t,x,y),\\ y'(t)= \varphi_2(t,x,y), & t\in[t_0;T]\\ x(t_0)=x_0,\\ y(t_0)=y_0. \end{cases} $$ avec \begin{align*} t_0&=0, & T&=4\pi, & \varphi_1(t,x,y)&= -y(t), & \varphi_2(t,x,y)&=x(t), & x_0&=1, & y_0&=0. \end{align*} 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 possède un invariant: $$ I(t)=x^2(t) + y^2(t) = 1, \qquad \forall t. $$ Dans une simulation numérique on aimerait que cette propriété soit préservée 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 $h=t_{n+1}-t_n=\pi/48$.

  1. 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$.
  2. Même exercice pour le schéma d'Euler implicite (la linéarité du second membre permet de rendre se schéma explicite). Que peut-on constater? Commenter les résultats et écrire l'expression analytique de $J_n$.
  3. Même exercice pour le schéma de Crank-Nicolson (la linéarité du second membre permet de rendre se schéma explicite). Que peut-on constater? Commenter les résultats et montrer analytiquement que $J_n=1$ pour tout $n$.
  4. Même exercice pour le schéma symplectique obtenu en écrivant le schéma d'Euler explicite pour la première équation et le schéma d'Euler implicite pour la seconde. Que peut-on constater? Commenter les résultats et donner l'expression analytique de $J_n$.
  5. Même exercice pour le schéma obtenu en écrivant le schéma d'Euler explicite pour la première équation et le schéma de Crank-Nicolson pour la seconde. Que peut-on constater? Commenter les résultats et donner l'expression analytique de $J_n$.

Solution exercice 4

In [22]:
t0 = 0.
tfinal = 4.*math.pi

x0 = 1.
y0 = 0.

def phi1(t,x,y):
	return -y

def phi2(t,x,y):
	return x

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 possède un invariant : $$ I(t)= x^2(t) + y^2(t) = 1, \qquad \forall t. $$ Dans une simulation numérique on aimerait que cette propriété soit préservée 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=u_n^2+w_n^2\approx I_n=I(t_n)$.

Choisissons $h=t_{n+1}-t_n=\pi/48$.

In [23]:
h = math.pi/96
tt = arange(t0,tfinal+h,h)

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

In [24]:
def euler_progressif(phi1,phi2,tt):
	uu = [x0]
	ww = [y0]
	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)(u_n^2+w_n^2) =(1+h)J_n $$ soit encore $$ J_n=(1+h)^nJ_0=(1+h)^n. $$

In [25]:
[uu_ep, ww_ep] = euler_progressif(phi1,phi2,tt)
JJ_ep = [(uu_ep[i])**2+(ww_ep[i])**2 for i in range(len(tt))]

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

figure()
plot(tt,JJ_ep)
xlabel('t')
title('Euler progressif - Invariant')

figure()
plot(uu_ep,ww_ep)
xlabel('x(t)')
ylabel('y(t)')
title('Euler progressif - y(x)') 
Out[25]:
Text(0.5,1,'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} $$

In [26]:
def euler_regressif(phi1,phi2,tt):
	uu = [x0]
	ww = [y0]
	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]

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}{(1+h^2)^2} =\frac{1+h}{(1+h^2)^2}J_n $$ soit encore $$ J_n=\left(\frac{1+h}{(1+h^2)^2}\right)^nJ_0=\left(\frac{1+h}{(1+h^2)^2}\right)^n. $$

In [27]:
[uu_er, ww_er] = euler_regressif(phi1,phi2,tt)
JJ_er = [(uu_er[i])**2+(ww_er[i])**2 for i in range(len(tt))]

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

figure()
plot(tt,JJ_er)
xlabel('t')
title('Euler regressif - Invariant')

figure()
plot(uu_er,ww_er)
xlabel('x(t)')
ylabel('y(t)')
title('Euler regressif - y(x)') 
Out[27]:
Text(0.5,1,'Euler regressif - y(x)')

Crank-Nicolson $$ \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+\dfrac{h}{2}\left(\varphi_1(t_{n+1},u_{n+1},w_{n+1})+\varphi_1(t_{n},u_{n},w_{n})\right) =u_n-\dfrac{h}{2}\left(w_{n+1}+w_n\right),\\ w_{n+1}=w_n+\dfrac{h}{2}\left(\varphi_2(t_{n+1},u_{n+1},w_{n+1})+\varphi_2(t_{n},u_{n},w_{n})\right) =w_n+\dfrac{h}{2}\left(u_{n+1}+u_n\right). \end{cases} $$

Si on multiplie la première équation par $(u_{n+1}-u_n)$ et la deuxième par $(w_{n+1}-w_n)$ on a $$ \begin{cases} u_{n+1}^2-u_n^2=-\dfrac{h}{2}(w_{n+1}+w_n)(u_{n+1}-u_n),\\ w_{n+1}^2-w_n^2= \dfrac{h}{2}(u_{n+1}+u_n)(w_{n+1}+w_n). \end{cases} $$ En sommant les deux équations on trouve $u_{n+1}^2-u_n^2+w_{n+1}^2-w_n^2=0$ soit encore $u_{n+1}^2+w_{n+1}^2=u_n^2+w_n^2$, autrement dit $J_{n+1}=J_n$ pour tout $n$: le schéma de Crank-Nicolson preserve l'invariant.

Les fonction $\varphi_1$ et $\varphi_2$ étant linéaires, on peut resoudre le système linéaire et on obtient une écriture explicite $$ \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=\dfrac{\left(1-\frac{1}{4}h^2\right)u_n-hw_{n}}{1+\frac{1}{4}h^2},\\ w_{n+1}=\dfrac{\left(1-\frac{1}{4}h^2\right)w_n+hu_{n}}{1+\frac{1}{4}h^2}. \end{cases} $$

In [28]:
def CrankNicolson(phi1,phi2,tt):
	uu = [x0]
	ww = [y0]
	for i in range(len(tt)-1):
		uu.append(((1-0.25*h**2)*uu[i]-h*ww[i])/(1+0.25*h**2))
		ww.append(((1-0.25*h**2)*ww[i]+h*uu[i])/(1+0.25*h**2))
	return [uu,ww]
In [29]:
[uu_cn, ww_cn] = CrankNicolson(phi1,phi2,tt)
JJ_cn = [(uu_cn[i])**2+(ww_cn[i])**2 for i in range(len(tt))]

figure()
plot(tt,uu_cn,tt,ww_cn)
xlabel('t')
legend(['x(t)','y(t)'])
title('Crank-Nicolson - x(t) et y(t)') 

figure()
plot(tt,JJ_cn)
xlabel('t')
title('Crank-Nicolson - Invariant')

figure()
plot(uu_cn,ww_cn)
xlabel('x(t)')
ylabel('y(t)')
title('Crank-Nicolson - y(x)') 
Out[29]:
Text(0.5,1,'Crank-Nicolson - y(x)')

Euler symplectique $$ \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+1},u_{n+1},w_{n+1})=w_n+hu_{n+1} \end{cases} $$

In [30]:
def euler_sympl(phi1,phi2,tt):
	uu = [x0]
	ww = [y0]
	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+1],uu[i+1],ww[i]))
	return [uu,ww]

Si on multiplie la première équation par $(u_{n+1}+u_n)$ et la deuxième par $(w_{n+1}+w_n)$ on a $$ \begin{cases} u_{n+1}^2-u_n^2=-h(u_{n+1}+u_n)w_n,\\ w_{n+1}^2-w_n^2= h(w_{n+1}+w_n)u_{n+1}. \end{cases} $$ En sommant les deux équations on trouve $u_{n+1}^2-u_n^2+w_{n+1}^2-w_n^2=h(u_{n+1}w_{n+1}-w_nu_n)$, soit encore $J_{n+1}-J_n=h(u_{n+1}w_{n+1}-w_nu_n)$ : l'invariant n'est pas conservé. Cependant on voit que $$ J_{N}-J_0=h(u_Nw_N-w_0u_0)=hu_Nw_N. $$

In [31]:
[uu_es, ww_es] = euler_sympl(phi1,phi2,tt)
JJ_es = [(uu_es[i])**2+(ww_es[i])**2 for i in range(len(tt))]

figure()
plot(tt,uu_es,tt,ww_es)
xlabel('t')
legend(['x(t)','y(t)'])
title('Euler symplectique - x(t) et y(t)') 

figure()
plot(tt,JJ_es)
xlabel('t')
title('Euler symplectique - Invariant')

figure()
plot(uu_es,ww_es)
xlabel('x(t)')
ylabel('y(t)')
title('Euler symplectique - y(x)') 
Out[31]:
Text(0.5,1,'Euler symplectique - y(x)')

Eler explicite et Crank-Nicolson $$ \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+\dfrac{h}{2}\left(\varphi_2(t_{n+1},u_{n+1},w_{n+1})+\varphi_2(t_{n},u_{n},w_{n})\right) =w_n+\dfrac{h}{2}\left(u_{n+1}+u_n\right) =hu_n+\left(1-\frac{h^2}{2}\right)w_n. \end{cases} $$

\begin{align*} u_{n+1}^2+w_{n+1}^2 &= \left(u_n-hw_n\right)^2+\left(hu_n+\left(1-\frac{h^2}{2}\right)w_n\right)^2 \\ &= u_n^2+h^2w_n^2-2hu_nw_n+h^2u_n^2+\left(1-\frac{h^2}{2}\right)^2w_n^2+2h\left(1-\frac{h^2}{2}\right)u_nw_n \\ &= \left(1+h^2\right)u_n^2+\left(1+\frac{h^4}{4}\right)w_n^2-h^3u_nw_n \\ &=u_n^2+w_n^2+h^2\left(u_n^2+\frac{h^2}{4}w_n^2-hu_nw_n\right) \\ &=u_n^2+w_n^2+h^2(u_n-\frac{h}{2}w_n)^2 \end{align*} autrement dit $J_{n+1}-J_n=h^2(u_n-\frac{h}{2}w_n)^2$ : le schéma ne preserve pas l'invariant.

In [32]:
def euler_CN(phi1,phi2,tt):
	uu = [x0]
	ww = [y0]
	for i in range(len(tt)-1):
		uu.append(uu[i]+h*phi1(tt[i],uu[i],ww[i]))
		ww.append(ww[i]+0.5*h*(phi2(tt[i+1],uu[i+1],ww[i])+phi2(tt[i],uu[i],ww[i])))
	return [uu,ww]
In [33]:
[uu_ecn, ww_ecn] = euler_CN(phi1,phi2,tt)
JJ_ecn = [(uu_ecn[i])**2+(ww_ecn[i])**2 for i in range(len(tt))]

figure()
plot(tt,uu_ecn,tt,ww_ecn)
xlabel('t')
legend(['x(t)','y(t)'])
title('Euler Crank Nicolson - x(t) et y(t)') 

figure()
plot(tt,JJ_ecn)
xlabel('t')
title('Euler Crank Nicolson - Invariant')

figure()
plot(uu_ecn,ww_ecn)
xlabel('x(t)')
ylabel('y(t)')
title('Euler Crank Nicolson - y(x)') 
Out[33]:
Text(0.5,1,'Euler Crank Nicolson - y(x)')