2022 CT

• ⏱ 2H30
• ✍ Deposer votre notebook .ipynb sur Moodle (dépôt de devoir)

Dans le sujet il y a plusieurs paramètres. Pour les fixer, remplacez mon nom et prénom par les votres dans les cellules indiquées et vous obtiendrez un choix pour chaque paramètre.

Exercice 1 : Méthode d'Euler exponentielle

On considère le problème de Cauchy $$ \begin{cases} y'(t)=ay(t)+b(t), & t\in[0;2]\\ y(t_0)=y_0 \end{cases} $$ La méthode d'Euler exponentielle appliquée à ce problème de Cauchy s'écrit comme suit: $$\text{(EEx) } \begin{cases} u_0=y_0,\\ u_{n+1}=e^{ha}u_n+hg(ha)b(t_n+h/2), & n=0,\dots, N-1 \end{cases} $$ où $g(z)=\dfrac{e^z-1}{z}$.

La donnée initiale $y_0$, la constante $a$, la fonction $t\mapsto b(t)$ et le nombre de points $N$ sont donnés dans la cellule ci-dessous (modifiér nom et prénom).

%reset -f
%matplotlib inline

from matplotlib.pylab import *
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))


nom="Faccanoni"
prenom="Gloria"

print("Paramètres de l'exercice 1")
S=sum([ord(c) for c in nom+prenom])
L=list(range(1,5))
idx=S%len(L)
y0=L[idx]
prlat("y_0 = ",y0)

L=[ (-10,1), (-10,3), (-10,5), (-5,1), (-5,3) ]
idx=S%len(L)
a,d = L[idx]
c,N = -a, 1-a
prlat("a = ",a)
prlat("b(t) = ",c,"t+",d)
prlat("N = ",N)

Paramètres de l'exercice 1

$\displaystyle y_0 = 1$

$\displaystyle a = -5$

$\displaystyle b(t) = 5t+3$

$\displaystyle N = 6$

Q1 [3 points] Comparer la solution exacte, la solution approchée obtenue par la méthode d'Euler explicite (EE) et la solution obtenue avec la méthode d'Euler exponentielle (EEx) avec $N+1$ points.

# ======================================================
# variables globales
# ======================================================
t0     = 0
tfinal = 2

# ======================================================
# solution exacte symbolique
# ======================================================
t = sympy.Symbol('t')
y = sympy.Function('y')

b = lambda t :  c*t+d
phi = lambda t,y :  a*y+b(t)

edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
display(edo)
solgen = sympy.dsolve(edo,y(t))
# display(solgen)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
# display(consts)
solpar=solgen.subs(consts).simplify().factor(sympy.exp(a*t))
display(solpar.simplify())

# ======================================================
# solution exacte numérique
# ======================================================
print(f'Choix: {a=}, b(t)={c}t+{d} donc la solution exacte est')

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

# ======================================================
# SCHÉMAS
# ======================================================
schemas = ['EEx','EE']
Nb_schemas = len(schemas)

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

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

# ======================================================
# CALCUL SOLUTION APPROCHEE
# ======================================================
tt = linspace(t0,tfinal,N+1)
yy = sol_exacte(tt)

uu = { schemas[s] : eval(schemas[s])(phi,tt,y0) for s in range(Nb_schemas) }
err= { schemas[s] : norm(uu[schemas[s]]-yy,inf) for s in range(Nb_schemas) }

figure(1,figsize=(28,6))
idx=0
for key in uu:
    subplot(1,Nb_schemas,idx+1)
    plot(tt,yy,'b-',tt,uu[key],'r-D')    
    legend(['Exacte', f'{key}, {N=}'])
    title( f'{key} - max(|err|)= {err[key]:g}' )
    idx+=1

$\displaystyle \frac{d}{d t} y{\left(t \right)} = 5 t - 5 y{\left(t \right)} + 3$

$\displaystyle y{\left(t \right)} = t + \frac{2}{5} + \frac{3 e^{- 5 t}}{5}$

Choix: a=-5, b(t)=5t+3 donc la solution exacte est

Q2 [2 points] Comparer ensuite empiriquement les ordres de convergence.

# ======================================================
# CALCUL ORDRE
# ======================================================
H = []
err= { schemas[s] : [] for s in range(Nb_schemas) }

figure(2,figsize=(12,6))
ax=[]
for i in range(Nb_schemas): 
    ax.append(subplot(1,Nb_schemas,i+1))
    
