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.
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
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
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)
Q1 [1 points] Écrire la suite définie par réccurence associée à ce schéma.
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$
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 :
Ordre 2 ssi (ordre 1 et = 1/2)
Ordre 3 ssi (ordre 2 et premier = 1/3 et deuxième 1/6)
Ordre 4 ssi (ordre 3 et premier = 1/4, deuxieme = 1/8, troisieme = 1/12 et quatrième = 1/24)
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 )
Notons $x=h\beta>0$. On trouve donc
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();
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 )
Avec nos paramètres on obtient bien:
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
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.
- Calculer la solution exacte (les données initiales sont fixées dans la cellule précédente).
- 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)
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() )
# 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)