None
⏱ 2h30, deposer le notebook dans Moodle
Q [4 points] Pour calculer des solutions approchées on doit utiliser un schéma numérique. Quelle caractéristique doit-il avoir ? Justifier la réponse.
🐸 Il s'agit d'un des problèmes de Cauchy proposés lors du CT 2023 de l'ECUE M61 Équations différentielles.
Bien qu'on ne puisse pas calculer analytiquement la solution exacte, nous pouvons tracer les lignes de courant pour avoir une idée du comportement de la solution en fonction de la donnée initiale. On en déduit que le problème de Cauchy est numériquement mal posé. Il faut alors utiliser un schéma d'ordre elevé.
Pour s'en convaincre nous pouvons aussi demander à odeint
du module scipy
de tracer la solution approchée pour différentes valeurs de $y_0$.
%reset -f
%matplotlib inline
from matplotlib.pylab import *
from scipy.integrate import odeint
# variables globales
t0 = 0
tfinal = 10
phi = lambda t,y: sin(pi*2*y)*(1-1/(4+t)-y)
N = 200
tt = linspace(t0, tfinal, N + 1)
figure(figsize=(10,5))
# On trace d'abord les lignes de courant
g1 = linspace(0,10,21)
g2 = linspace(0.5,1,21)
T,Y = meshgrid(g1,g2)
DT, DY = 1, phi(T,Y) # compute growth rate on the grid
M = hypot(DT,DY) # = sqrt(DT**2+DY**2) # norm growth rate
streamplot(T,Y, DT/M, DY/M, color=M, linewidth=0.5, cmap=cm.viridis, density=2, arrowstyle='->', arrowsize=1.5)
grid()
xlabel('t')
ylabel('y')
title('Lignes de courant');
# On trace quelque solution approchée
for y0 in linspace(0.7,0.9,7):
uu = odeint(phi,y0,tt, tfirst=True)
plot(tt,uu,lw=4,label=f'$y_0={y0:g}$')
legend(bbox_to_anchor=(1.04, 1), loc='upper left');
On remarque qu'il s'agit de la méthode de Crank-Nicolson sur l'intervalle $[t_{n-1},t_{n+1}]$.
Q2 [1 points] Étudier la zéro-stabilité.
Le schéma est de la forme $$ u_{n+1} = a_0u_n+a_1u_{n-1}+hb_0\varphi_n+hb_1\varphi_{n-1}+hb_{-1}\varphi_{n+1}. $$ On a
Le premier polynôme caractéristique est $$ \varrho(r)= %r^{p+1}-\sum_{j=0}^p a_jr^{p-j}= r^2-a_0r-a_1= r^2-1 =(r-1)(r+1). $$
Les racines sont $$ r_0=1,\quad r_{1}=-1. $$ La méthode est donc zéro-stable car $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \text{ multiplicité de $r_j$ > }1 \implies |r_j|<1 \end{cases} $$
Q2 [3 points] La méthode est elle consistente? Est elle convergente? Si oui, étudier théoriquement l'ordre de convergence.
La méthode étant zéro-stable, si elle est consistante, alors elle est convergente. Pour étudier la convergence introduisons la quantité $$ \begin{align*} \xi(i) &= \sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j \\ &= (-0)^{i}a_0+(-1)^ia_1+i \Big( (1)^{i-1}b_{-1}+(-0)^{i-1}b_0+(-1)^{i-1}b_1 \Big) \end{align*} $$ avec $(-0)^0=1$.
Pour que la méthode soit consistante il faut que $\xi(0)=\xi(1)=1$. De plus, pour que la méthode soit d'ordre $\omega$, il faut que $\xi(i)=1$ pour tout $i\le\omega$.
De plus, $\omega$ vérifie la première barrière de Dahlquist: le schéma étant implicite à $q=2$ (pair) pas consistante et zéro-stable, $\omega=q+2\le4$.
En calculant $\xi$ on trouve que la méthode est d'ordre 2 (cf. calculs avec sympy
ci-dessous).
%reset -f
import sympy
aa_eval = [0,1]
bb_eval = [0,1]
bm1_eval = 1
###########################################################################
from IPython.display import display, Math
def xi(i,aa,bb,bm1):
sa = sum( [ (-j)**i * aa[j] for j in range(q) ] )
sb = bm1+sum( [(-j)**(i-1)*bb[j] for j in range(1,q)] )
if i==1:
sb += bb[0]
# return (sa).factor()+(i*sb).factor()
return (sa)+(i*sb)
q = len(aa_eval) # nb de pas
p = q-1 # j = 0...p # coeffs du schema
omega = q+2 if q%2==0 else q+1
aa = sympy.symbols(f'a_0:{q}')
bb = sympy.symbols(f'b_0:{q}')
bm1 = sympy.Symbol('b_{-1}')
expimp = "explicite" if bm1==0 else "implicite"
print(f"{'★'*70}\nC'est une méthode {expimp} à q = {q} pas d'ordre ω ≤ {omega} avec")
texte = ""
for j in range(p+1):
texte +=f"{sympy.latex(aa[j])}={sympy.latex(aa_eval[j])}\qquad "
texte += f"{sympy.latex(bm1)}={sympy.latex(bm1_eval)}\qquad "
for j in range(p+1):
texte += f"{sympy.latex(bb[j])}={sympy.latex(bb_eval[j])}\qquad "
display(Math( texte ))
print(f"{'★'*70}")
xi_eval = [xi(i,aa_eval,bb_eval,bm1_eval) for i in range(omega+1)]
for i in range(omega+1):
display(Math(f"ξ({i}) = {sympy.latex(xi(i,aa,bb,bm1))} = {sympy.latex(xi_eval[i])}") )
if xi_eval[i]!=1:
print(f"La méthode est d'ordre ω = {i-1}\n{'★'*70}")
break
Q3 [3 points] Vérifier empiriquement l'ordre de convergence sur le problème de Cauchy $$\begin{cases} y'(t) = (y(t)+2)\dfrac{2t^2+2t-1}{1+t}, &\forall t \in I=[0,2],\\ y(0) = -1, \end{cases} $$
%matplotlib inline
from matplotlib.pylab import *
rcParams.update({'font.size': 12})
from scipy.optimize import fsolve
t0 = 0
tfinal = 2
y0 = -1
phi = lambda t,y: (y+2)*(2*t**2+2*t-1)/(1+t)
##############################################
# solution exacte
##############################################
t = sympy.Symbol('t')
y = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
solgen = sympy.dsolve(edo)
# display(solgen)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
# display(consts)
solpar = solgen.subs(consts)
display(solpar)
sol_exacte = sympy.lambdify(t,solpar.rhs,'numpy')
##############################################
# schéma (initialisation avec sol exacte pour ordre de convergence)
##############################################
def multipas(phi, tt, sol_exacte):
h = tt[1] - tt[0]
uu = [sol_exacte(tt[i]) for i in range(2)]
for i in range(1,len(tt) - 1):
eq = lambda x : -x+uu[i-1]+h*(phi(tt[i+1],x)+phi(tt[i-1],uu[i-1]))
temp = fsolve( eq ,uu[i])
uu.append( temp[0] )
return uu
##############################################
# ordre
##############################################
H = []
err = []
N = 10
figure(figsize=(16,5))
ax1 = subplot(1,2,1)
for k in range(6):
N += 20
tt = linspace(t0, tfinal, N + 1)
H.append( tt[1] - tt[0] )
yy = sol_exacte(tt)
uu = multipas(phi, tt, sol_exacte)
err.append( norm(uu-yy,inf) )
ax1.plot(tt,uu,label=f'Approchée avec N={N}')
ax1.plot(tt,yy,label='Exacte')
xlabel('$t$')
ylabel('$y$')
ax1.grid(True)
ax1.legend()
ax2 = subplot(1,2,2)
ax2.plot( log(H), log(err), 'r-o')
xlabel(r'$\ln(h)$')
ylabel(r'$\ln(e)$')
title(f'Multipas ordre = {polyfit(log(H),log(err),1)[0]:1.2f}')
ax2.grid(True)
%reset -f
%matplotlib inline
from matplotlib.pylab import *
from scipy.optimize import fsolve
import sympy
sympy.init_printing()
Q1 [5 points] Pour quelles valeurs des paramètres le schéma est d'ordre 3? Et d'ordre 4?
Soit $\omega$ l'ordre de la méthode.
C'est un schéma implicite à $2$ étages ($s=2$) donc $\omega\le2s=4$.
sigma = sympy.Symbol(r"\sigma")
theta = sympy.Symbol(r"\vartheta")
gamma = sympy.Symbol(r"\gamma")
from IPython.display import display, Math
prlat = lambda *args : display(Math(''.join( sympy.latex(arg) for arg in args )))
c = [sigma,1]
b = [gamma,1-gamma]
A = [[sigma,0],[theta,1-theta]]
s = len(c)
# Affichage matrice de Butcher
print(f'{"Matrice de Butcher":_^80}')
But = sympy.Matrix(A)
But=But.col_insert(0,sympy.Matrix(c))
last=[0]
last.extend(b)
But=But.row_insert(s,sympy.Matrix(last).T)
prlat(sympy.Matrix(But))
#################
print(f'{r"Conditions pour être d ordre 1" :_^80}')
texte = "\sum_{j=1}^s b_j =" + f"{sum(b)}"
texte += r"\text{ doit être égale à }1"
display(Math(texte))
for i in range(s):
texte = f'i={i},\quad' +r'\sum_{j=1}^s a_{ij}=' + sympy.latex(sum(A[i]))
texte += r"\text{ doit être égale à }c_i="+sympy.latex(c[i])
display(Math( texte ))
print(f'{r"Conditions à ajouter pour être d ordre 2" :_^80}')
texte = '\sum_{j=1}^s b_j c_j=' +sympy.latex(sum([b[i]*c[i] for i in range(s)]))
texte += r"\text{ doit être égale à }\frac{1}{2}"
display(Math(texte))
print(f'{r"Conditions à ajouter pour être d ordre 3" :_^80}')
texte = '\sum_{j=1}^s b_j c_j^2='
texte += sympy.latex( sum([b[i]*c[i]**2 for i in range(s)]).simplify() )
texte += r"\text{ doit être égale à }\frac{1}{3}"
display(Math(texte))
texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j='
texte = texte + sympy.latex(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]))
texte += r"\text{ doit être égale à }\frac{1}{6}"
display(Math(texte))
On commence à resoudre le système dont les inconnues sont les paramètres $\sigma, \gamma, \vartheta$ pour avoir au moins l'ordre 3:
cnst = sympy.solve( [ sum([b[i]*c[i] for i in range(s)])-1/sympy.S(2),
sum([b[i]*c[i]**2 for i in range(s)])-1/sympy.S(3),
sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])-1/sympy.S(6)
] )[0]
display(cnst)
print(f'{"La matrice de Butcher devient":_^80}')
But = But.subs(cnst)
prlat(sympy.Matrix(But))
Le schéma n'est pas d'ordre 4 car:
print(f'{r"Conditions à ajouter pour être d ordre 4" :_^80}')
texte = '\sum_{j=1}^s b_j c_j^3='
texte = texte + sympy.latex( sum([b[i]*c[i]**3 for i in range(s)]).subs(cnst) )
texte += r"\text{ doit être égale à }\frac{1}{4}"
display(Math(texte))
texte = '\sum_{i,j=1}^s b_i c_i a_{ij} c_j='
texte = texte + sympy.latex( sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).subs(cnst) )
texte += r"\text{ doit être égale à }\frac{1}{8}"
display(Math(texte))
texte = '\sum_{i,j=1}^s b_i a_{ij} c_j^2='
texte = texte + sympy.latex( sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).subs(cnst) )
texte += r"\text{ doit être égale à }\frac{1}{12}"
display(Math(texte))
texte = '\sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k='
texte = texte + sympy.latex( 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)]).subs(cnst) )
texte += r"\text{ doit être égale à }\frac{1}{24}"
display(Math(texte))
Q2 [3 points] Étudier empiriquement l'ordre de convergence sur le problème de Cauchy $$ \begin{cases} y'(t)=\dfrac{t^2+y^2(t)}{ty(t)}, & t\in[1;5]\\ y(1)=1 \end{cases} $$
# variables globales
t0 = 1
tfinal = 5
y0 = 1
phi = lambda t,y: (t**2+y**2)/(t*y)
##############################################
# solution exacte
##############################################
import sympy
sympy.init_printing()
t = sympy.Symbol('t')
y = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
solgenLIST = sympy.dsolve(edo,y(t))
display(solgenLIST)
for solgen in solgenLIST:
#display(solgen)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)
#display(consts)
if consts!=[]:
solpar = solgen.subs(consts[0]).simplify()
display(solgen)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)
solpar = solgen.subs(consts).simplify()
display(solpar)
##############################################
# solution approchée
##############################################
sol_exacte = sympy.lambdify(t,solpar.rhs,'numpy')
# Je note sigma_e etc pour sigma evaluée comme un float
sigma_e = 1/3
gamma_e = 1/(2*(1-sigma_e)) # pour ordre 2
theta_e = 1
def RK(tt,phi,y0):
uu = [y0]
h = tt[1]-tt[0]
for i in range(len(tt)-1):
sys = lambda X : [ -X[0]+phi( tt[i]+h*sigma_e , uu[i]+h*X[0]*sigma_e ) ,
-X[1]+phi( tt[i]+h*theta_e , uu[i]+h*(X[0]*theta_e+X[1]*(1-theta_e)) )
]
K_start = phi(tt[i],uu[i])
K1, K2 = fsolve( sys , [K_start,K_start ] )
uu.append( uu[i]+h*(gamma_e*K1+(1-gamma_e)*K2) )
return uu
##############################################
# ordre
##############################################
H = []
err = []
figure(figsize=(16,5))
ax1 = subplot(1,2,1)
ax2 = subplot(1,2,2)
N = 50
for k in range(6):
N += 50
tt = linspace(t0, tfinal, N + 1)
H.append( tt[1] - tt[0] )
yy = sol_exacte(tt)
uu = RK(tt,phi,y0)
err.append( norm(uu-yy,inf) )
ax1.plot(tt,uu,label=f'Approchée avec N={N}')
ax1.plot(tt,yy,label='Exacte')
xlabel('$t$')
ylabel('$y$')
ax1.grid(True)
ax1.legend()
ax2.plot( log(H), log(err), 'r-o')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
title(f'RK : ordre = {polyfit(log(H),log(err),1)[0]:1.2f}')
ax2.grid(True)
Q3 [2 points] Étudier théoriquement la A-stabilité du schéma.
La méthode appliquée à l'EDO $y'(t)=-\beta y(t)$ avec $\beta>0$ s'écrit: $$\begin{cases} u_0 = y_0 \\ K_1 = -\beta u_n-(\beta h) \sigma K_1 \leadsto K_1=\frac{-\beta}{1+(\beta h) \sigma}u_n\\ K_2 = -\beta\left(u_n+h K_1 \vartheta + h K_2 (1-\vartheta)\right)\leadsto K_2=\frac{1+\frac{(\beta h)\vartheta}{1+(\beta h) \sigma}}{1+(\beta h)(1-\vartheta)}(-\beta u_n)\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_1+K_2\right) & n=0,1,\dots N-1 \end{cases} $$
En particulier, dans le cas d'ordre 3, on a $$\begin{cases} u_0 = y_0 \\ K_1 = -\beta u_n-\beta \frac{h}{3}K_1 \leadsto K_1=\frac{-3\beta}{3+\beta h}u_n\\ K_2 = -\beta\left(u_n+hK_1\right)\leadsto K_2=\frac{2\beta^2h-3\beta}{3+\beta h}u_n\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_1+K_2\right) & n=0,1,\dots N-1 \end{cases} $$ L'étude ci-dessous montre que, dans ce cas, le schéma est A-stable ssi $h<\frac{6}{\beta}$.
sympy.var('u_n, h, β, K1, K2,')
g = lambda y : -β*y
eq1 = sympy.Eq( K1 , g(u_n+h*K1*sigma))
eq2 = sympy.Eq( K2 , g(u_n+h*(K1*theta+K2*(1-theta))) )
sol = sympy.solve([eq1,eq2],[K1,K2])
sol = { key:value.factor() for key,value in sol.items() }
display(sol)
RHS = u_n+h*(gamma*K1+(1-gamma)*K2).factor()
RHS = RHS.subs(sol).simplify().factor()
texte = 'u_{n+1}=' + sympy.latex(RHS)
display(Math(texte))
print(r"Notons x=hβ>0. On trouve donc")
x = sympy.Symbol('x',positive=True)
RHS = RHS.subs(h*β,x).simplify()
texte = 'u_{n+1}=' + sympy.latex(RHS)
display(Math(texte))
print("Considérons le schéma d'ordre 3:")
display(cnst)
RHS = RHS.subs(cnst)
texte = r"\text{La méthode est A-stable ssi }|q(x)|<1\text{ avec }q(x)=\frac{u_{n+1}}{u_n}="
texte += sympy.latex(RHS/u_n)
display(Math(texte))
sympy.solve([-1<RHS/u_n,RHS/u_n<1])