for k in range(6):
    N += 20
    h = (tfinal-t0)/N
    H.append(h)
    tt = linspace(t0,tfinal,N+1)
    yy = sol_exacte(tt)
    uu = { schemas[s] : eval(schemas[s])(phi,tt,y0) for s in range(Nb_schemas) }
    idx=0
    for key in uu:
        err[key].append( norm(uu[key]-yy,inf) )
        ax[idx].plot(tt,uu[key],label=f'{key}, {N=}')    
        idx+=1

for i in range(Nb_schemas): 
    ax[i].plot(tt,yy,'b-',label='Exacte')
    ax[i].legend()
    ax[i].set_title(schemas[i])

ordre = {  key : polyfit(log(H),log(err[key]),1)[0] for key in err  }
    
figure(3,figsize=(7,7))
markers=['^', 's', 'p', 'h', '8', 'D', '>', '<', '*','+','o']

for s in range(Nb_schemas):
    loglog(H,err[schemas[s]], label=f'{schemas[s]}, ordre {ordre[schemas[s]]:1.2f}',marker=markers[s])
xlabel('$h$')
ylabel('$e$')
legend(bbox_to_anchor=(1.1, 1.05))
grid(True)

Exercice 2 : étude d'un schéma Runge-Kutta

Considérons le schéma dont la matrice de Butcher est donnée ci-dessous: $$ \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \dfrac{1}{2} & \dfrac{5}{24} & \dfrac{1}{3}&\dfrac{-1}{24}\\ 1 & 0 &1&0\\ \hline & \dfrac{1}{6} & \dfrac{2}{3} & \dfrac{1}{6} \end{array} $$

Q1 [1 points] Écrire la suite définie par réccurence associée à ce schéma.

$$\begin{cases} u_0 = y_0 \\ K_1 = \varphi\left(t_n,u_n\right)\\ K_2 = \varphi\left(t_n+\frac{h}{2},u_n+h\left(\frac{5}{24}K_1+\frac{1}{3}K_2-\frac{1}{24}K_3\right)\right)\\ K_3 = \varphi\left(t_{n+1},u_n+hK_2\right)\\ u_{n+1} = u_n + \frac{h}{6}\left(K_1+4K_2+K_3\right) & n=0,1,\dots N-1 \end{cases} $$

Q2 [3 points] Étudier théoriquement l'ordre du schéma.

Soit $\omega$ l'ordre de la méthode.

• C'est un schéma implicite à $3$ étages ($s=3$) donc $\omega\le2s=6$

• Si $\begin{cases} \displaystyle\sum_{j=1}^s b_{i}=1& \\ \displaystyle c_i=\sum_{j=1}^s a_{ij}&i=1,\dots,s \end{cases}$ alors $\omega\ge1$

• Si, de plus, $\displaystyle\sum_{j=1}^s b_j c_j=\frac{1}{2}$ alors $\omega\ge2$

• Si, de plus,
  $\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^2=\frac{1}{3}\\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j=\frac{1}{6} \end{cases}$ alors $\omega\ge3$

• Si, de plus, $\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^3=\frac{1}{4}& \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i c_i a_{ij} c_j=\frac{1}{8} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j^2=\frac{1}{12} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s\sum_{k=1}^s b_i a_{ij}a_{jk} c_k=\frac{1}{24} \end{cases}$ alors $\omega\ge4$

Calculons donc toutes ces sommes:

%reset -f
%matplotlib inline

from matplotlib.pylab import *
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))

c=[0,sympy.Rational(1,2),1]
b=[sympy.Rational(1,6),sympy.Rational(2,3),sympy.Rational(1,6)]
A=[[0,0,0],[sympy.Rational(5,24),sympy.Rational(1,3),-sympy.Rational(1,24)],[0,1,0]]
s=len(c)

print('Consistance i.e. ordre 1 ssi tous = 0 :')
prlat("\sum_{j=1}^s b_j-1 =",sum(b)-1)
for i in range(s):
    prlat(f'i={i}','\quad \sum_{j=1}^s a_{ij}-c_i=',(sum(A[i])-c[i]))

print('Ordre 2  ssi (ordre 1 et = 1/2)')
prlat('\sum_{j=1}^s b_j c_j=',sum([b[i]*c[i] for i in range(s)]))

