Code
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from scipy.optimize import fsolve
from scipy.integrate import solve_ivp, odeint
from IPython.display import display, Markdown, Math2023 Contrôle Terminal session 1
⏱ 2h30, deposer le notebook dans Moodle
EXERCICE 1 : étude d’un système
On s’intèresse au problème de Cauchy \( \begin{cases} y'(t)=\sin(y(t))-2z(t)e^{y(t)},\\ z'(t)=z(t)\Big(z(t)e^{y(t)}-\cos(y(t))\Big),\\ y(0)=\frac{1}{2},\\ z(0)=\frac{1}{4}, \end{cases} \qquad t\in[0,30] \)
Q1 [1 point] Notons \(E(y,z) = z\sin(y)-z^2e^y\) et \(\mathcal{E}(t)=E(y(t),z(t))\). Montrer que \(\mathcal{E}\) est un invariant.
On remarque que \(\frac{\partial E}{\partial y}(y(t),z(t))=-z'(t)\) et \(\frac{\partial E}{\partial z}(y(t),z(t))=y'(t)\). On a alors \( \mathcal{E}'(t)= \frac{d }{dt}E(y(t),z(t))= \frac{\partial E}{\partial y}y'(t)+\frac{\partial E}{\partial z}z'(t)= -z'(t)y'(t)+y'(t)z'(t)=0. \) Donc \(\mathcal{E}(t)=\mathcal{E}(0)=0\).
Q2 [2 points] Calculer la solution approchée obtenue avec le schéma d’Euler Explicite avec \(N=100\). Dans trois repères différents afficher
- dans le premier repère \(t\mapsto y(t)\) et \(t\mapsto z(t)\) dans le même,
- dans le deuxième repère \(y\mapsto z(y)\),
- dans le troisième repère \(t\mapsto\mathcal{E}(t)\).
Code
def affichage(tt,uu,ww,HH_num,s):
plt.suptitle(f'Schéma {s}')
#plt.tight_layout()
plt.subplot(1,3,1)
plt.plot(tt,uu,'o--',label=r'$y(t)$')
plt.plot(tt,ww,'d--',label=r'$z(t)$')
plt.xlabel('t')
plt.legend()
plt.grid()
plt.subplot(1,3,2)
plt.plot(uu,ww,'o--')
plt.xlabel('y')
plt.ylabel('z')
plt.grid()
Y1,Y2 = np.meshgrid(np.linspace(min(uu),max(uu),21),np.linspace(min(ww),max(ww),21))
V1,V2 = pphi(tt,[Y1,Y2])
r1 = np.sqrt(1+V1**2)
r2 = np.sqrt(1+V2**2)
plt.quiver(Y1, Y2, V1/r1, V2/r2, cmap=plt.cm.viridis, scale_units='xy')
# plt.axis('equal')
plt.subplot(1,3,3)
plt.plot(tt,HH_num,'o--')
plt.xlabel('t')
plt.ylabel(r'$\mathcal{E}$')
plt.grid();
t0 = 0
y0, z0 = 0.5, 0.25
tfinal = 30
N = 100
phi1 = lambda t,y,z : np.sin(y)-2*z*np.exp(y)
phi2 = lambda t,y,z : z*(z*np.exp(y)-np.cos(y))
Ham = lambda y,z : z*np.sin(y)-z**2*np.exp(y)
tt = np.linspace(t0,tfinal,N+1)
# pour odeint
pphi = lambda t,yy : [phi1(t,yy[0],yy[1]),phi2(t,yy[0],yy[1])]
yy0 = [y0,z0]Code
def EE(phi1,phi2,tt,y0,z0):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
ww = np.zeros_like(tt)
uu[0] = y0
ww[0] = z0
for i in range(len(tt)-1):
uu[i+1] = uu[i]+h*phi1(tt[i],uu[i],ww[i])
ww[i+1] = ww[i]+h*phi2(tt[i],uu[i],ww[i])
return [uu,ww]
[uu, ww] = EE(phi1,phi2,tt,y0,z0)
HH_num = Ham(uu,ww)
plt.figure(figsize=(18,5))
affichage(tt,uu,ww,HH_num,"EE")
Q3 [3 points]
Même exercice avec le schéma de Cranck Nicolson.
Code
def CN(phi1,phi2,tt,y0,z0):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
ww = np.zeros_like(tt)
uu[0] = y0
ww[0] = z0
for i in range(len(tt)-1):
sys = lambda vv : [-vv[0]+uu[i]+h/2*(phi1(tt[i],uu[i],ww[i])+phi1(tt[i+1], vv[0], vv[1])),
-vv[1]+ww[i]+h/2*(phi2(tt[i],uu[i],ww[i])+phi2(tt[i+1], vv[0], vv[1])) ]
uu[i+1], ww[i+1] = fsolve( sys, [uu[i],ww[i]] )
return [uu,ww]
[uu, ww] = CN(phi1,phi2,tt,y0,z0)
HH_num = Ham(uu,ww)
plt.figure(figsize=(18,5))
affichage(tt,uu,ww,HH_num,"CN")
Q4 [1 points]
Même exercice avecodeint.
Code
uu,ww = odeint(pphi,yy0,tt,tfirst=True).T
HH_num = Ham(uu,ww)
plt.figure(figsize=(18,5))
affichage(tt,uu,ww,HH_num,"scipy.odeint")
EXERCICE 2 : étude d’un schéma multipas
On notera \(\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)\). Considérons la méthode multipas \( u_{n+1}=\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+h\frac{2}{3}\varphi_{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}+h\left(b_0\varphi_n+b_1\varphi_{n-1}\right)+hb_{-1}\varphi_{n+1}. \) On a
- \(p=1\): c’est un méthode à \(q=p+1=2\) pas;
- \(a_0=\frac{4}{3}\) et \(a_1=-\frac{1}{3}\),
- \(b_0=0\) et \(b_1=0\),
- \(b_{-1}=\frac{2}{3}\neq0\) (\(\implies\) la méthode est implicite).
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-\frac{4}{3}r+\frac{1}{3} =(r-1)\left(r-\frac{1}{3}\right). \)
Les racines sont \( r_0=1,\quad r_{1}=\frac{1}{3}. \) 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.
On a prouvé que la méthode est zéro-stable. Ainsi, si elle est consistante, alors elle est convergente.
Introduisons la quantité \( \begin{aligned} \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) \\ &= (-0)^{i}\frac{4}{3}-(-1)^i\frac{1}{3}+ (1)^{i-1}i \frac{2}{3} \end{aligned} \) 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:
\(\xi(0) = (-0)^{0}\frac{4}{3}-(-1)^0\frac{1}{3}+ (1)^{0-1}0 \frac{2}{3} = \frac{4}{3}-\frac{1}{3} = 1\)
\(\xi(1) = (-0)^{1}\frac{4}{3}-(-1)^1\frac{1}{3}+ (1)^{1-1}1 \frac{2}{3} = \frac{1}{3}+ \frac{2}{3}=1\)
\(\xi(2) = (-0)^{2}\frac{4}{3}-(-1)^2\frac{1}{3}+ (1)^{2-1}2 \frac{2}{3} = -\frac{1}{3}+ 2 \frac{2}{3} = 1\)
\(\xi(3)=\frac{1}{3}+ 3 \frac{2}{3}=\frac{7}{3}\)
(vérification avec sympy ci-dessous).
Code
aa_eval = [sp.Rational(4,3), -sp.Rational(1,3)]
bb_eval = [0,0]
bm1_eval = sp.Rational(2,3)
#######################################################
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()
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 = sp.symbols(f'a_0:{q}')
bb = sp.symbols(f'b_0:{q}')
bm1 = sp.Symbol('b_{-1}')
display(Markdown(f"C'est une méthode à $q={p + 1}$ pas d'ordre $\\omega\\le{omega}$ avec"))
texte = ''.join([ f"{sp.latex(aa[j])}={sp.latex(aa_eval[j])},"+r"\quad " for j in range(p+1)])
texte += ''.join([ f"{sp.latex(bb[j])}={sp.latex(bb_eval[j])},"+r"\quad " for j in range(p+1)])
texte += f"{sp.latex(bm1)}={sp.latex(bm1_eval)}"
display(Math( texte ))
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"\\xi({i}) = {sp.latex(xi(i,aa,bb,bm1))} = {sp.latex(xi_eval[i])}" ))
if xi_eval[i]!=1:
display(Markdown( f"La méthode est donc d'ordre $\\omega = {i-1}$" ))
break
C’est une méthode à \(q=2\) pas d’ordre \(\omega\le4\) avec
\(\displaystyle a_{0}=\frac{4}{3},\quad a_{1}=- \frac{1}{3},\quad b_{0}=0,\quad b_{1}=0,\quad b_{-1}=\frac{2}{3}\)
\(\displaystyle \xi(0) = a_{0} + a_{1} = 1\)
\(\displaystyle \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = 1\)
\(\displaystyle \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) = 1\)
\(\displaystyle \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) = \frac{7}{3}\)
La méthode est donc d’ordre \(\omega = 2\)
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^2}, &\forall t \in I=[2,5],\\ y(2) = 1, \end{cases} \)
Code
# variables globales
t0 = 2
tfinal = 5
y0 = 1
phi = lambda t,y: (y+2)*(2*t**2+2*t-1)/(1-t**2)
##############################################
# solution exacte
##############################################
t = sp.Symbol('t')
y = sp.Function('y')
edo = sp.Eq( sp.diff(y(t),t) , phi(t,y(t)) )
solgen = sp.dsolve(edo)
# display(solgen)
consts = sp.solve( sp.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
# display(consts)
solpar = solgen.subs(consts)
display(solpar)
sol_exacte = sp.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 = np.zeros_like(tt)
uu[0:2] = sol_exacte(tt[0:2])
# uu = [sol_exacte(tt[i]) for i in range(2)]
for i in range(1,len(tt) - 1):
eq = lambda x : -x+4/3*uu[i]-1/3*uu[i-1]+h*2/3*phi(tt[i+1],x)
temp = fsolve( eq ,uu[i])
uu[i+1] = temp[0]
return uu
##############################################
# ordre
##############################################
H = []
err = []
N = 10
plt.figure(figsize=(16,5))
ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
for k in range(6):
N += 20
tt = np.linspace(t0, tfinal, N + 1)
H.append( tt[1] - tt[0] )
yy = sol_exacte(tt)
uu = multipas(phi, tt, sol_exacte)
err.append( np.linalg.norm(uu-yy,np.inf) )
ax1.plot(tt,uu,label=f'Approchée avec N={N}')
ax1.plot(tt,yy,label='Exacte') # sur la grille la plus fine
ax1.grid(True)
ax1.legend()
ax2 = plt.subplot(1,2,2)
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Multipas ordre = {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)\(\displaystyle y{\left(t \right)} = -2 + \frac{3 \sqrt{3} e^{4} e^{- 2 t}}{\left(t - 1\right)^{\frac{3}{2}} \sqrt{t + 1}}\)
EXERCICE 3 : étude d’un schéma RK
Considérons le schéma dont la matrice de Butcher est donnée ci-dessous: \( \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \eta & \eta & 0 & 0\\ \sigma & 0 & \sigma & 0\\ \hline & 1-\gamma & 0 & \gamma \end{array} \)
Dans le sujet de cet exercice, il y a un paramètre \(\eta\). Pour le fixer, remplacez mon nom et prénom par les vôtres et vous obtiendrez la valeur pour ce paramètre dans tout l’exercice.
Code
# =======================================================
nom = "Faccanoni"
prenom = "Gloria"
# =======================================================
L = list(range(2,5))
idx = sum([ord(c) for c in nom+prenom])%len(L)
print("Paramètre de l'exercice ")
display(Math(f"\\eta = {sp.latex(sp.Rational(1,L[idx]))}"))Paramètre de l'exercice \(\displaystyle \eta = \frac{1}{3}\)
Q1 [3 points] Pour quelles valeurs des paramètres \(\sigma\) et \(\gamma\) le schéma est d’ordre 2? Et d’ordre 3?
Soit \(\omega\) l’ordre de la méthode.
- C’est un schéma explicite à \(3\) étages (\(s=3\)) donc \(\omega\le s=3\).
- Il est consistant (i.e. d’ordre 1) pour tout \(\sigma\) et tout \(\gamma\).
Les conditions pour être d’ordre 2 et 3 sont rappelées ci-dessous. On trouve que :
- il est d’ordre 2 ssi \(\sigma\gamma=\frac{1}{2}\),
- il est d’ordre 3 ssi, de plus, \(\sigma^2\gamma=\frac{1}{3}\) et \(\sigma\gamma\eta=\frac{1}{6}\); autrement dit, ssi \(\sigma\gamma=\frac{1}{2}\), \(\sigma=\frac{2}{3}\) et \(\eta=\frac{1}{3}\).
Conclusion : il est au moins d’ordre 2 ssi \(\sigma\gamma=\frac{1}{2}\), et d’ordre 3 ssi \(\sigma=\frac{2}{3}\), \(\gamma=\frac{3}{4}\) et \(\eta=\frac{1}{3}\).
Ci-dessous la vérification avec sympy.
Code
sigma = sp.Symbol(r"\sigma")
gamma = sp.Symbol(r"\gamma")
# =======================================================
eta = sp.Rational(1,3)
# eta = sp.Symbol(r"\eta")
# =======================================================
def ordre_RK(s, A=None, b=None, c=None):
j = sp.symbols('j')
if A is None: A = sp.MatrixSymbol('a', s, s)
else: A = sp.Matrix(A)
if c is None: c = sp.symbols(f'c_0:{s}')
if b is None: b = sp.symbols(f'b_0:{s}')
display(Markdown("**Matrice de Butcher**"))
matrice_Butcher(s, A, b, c)
display(Markdown(f"**On a {s+1} conditions pour avoir consistance = pour être d'ordre 1**"))
ordre_1(s, A, b, c, j)
display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
ordre_2(s, A, b, c, j)
display(Markdown("**On doit ajouter 2 conditions pour être d'ordre 3**"))
ordre_3(s, A, b, c, j)
# display(Markdown("**On doit ajouter 4 conditions pour être d'ordre 4**"))
# ordre_4(s, A, b, c, j)
return None
def matrice_Butcher(s, A, b, c):
But = sp.Matrix(A)
But = But.col_insert(0, sp.Matrix(c))
last = [sp.Symbol(" ")]
last.extend(b)
But = But.row_insert(s,sp.Matrix(last).T)
display(Math(sp.latex(sp.Matrix(But))))
return None
def ordre_1(s, A, b, c, j):
texte = r"\sum_{j=1}^s b_j =" + f"{sum(b).simplify()}"
texte += r"\text{ doit être égale à }1"
display(Math(texte))
for i in range(s):
somma = sp.summation(A[i,j],(j,0,s-1)).simplify()
texte = r'\sum_{j=1}^s a_{'+str(i)+r'j}=' + sp.latex( somma )
texte += r"\text{ doit être égale à }"+sp.latex(c[i])
display(Math( texte ))
return None
def ordre_2(s, A, b, c, j):
texte = r'\sum_{j=1}^s b_j c_j=' +sp.latex(sum([b[i]*c[i] for i in range(s)]).simplify())
texte += r"\text{ doit être égale à }\frac{1}{2}"
display(Math(texte))
return None
def ordre_3(s, A, b, c, j):
texte = r'\sum_{j=1}^s b_j c_j^2='
texte += sp.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='
somma = sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
texte = texte + sp.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{6}"
display(Math(texte))
return None
def ordre_4(s, A, b, c, j):
texte = r'\sum_{j=1}^s b_j c_j^3='
texte += sp.latex( sum([b[i]*c[i]**3 for i in range(s)]).simplify() )
texte += r"\text{ doit être égale à }\frac{1}{4}"
display(Math(texte))
texte = r'\sum_{i,j=1}^s b_i c_i a_{ij} c_j='
somma = sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
texte = texte + sp.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{8}"
display(Math(texte))
texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j^2='
somma = sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify()
texte = texte + sp.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{12}"
display(Math(texte))
texte = r'\sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k='
somma = sum([b[i]*A[i,j]*A[j,k]*c[k] for k in range(s) for j in range(s) for i in range(s)]).simplify()
texte = texte + sp.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{24}"
display(Math(texte))
return None
# APPLICATION A NOTRE CAS
c = [0,eta,sigma]
b = [1-gamma,0,gamma]
A = [[0,0,0],[eta,0,0],[0,sigma,0]]
s = len(c)
display(Markdown(f"La méthode de Runge-Kutta est à {s} étages"))
ordre_RK(s,A,b,c)
La méthode de Runge-Kutta est à 3 étages
Matrice de Butcher
\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\\frac{1}{3} & \frac{1}{3} & 0 & 0\\\sigma & 0 & \sigma & 0\\ & 1 - \gamma & 0 & \gamma\end{matrix}\right]\)
On a 4 conditions pour avoir consistance = pour être d’ordre 1
\(\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1\)
\(\displaystyle \sum_{j=1}^s a_{0j}=0\text{ doit être égale à }0\)
\(\displaystyle \sum_{j=1}^s a_{1j}=\frac{1}{3}\text{ doit être égale à }\frac{1}{3}\)
\(\displaystyle \sum_{j=1}^s a_{2j}=\sigma\text{ doit être égale à }\sigma\)
On doit ajouter 1 condition pour être d’ordre 2
\(\displaystyle \sum_{j=1}^s b_j c_j=\gamma \sigma\text{ doit être égale à }\frac{1}{2}\)
On doit ajouter 2 conditions pour être d’ordre 3
\(\displaystyle \sum_{j=1}^s b_j c_j^2=\gamma \sigma^{2}\text{ doit être égale à }\frac{1}{3}\)
\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{\gamma \sigma}{3}\text{ doit être égale à }\frac{1}{6}\)
Q2 [3 points] Étudier empiriquement l’ordre de convergence sur le problème de Cauchy \( \begin{cases} y'(t)=\left(\cos(2\pi t)-4t^2\right)y(t), & t\in[0;5]\\ y(0)=4 \end{cases}\) On affichera les courbes de convergence pour un schéma d’ordre 1 (à vous de donner un choix de \(\sigma\) et \(\gamma\) pour obtenir un tel schéma) et pour un schéma d’ordre 2 (ou 3 si possible).
💭 ne pas confondre
sp.cosousp.pi… (du modulesympypour le calcul de la solution exacte) etnp.cosounp.pi… (du modulenumpypour le calcul numérique).
Code
# variables globales
t0 = 0
tfinal = 5
y0 = 4
##############################################
# solution exacte
##############################################
phi_sym = lambda t,y: (4*sp.cos(2*sp.pi*t)-4*t**2)*y
t = sp.Symbol('t')
y = sp.Function('y')
edo = sp.Eq( sp.diff(y(t),t) , phi_sym(t,y(t)) )
sol = sp.dsolve(edo,y(t), ics={y(t0):y0})
display(sol)
##############################################
# solution approchée
##############################################
sol_exacte = sp.lambdify(t,sol.rhs,'numpy')
phi = lambda t,y: (4*np.cos(2*np.pi*t)-4*t**2)*y
def RK(tt,phi,y0,sigma,gamma,eta):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for i in range(len(tt)-1):
K1 = phi( tt[i] , uu[i] )
K2 = phi( tt[i]+h*eta , uu[i]+h*eta*K1 )
K3 = phi( tt[i]+h*sigma , uu[i]+h*sigma*K2 )
uu[i+1] = uu[i]+h*((1-gamma)*K1+gamma*K3)
return uu
##############################################
# ordre
##############################################
H = []
err_1, err_2, err_3 = [], [], []
plt.figure(figsize=(16,10))
ax11 = plt.subplot(3,2,1)
ax12 = plt.subplot(3,2,2)
ax21 = plt.subplot(3,2,3)
ax22 = plt.subplot(3,2,4)
ax31 = plt.subplot(3,2,5)
ax32 = plt.subplot(3,2,6)
ax11.set_title("Schéma d'ordre 1")
ax21.set_title("Schéma d'ordre 2")
ax31.set_title("Schéma d'ordre 3")
N = 0
for k in range(6):
N += 100
tt = np.linspace(t0, tfinal, N + 1)
H.append( tt[1] - tt[0] )
yy = sol_exacte(tt)
sigma, gamma, eta = 1/2, 1/2, 1/3 # cas ordre 1
uu_1 = RK(tt,phi,y0,sigma,gamma,eta)
err_1.append( np.linalg.norm(uu_1-yy,np.inf) )
ax11.plot(tt,uu_1,label=f'Approchée avec N={N}')
sigma, gamma, eta = 1/4, 2, 1/6 # cas ordre 2
uu_2 = RK(tt,phi,y0,sigma,gamma,eta)
err_2.append( np.linalg.norm(uu_2-yy,np.inf) )
ax21.plot(tt,uu_2,label=f'Approchée avec N={N}')
sigma, gamma, eta = 2/3, 3/4, 1/3 # cas ordre 3
uu_3 = RK(tt,phi,y0,sigma,gamma,eta)
err_3.append( np.linalg.norm(uu_3-yy,np.inf) )
ax31.plot(tt,uu_3,label=f'Approchée avec N={N}')
ax11.plot(tt,yy,label='Exacte')
plt.xlabel('$t$')
plt.ylabel('$y$')
ax11.grid(True)
ax11.legend()
ax12.plot( np.log(H), np.log(err_1), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
ax12.set_title(f'RK : ordre = {np.polyfit(np.log(H),np.log(err_1),1)[0]:1.2f}')
ax12.grid(True)
ax21.plot(tt,yy,label='Exacte')
plt.xlabel('$t$')
plt.ylabel('$y$')
ax21.grid(True)
ax21.legend()
ax22.plot( np.log(H), np.log(err_2), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
ax22.set_title(f'RK : ordre = {np.polyfit(np.log(H),np.log(err_2),1)[0]:1.2f}')
ax22.grid(True)
ax31.plot(tt,yy,label='Exacte')
plt.xlabel('$t$')
plt.ylabel('$y$')
ax31.grid(True)
ax31.legend()
ax32.plot( np.log(H), np.log(err_3), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
ax32.set_title(f'RK : ordre = {np.polyfit(np.log(H),np.log(err_3),1)[0]:1.2f}')
ax32.grid(True)\(\displaystyle y{\left(t \right)} = 4 e^{- \frac{4 t^{3}}{3} + \frac{2 \sin{\left(2 \pi t \right)}}{\pi}}\)
Q3 [2 points] Étudier théoriquement la A-stabilité du schéma d’ordre maximal (dans la correction ci-dessous \(\eta=\frac{1}{3}\) et donc l’ordre maximal est 3).
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\\ K_2 = -\beta \left(u_n+\frac{h}{3}K_1\right)\leadsto K_2=\left(1-\frac{\beta h}{3}\right)K_1\\ K_3 = -\beta \left(u_n+\frac{2h}{3}K_2\right)\leadsto K_3=\left(1-\frac{2\beta h}{3}\left(1-\frac{\beta h}{3}\right)\right)K_1\\ u_{n+1} = u_n + \frac{h}{4}\left(K_1+4K_3\right) & n=0,1,\dots N-1 \end{cases} \) L’étude ci-dessous montre que le schéma est A-stable ssi \(h<\frac{\vartheta}{\beta}\) avec \(\vartheta\approx2.7\).
Code
u_n = sp.Symbol('u_n')
h = sp.Symbol('h')
K1 = sp.Symbol('K1')
K2 = sp.Symbol('K2')
beta = sp.Symbol(r'\beta')
sigma = sp.Symbol(r"\sigma")
gamma = sp.Symbol(r"\gamma")
eta = sp.Symbol(r"\eta")
g = lambda y : -beta*y
K1 = g(u_n)
K2 = g(u_n+h*eta*K1)
K3 = g(u_n+h*sigma*K2)
RHS = u_n+h*((1-gamma)*K1+gamma*K3)
RHS = RHS.factor()
texte = 'u_{n+1}=' + sp.latex(RHS)
display(Math(texte))
display(Markdown(r"Notons $x=h\beta>0$. On trouve donc"))
x = sp.Symbol('x',positive=True)
RHS = RHS.subs(h*beta,x).simplify()
texte = 'u_{n+1}=' + sp.latex(RHS)
display(Math(texte))
texte = r"La méthode est A-stable ssi $|q(x)|<1$ avec $q(x)=\dfrac{u_{n+1}}{u_n}="
texte += sp.latex(RHS/u_n) + "$"
display(Markdown(texte))
# ======== On remplace les paramètres par les valeurs optimales ==========
RHS = RHS.subs(gamma,sp.Rational(4,3)).subs(sigma,sp.Rational(2,3)).subs(eta,sp.Rational(1,3)).simplify()
q = sp.apart(RHS/u_n)
# ======================================================================
texte = r"Dans le cas optimal on a $\sigma=\dfrac{2}{3}$ et $\gamma=\dfrac{4}{3}$ ainsi \(q(x)=\dfrac{u_{n+1}}{u_n}="
texte += sp.latex(RHS/u_n) + "\)"
display(Markdown(texte))
display(Markdown(r"Étudions la fonction $q(x)$ pour $x>0$:"))
sp.solve( [ -1<RHS/u_n , RHS/u_n<1 ] )
display(Markdown(f"$q(0)={q.subs(x,0)}$"))
ell = sp.Limit(q,x,sp.oo)
display(Math(f"{sp.latex(ell)}={sp.latex(ell.doit())}"))
dq = sp.diff(q,x).simplify()
display(Markdown(f"$\\displaystyle q'(x)={sp.latex(dq)}$"))
sol=sp.solveset(dq,domain=sp.S.Reals)
display(Markdown(f"$q'(x)<0$ pour tout $x>0$ et l'on a"))
sol = sp.solveset(q+1,domain=sp.S.Reals)
texte = r"q(x)=-1 \iff x= " + sp.latex(sol) + r"\approx" + sp.latex(sol.evalf())
display( Math( texte ) )
# p1 = sp.plot(q,(x,0,5),ylim=(-2,2),xlabel=r'$x$',ylabel=r'$q(x)$',show=False,grid=True)
# p2 = sp.plot(1,(x,0,5),ylim=(-2,2),show=False,line_color='r')
# p3 = sp.plot(-1,(x,0,5),ylim=(-2,2),show=False,line_color='r')
# p1.extend(p2)
# p1.extend(p3)
# p1.show()
plt.figure(figsize=(12,6))
q = sp.lambdify(x,RHS/u_n,'numpy')
xx = np.linspace(0,5,100)
plt.plot(xx,q(xx),label=r'$q(x)$')
plt.hlines(1,0,5,'r')
plt.hlines(-1,0,5,'r')
plt.axis([0,5,-2,2])
plt.grid();\(\displaystyle u_{n+1}=- u_{n} \left(\beta^{3} \eta \gamma \sigma h^{3} - \beta^{2} \gamma \sigma h^{2} + \beta h - 1\right)\)
Notons \(x=h\beta>0\). On trouve donc
\(\displaystyle u_{n+1}=u_{n} \left(- \eta \gamma \sigma x^{3} + \gamma \sigma x^{2} - x + 1\right)\)
La méthode est A-stable ssi \(|q(x)|<1\) avec \(q(x)=\dfrac{u_{n+1}}{u_n}=- \eta \gamma \sigma x^{3} + \gamma \sigma x^{2} - x + 1\)
Dans le cas optimal on a \(\sigma=\dfrac{2}{3}\) et \(\gamma=\dfrac{4}{3}\) ainsi \(q(x)=\dfrac{u_{n+1}}{u_n}=- \frac{8 x^{3}}{27} + \frac{8 x^{2}}{9} - x + 1\)
Étudions la fonction \(q(x)\) pour \(x>0\):
\(q(0)=1\)
\(\displaystyle \lim_{x \to \infty}\left(- \frac{8 x^{3}}{27} + \frac{8 x^{2}}{9} - x + 1\right)=-\infty\)
\(\displaystyle q'(x)=- \frac{8 x^{2}}{9} + \frac{16 x}{9} - 1\)
\(q'(x)<0\) pour tout \(x>0\) et l’on a
\(\displaystyle q(x)=-1 \iff x= \left\{- \frac{1}{24 \sqrt[3]{\frac{43}{432} + \frac{\sqrt{822}}{288}}} + 1 + 3 \sqrt[3]{\frac{43}{432} + \frac{\sqrt{822}}{288}}\right\}\approx\left\{2.68038081801272\right\}\)