print('Ordre 3 ssi (ordre 2 et premier = 1/3 et deuxième 1/6)')
prlat('\sum_{j=1}^s b_j c_j^2=',sum([b[i]*c[i]**2 for i in range(s)]).simplify())
prlat('\sum_{i,j=1}^s b_i a_{ij} c_j=',sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]))
print('Ordre 4 ssi (ordre 3 et premier = 1/4, deuxieme = 1/8, troisieme = 1/12 et quatrième = 1/24)')
prlat('\sum_{j=1}^s b_j c_j^3=',sum([b[i]*c[i]**3 for i in range(s)]).simplify())
prlat('\sum_{i,j=1}^s b_i c_i a_{ij} c_j=',sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]))
prlat('\sum_{i,j=1}^s b_i a_{ij} c_j^2=',sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]))
prlat('\sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k=',sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]))

Consistance i.e. ordre 1 ssi tous = 0 :

$\displaystyle \sum_{j=1}^s b_j-1 =0$

$\displaystyle i=0\quad \sum_{j=1}^s a_{ij}-c_i=0$

$\displaystyle i=1\quad \sum_{j=1}^s a_{ij}-c_i=0$

$\displaystyle i=2\quad \sum_{j=1}^s a_{ij}-c_i=0$

Ordre 2  ssi (ordre 1 et = 1/2)

$\displaystyle \sum_{j=1}^s b_j c_j=\frac{1}{2}$

Ordre 3 ssi (ordre 2 et premier = 1/3 et deuxième 1/6)

$\displaystyle \sum_{j=1}^s b_j c_j^2=\frac{1}{3}$

$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{1}{6}$

Ordre 4 ssi (ordre 3 et premier = 1/4, deuxieme = 1/8, troisieme = 1/12 et quatrième = 1/24)

$\displaystyle \sum_{j=1}^s b_j c_j^3=\frac{1}{4}$

$\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=\frac{1}{8}$

$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=\frac{5}{72}$

$\displaystyle \sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k=\frac{5}{144}$

D'après ces résultats, le schéma est d'ordre $3$.

Q3 [4 points] Étudier théoriquement la A-stabilité du schéma.

La méthode appliquée à l'EDO $y'(t)=-\beta y(t)$ s'écrit: $$\begin{cases} u_0 = y_0 \\ K_1 = -\beta\left(u_n\right)\\ K_2 = -\beta\left(u_n+\frac{5h}{24}K_1+\frac{h}{3}K_2+\frac{-h}{24}K_3\right)\\ K_3 = -\beta\left(u_n+hK_2\right)\\ u_{n+1} = u_n + \frac{h}{6}\left(K_1+4K_2+K_3\right) & n=0,1,\dots N-1 \end{cases}$$

sympy.var('u_n, h, β, K1, K2, K3')  
g = lambda y : -β*y
eq1 = sympy.Eq( K1 , g(u_n))
eq2 = sympy.Eq( K2 , g(u_n+5*h/24*K1+h/3*K2-h/24*K3) )
eq3 = sympy.Eq( K3 , g(u_n+h*K2) )
sol = sympy.solve([eq1,eq2,eq3],[K1,K2,K3])
display(sol)
RHS = u_n+h/6*(K1+4*K2+K3)
RHS = RHS.subs(sol).simplify()
prlat( 'u_{n+1}=' , RHS )
print(r"Notons $x=h\beta>0$. On trouve donc")
x = sympy.Symbol('x',positive=True)
RHS = RHS.subs(h*β,x).simplify()
prlat( 'u_{n+1}=' , RHS )

$\displaystyle \left\{ K_{1} : - u_{n} β, \ K_{2} : \frac{4 u_{n} β \left(h β - 6\right)}{h^{2} β^{2} + 8 h β + 24}, \ K_{3} : \frac{u_{n} β \left(- 5 h^{2} β^{2} + 16 h β - 24\right)}{h^{2} β^{2} + 8 h β + 24}\right\}$

$\displaystyle u_{n+1}=- \frac{u_{n} \left(h^{3} β^{3} - 5 h^{2} β^{2} + 16 h β - 24\right)}{h^{2} β^{2} + 8 h β + 24}$

Notons $x=h\beta>0$. On trouve donc

$\displaystyle u_{n+1}=- \frac{u_{n} \left(x^{3} - 5 x^{2} + 16 x - 24\right)}{x^{2} + 8 x + 24}$

Par induction on obtient $$ u_{n}=\left(-\frac{x^3-5x^2+16x-24}{x^2+8x+24}\right)^nu_0. $$ Par conséquent, $\lim\limits_{n\to+\infty}u_n=0$ si et seulement si $$ \left|\frac{x^3-5x^2+16x-24}{x^2+8x+24}\right|<1. $$

Étudions la fonction $x\mapsto\frac{x^3-5x^2+16x-24}{x^2+8x+24}$ pour $x>0$:

q = sympy.apart(-RHS/u_n)
prlat( r'-\frac{u_{n+1}}{u_n}=' , -RHS/u_n , "=" , q , "=:q(x)" )
prlat( r'q(0)=' , q.subs(x,0) )
dq=sympy.diff(q,x).simplify()
prlat( r"q'(x)=" , dq )
prlat( r"q'(x)>0 \text{ ? } " ,sympy.solve(dq>0) )
prlat( r"q(x)=1 \iff x= " ,sympy.solve(q-1) )
# prlat( r"q(x)=0 \iff x= " ,[ xi.evalf() for xi in sympy.solve(q)] )

xx=linspace(0,8,100)
yy=(xx**3-5*xx**2+16*xx-24)/(xx**2+8*xx+24)
plot(xx,yy)
plot([xx[0],xx[-1]],[-1,-1])
plot([xx[0],xx[-1]],[1,1])
grid();

$\displaystyle -\frac{u_{n+1}}{u_n}=\frac{x^{3} - 5 x^{2} + 16 x - 24}{x^{2} + 8 x + 24}=x + \frac{96 \left(x + 3\right)}{x^{2} + 8 x + 24} - 13=:q(x)$

$\displaystyle q(0)=-1$

$\displaystyle q'(x)=\frac{x^{4} + 16 x^{3} + 16 x^{2} - 192 x + 576}{x^{4} + 16 x^{3} + 112 x^{2} + 384 x + 576}$

$\displaystyle q'(x)>0 \text{ ? } \text{True}$

$\displaystyle q(x)=1 \iff x= \left[ 6\right]$

Elle est strictement croissante sur $[0,+\infty[$ et $q(0)=-1$. Il existe donc un et un seul $x_0$ tel que $q(x_0)=1$. En resolvant l'équation on trouve $x_0=6$.

La condition de stabilité est donc satisfaite ssi $0<x<6$, donc ssi $$ h<\frac{6}{\beta}. $$

Q4 [3 points] Considérons le problème de Cauchy $y'(t)=-50y(t)$ pour $t\in]0,1]$ et $y(0)=1$.
Vérifier empiriquement qu'avec $N=8$ (i.e. 9 points de discrétisation) le schéma est instable et avec $N=21$ il est stable.

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

phi = lambda t,y: -β*y

##############################################
# solution exacte
##############################################
sol_exacte = lambda t : y0*exp(-β*t)

##############################################
# solution approchée
##############################################
def RK(tt,phi,y0):
    uu = [y0]
    h = tt[1]-tt[0]
    for i in range(len(tt)-1):
        K1=phi(tt[i],uu[i])
        sys = lambda X : [-X[0]+phi(tt[i]+h/2,uu[i]+h*(5/24*K1+X[0]/3-X[1]/24)) , -X[1]+phi(tt[i+1],uu[i]+h*X[0])]
        K2, K3 = fsolve( sys , [K1,K1] )
        uu.append( uu[i]+h/6*(K1+4*K2+K3) )
    return uu

##############################################
# ordre
##############################################
# 6/β>h=(tfinal-t0)/N
# N>(tfinal-t0)*β/6

figure(figsize=(16,5))

subplot(1,2,1)
N = int((tfinal-t0)*β/6)
tt = linspace(t0, tfinal, N + 1)
h  = tt[1] - tt[0]
yy = array([sol_exacte(t) for t in tt])
uu = RK(tt, phi, y0)
err = norm(uu-yy,inf)
plot(tt,uu,label=f'Approchée avec N={N}')
plot(tt,yy,label='Exacte')
xlabel('$t$')
ylabel('$y$')
grid(True)
legend();

subplot(1,2,2)
N = int(int((tfinal-t0)*β/6)*2.46115761806205+2)
tt = linspace(t0, tfinal, N + 1)
h  = tt[1] - tt[0]
yy = array([sol_exacte(t) for t in tt])
uu = RK(tt, phi, y0)
err = norm(uu-yy,inf)
plot(tt,uu,label=f'Approchée avec N={N}')
plot(tt,yy,label='Exacte')
xlabel('$t$')
ylabel('$y$')
grid(True)
legend();

Annexe: calculons la suite géométrique pour étudier la A-stabilité dans le cas générale d'un schéma dont la matrice de Butcher est $$ \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ c_2 & a_{21} & a_{22}&a_{23}\\ 1 & 0 &1&0\\ \hline & b_1 & b_2 & b_3 \end{array} $$

sympy.var('u_n, h, β, K1, K2, K3')  
sympy.var('b_1, b_2, b_3, a_21, a_22, a_23, x')  
g = lambda y : -β*y
eq1 = sympy.Eq( K1 , g(u_n))
eq2 = sympy.Eq( K2 , g(u_n+h*a_21*K1+h*a_22*K2+h*a_23*K3) )
eq3 = sympy.Eq( K3 , g(u_n+h*K2) )
sol = sympy.solve([eq1,eq2,eq3],[K1,K2,K3])
display(sol)
RHS = u_n+h*(b_1*K1+b_2*K2+b_3*K3)
RHS = RHS.subs(sol).simplify()
prlat( 'u_{n+1}=' , RHS )
RHS = RHS.subs(h,x/β).simplify()
prlat( 'u_{n+1}=' , RHS )
print("Avec nos paramètres on obtient bien:")
RHS = RHS.subs({a_21:5/sympy.S(24),a_22:1/sympy.S(3),a_23:-1/sympy.S(24),b_1:1/sympy.S(6),b_2:4/sympy.S(6),b_3:1/sympy.S(6)}).simplify()
prlat( 'u_{n+1}=' , RHS )

$\displaystyle \left\{ K_{1} : - u_{n} β, \ K_{2} : \frac{u_{n} β \left(a_{21} h β + a_{23} h β - 1\right)}{a_{22} h β - a_{23} h^{2} β^{2} + 1}, \ K_{3} : \frac{u_{n} β \left(- a_{21} h^{2} β^{2} - a_{22} h β + h β - 1\right)}{a_{22} h β - a_{23} h^{2} β^{2} + 1}\right\}$

$\displaystyle u_{n+1}=\frac{u_{n} \left(a_{22} h β - a_{23} h^{2} β^{2} - h β \left(b_{1} \left(a_{22} h β - a_{23} h^{2} β^{2} + 1\right) - b_{2} \left(a_{21} h β + a_{23} h β - 1\right) + b_{3} \left(a_{21} h^{2} β^{2} + a_{22} h β - h β + 1\right)\right) + 1\right)}{a_{22} h β - a_{23} h^{2} β^{2} + 1}$

$\displaystyle u_{n+1}=\frac{u_{n} \left(a_{22} x - a_{23} x^{2} - x \left(b_{1} \left(a_{22} x - a_{23} x^{2} + 1\right) - b_{2} \left(a_{21} x + a_{23} x - 1\right) + b_{3} \left(a_{21} x^{2} + a_{22} x - x + 1\right)\right) + 1\right)}{a_{22} x - a_{23} x^{2} + 1}$

Avec nos paramètres on obtient bien:

$\displaystyle u_{n+1}=- \frac{u_{n} \left(x^{3} - 5 x^{2} + 16 x - 24\right)}{x^{2} + 8 x + 24}$

Exercice 3 : étude d'un système Hamiltonien

On s'intèresse au problème de Cauchy $$ \begin{cases} y'(t)=\varphi_1(t,y(t),z(t)),\\ z'(t)=\varphi_2(t,y(t),z(t)),\\ y(0)=y_0,\\ z(0)=z_0 \end{cases} $$ Soit $H\colon (y,z)\mapsto H(y,z)$ une fonction de deux variables régulière et posons $\varphi_1(t,y,z)=\dfrac{\partial H}{\partial z}$ et $\varphi_2(t,y,z)=-\dfrac{\partial H}{\partial y}$.

Q1 [1 point] Notons $\mathcal{H}\colon t \mapsto H(y(t),z(t))$. Montrer que $\mathcal{H}$ est constante.

On a $$ \mathcal{H}'(t)= \frac{d }{dt}H(y(t),z(t))= \frac{\partial H}{\partial y}y'(t)+\frac{\partial H}{\partial z}z'(t)= -z'(t)y'(t)+y'(t)z'(t)=0. $$

%reset -f
%matplotlib inline

from matplotlib.pylab import *
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))


nom="Faccanoni"
prenom="Gloria"

print("Paramètres de l'exercice 3")
L=[ (-1,0),(1,0),(0,1),(0,-1) ]
idx=sum([ord(c) for c in nom+prenom])%len(L)
y0,z0=L[idx]
prlat("y_0 = ",y0)
prlat("z_0 = ",z0)

L=[ (1,2),(2,1),(3,2),(2,3),(3,1),(1,3) ]
idx=sum([ord(c) for c in nom+prenom])%len(L)
(a,b)=L[idx]
prlat("a = ",a)
prlat("b = ",b)

Paramètres de l'exercice 3

$\displaystyle y_0 = -1$

$\displaystyle z_0 = 0$

$\displaystyle a = 3$

$\displaystyle b = 1$

Q2 [2 point] Posons $H(y,z)=\dfrac{(ay)^2+(bz)^2}{2}$ où $a$ et $b$ sont des constantes fixées dans la cellule précédente.

1. Calculer la solution exacte (les données initiales sont fixées dans la cellule précédente).
2. Dans la suite on posera $T=4\pi/(ab)$ et, pour les affichages, on chosira $N+1$ points avec $N=48$.
   Dans trois répères différents afficher
   • $t\mapsto y(t)$ et $t\mapsto z(t)$ dans le même répère,
   • $y\mapsto z(y)$,
   • $t\mapsto\mathcal{H}(t)$.

Le problème de Cauchy devient $$ \begin{cases} y'(t)=b^2z(t),\\ z'(t)=-a^2y(t),\\ y(0)=y_0,\\ z(0)=z_0 \end{cases} $$

t0 = 0

tfinal = 4*pi/(a*b)
N = 48

phi1 = lambda t,y,z : b**2*z
phi2 = lambda t,y,z : -a**2*y
Ham = lambda y,z : ((a*y)**2+(b*z)**2)/2

La solution exacte est

_h = sympy.Symbol('h',positive=True)

_phi1 = lambda t,y,z :  sympy.diff(Ham(y,z),z)
_phi2 = lambda t,y,z : -sympy.diff(Ham(y,z),y)

t = sympy.Symbol('t')
y = sympy.Function('y')
z = sympy.Function('z')

edo1 = sympy.Eq( sympy.diff(y(t),t) , _phi1(t,y(t),z(t)) )
edo2 = sympy.Eq( sympy.diff(z(t),t) , _phi2(t,y(t),z(t)) )
# display(edo1)
# display(edo2)
solgen = sympy.dsolve([edo1,edo2],[y(t),z(t)])
# display(solgen)
consts = sympy.solve( [ sympy.Eq( y0, solgen[0].rhs.subs(t,t0)) , sympy.Eq( z0, solgen[1].rhs.subs(t,t0)) ] , dict=True)[0]
# display(consts)
solpar_1=solgen[0].subs(consts)
solpar_2=solgen[1].subs(consts)
display(solpar_1,solpar_2)

$\displaystyle y{\left(t \right)} = - \cos{\left(3 t \right)}$

$\displaystyle z{\left(t \right)} = 3 \sin{\left(3 t \right)}$

Les trois graphes sont

func_1 = sympy.lambdify(t,solpar_1.rhs,'numpy')
func_2 = sympy.lambdify(t,solpar_2.rhs,'numpy')

tt = linspace(t0,tfinal,N+1)
h  = tt[1]-tt[0]

figure(figsize=(18,5))
yy = func_1(tt)
zz = func_2(tt)
HH = [ Ham(y,z) for (y,z) in zip(yy,zz) ]

subplot(1,3,1)
plot(tt,yy,'-*',tt,zz,'-*')
legend([r'$t\mapsto y(t)$',r'$t\mapsto z(t)$'])
xlabel(r'$t$')

subplot(1,3,2)
plot(yy,zz,'-*')
xlabel(r'$y$')
ylabel(r'$z$')
axis('equal')

subplot(1,3,3)
plot(tt,HH,'-*')
xlabel(r'$t$')
ylabel(r'$\mathcal{H}$')
axis('equal');

Q3 [4 points]
Dans la suite on notera $u_n\approx y(t_n)$, $w_n\approx z(t_n)$ et $J_n\approx\mathcal{H}(t_n)$.

Calculer la solution approchée obtenue par la méthode suivante appelée méthode d'Euler symplectique.
$$\text{(ES) }\begin{cases} u_0=y_0,\\ w_0=z_0,\\ u_{n+1}=u_n+h\varphi_1\left(t_n,u_n,w_n\right),\\ w_{n+1}=w_n+h\varphi_2\left(t_{n+1},u_{n+1},w_{n+1}\right). \end{cases}$$ Nota bene: le choix particulier de $\varphi_1$ et $\varphi_2$ permet de réécrire le schéma sous forme explicite.

Dans trois répère différents afficher, en comparant solution exacte et solution approchée,

• $\{(t_i,u_i)\}$ et $\{(t_i,w_i)\}$ dans un même répère,
• $\{(u_i,w_i)\}$,
• $\{(t_i,J_i)\}$ où $J_i=H(u_i,w_i)$.

Pour l'affichage des solutions exactes VS approchées, on peut écrire une fonction qui prend en paramètres la discrétisation tt, les deux listes solution exacte yy et zz, les deux listes solution approchée uu et ww, puis une chaîne de caractères s qui contient le nom du schéma, et pour finir HH la liste qui contient le calcul de l'Hamiltonien exact et HH_num celui approché.

def affichage(tt,yy,zz,uu,ww,s,HH,HH_num):

    subplot(1,3,1)
    plot(tt,uu,'o--',tt,ww,'d--')
    plot(tt,yy,tt,zz)
    xlabel('t')
    legend([r'$u(t)$ approchée',r'$w(t)$ approchée','$y(t)$ exacte','$z(t)$ exacte'])
    title(f'{s} - y(t) et z(t)') 
    grid()

    subplot(1,3,2)
    plot(uu,ww,'o--')
    plot(yy,zz)
    xlabel('y')
    ylabel('z')
    legend(['Approchée','Exacte'])
    title(f'{s} - z(y)')
    grid()
    axis('equal')

    subplot(1,3,3)
    plot(tt,HH_num,'o--')
    plot(tt,HH)
    xlabel('t')
    ylabel(r'$\mathcal{H}$')
    legend(['Approchée','Exacte'])
    title(f'{s} - '+r'$\mathcal{H}(t)$')
    grid()
    axis('equal');

Le schéma se réécrit $$\text{(ES) }\begin{cases} u_0=y_0,\\ w_0=z_0,\\ u_{n+1}=u_n+h b^2 w_n,\\ w_{n+1}=-a^2h u_n + (1-(abh)^2)w_n. \end{cases}$$ et $$ J_{n+1}= \frac{(au_{n+1})^2+(bw_{n+1})^2}{2} $$ Vérifions nos calculs:

_h=sympy.Symbol('h')
t_n=sympy.Symbol('t_n')
u_n=sympy.Symbol('u_n')
w_n=sympy.Symbol('w_n')
t_np1=t_n+h
u_np1=sympy.Symbol('u_{n+1}')
w_np1=sympy.Symbol('w_{n+1}')
eq1 = sympy.Eq( u_np1 , u_n+_h*_phi1(t_n,u_n,w_n))
eq2 = sympy.Eq( w_np1 , w_n+_h*_phi2(t_np1,u_np1,w_np1))
sol = sympy.solve([eq1,eq2],[u_np1,w_np1])
prlat( 'u_{n+1}=' , sol[u_np1] )
prlat( 'w_{n+1}=' , sol[w_np1].factor(w_n) )
# prlat( 'H_{n+1}=' , _Ham(sol[u_np1],sol[w_np1]).expand().simplify() )

$\displaystyle u_{n+1}=h w_{n} + u_{n}$

$\displaystyle w_{n+1}=- 9 h u_{n} - w_{n} \left(9 h^{2} - 1\right)$

# x0, y0, h variables globales
def ES(phi1,phi2,tt):
    uu = [y0]
    ww = [z0]
    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],0)) 
        # idem que
        uu.append( uu[i]+h*b**2*ww[i]  )
        ww.append( -a**2*h*uu[i]+(1-(a*b*h)**2)*ww[i]  )
    return [uu,ww]

[uu, ww] = ES(phi1,phi2,tt)
HH_num = [ Ham(u,w) for (u,w) in zip(uu,ww) ]

figure(figsize=(18,5))
affichage(tt,yy,zz,uu,ww,"Euler symplectique",HH,HH_num)