Code
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from scipy.optimize import fsolve
from IPython.display import display, Markdown, Math2024 Contrôle Terminal session 1
⏱ 2h30, deposer le notebook dans Moodle
EXERCICE 1 : Choix d’un schéma
Soient \(\beta, \varepsilon\in\mathbb{R}\) et considérons le problème de Cauchy suivant
\( \begin{cases} y'(t) = \beta \Big( \sin(t)-y(t)\Big) + \cos(t) &\text{pour } t>0, \\ y(0) = \varepsilon. \end{cases} \)
- Calculer la solution exacte en fonction de \(\beta\) et \(\varepsilon\).
- Si \(\beta = -5\) et \(\varepsilon=0\), quel type de schéma doit-on privilégier pour résoudre ce problème de Cauchy ?
Méthode 1 On peut résoudre ce problème de Cauchy directement en remarquant que l’edo se réécrit comme
\(
y'(t) - \cos(t) = -\beta \big(y(t)-\sin(t)\big),
\)
soit encore
\(
\big(y(t)-\sin(t)\big)' = -\beta \big(y(t)-\sin(t)\big).
\)
Si on pose \(z(t)=y(t)-\sin(t)\) on a le problème de Cauchy
\(
z'(t) = -\beta z(t)
\)
avec \(z(0) = y(0)-0=\varepsilon\). La solution est donc \(z(t) = z(0) e^{-\beta t}= \varepsilon e^{-\beta t}\) et donc
\(
\boxed{y(t)} = z(t)+\sin(t) \boxed{= \sin(t)+\varepsilon e^{-\beta t}}
\)
Méthode 2 On peut résoudre le problème de Cauchy comme en L1 : la solution générale de l’EDO \(a(t)y'(t)+b(t)y(t)=g(t)\) est donnée par
\(
y(t) = C e^{-A(t)} + B(t) e^{-A(t)}
\)
où \(A(t) = \int^t \frac{b(s)}{a(s)}ds\) et \(B(t) = \int^t \frac{g(s)}{a(s)}e^{A(s)}ds\). Dans notre cas, on a \(y'(t) +\beta y(t) = \beta \sin(t) + \cos(t)\), donc \(A(t) = \int^t \beta ds = \beta t\) et \(B(t) = \int^t e^{\beta s}(\beta \sin(s) + \cos(s))ds=e^{\beta t}\sin(t)\). La solution générale est alors \(y(t) = C e^{-\beta t} + \sin(t)\). En utilisant la condition initiale \(y(0)=\varepsilon\), on trouve \(C=\varepsilon\) et enfin \(y(t) = \sin(t)+\varepsilon e^{-\beta t}\).
Méthode 3 On peut calculer la solution en utilisant sympy :
Code
t = sp.symbols('t')
y = sp.Function('y')(t)
beta = sp.symbols(r'\beta', real=True)
epsilon = sp.symbols(r'\varepsilon', real=True)
dydt = y.diff(t)
phi_sp = lambda t,z : beta*(sp.sin(t)-z)+sp.cos(t)
ode = sp.Eq(dydt, phi_sp(t,y) )
sol = sp.dsolve(ode, y, ics={y.subs(t, 0) : epsilon})
sol\(\displaystyle y{\left(t \right)} = \varepsilon e^{- \beta t} + \sin{\left(t \right)}\)
Si \(\beta<0\), il s’agit d’un problème de Cauchy numériquement mal posé (c’est-à-dire, il est très sensible aux perturbations sur la donnée initiale). En effet,
- si \(\varepsilon>0\), \(y(t) \sim \varepsilon e^{-\beta t} \xrightarrow[t\to\infty]{} +\infty\),
- si \(\varepsilon<0\), \(y(t) \sim \varepsilon e^{-\beta t} \xrightarrow[t\to\infty]{} -\infty\),
- si \(\varepsilon=0\), \(y(t) = \sin(t) \in [-1;1]\) pour tout \(t\).
Une toute (petite) erreur de calcul aura le même effet qu’une perturbation de la condition initiale : on “réveille” le terme exponentielle. On doit donc choisir un schéma d’ordre élevé car les erreurs d’approximation amènent toujours à résoudre un problème perturbé.
Code
t0 = 0
tfinal = 10
phi = lambda t,y: -5*(np.sin(t)-y) + np.cos(t)
N = 200
tt = np.linspace(t0, tfinal, N + 1)
func = sp.lambdify(t, sol.rhs.subs({epsilon:0.5, beta:5}), 'numpy')
yy = func(tt)
plt.figure(figsize=(10,5))
g1 = np.linspace(0,10,21)
g2 = np.linspace(-3,3,21)
T,Y = np.meshgrid(g1,g2)
DT, DY = 1, phi(T,Y)
plt.plot(tt,yy,color='red')
M = np.hypot(DT,DY)
plt.streamplot(T,Y, DT/M, DY/M, color=M, linewidth=0.5, cmap=plt.cm.viridis, density=2, arrowstyle='->', arrowsize=1.5)
plt.grid()
plt.xlabel('t')
plt.ylabel('y')
plt.title('Lignes de courant');
Pour bien comprendre, allons voir l’influence de l’ordre du schéma en comparant les approximations obtenues avec le schéma d’Euler explicite et un schéma RK aussi explicite d’ordre 4 dans le cas \(\varepsilon=0\) et \(\beta=-5\). On remarque qu’initialement les approximations sont proches de la solution exacte, mais l’erreur d’approximation du schéma d’Euler explicite fait que la solution approchée correspond plutôt à une solution ayant \(\varepsilon>0\). Le schéma RK est plus précis mais finit par diverger lui aussi.
Code
tt = np.linspace(0, 5, 300)
fig, (ax1, ax2) = plt.subplots(nrows=1,ncols=2,figsize=(16, 6))
# Solutions exactes pour différentes valeurs de epsilon et beta = -5
for epsi in [-1e-3, -1e-8, 0, 1e-8, 1e-3]:
y_sol = sp.lambdify(t, sol.rhs.subs({beta:-5,epsilon:epsi}), 'numpy')
ax1.plot(tt, y_sol(tt), label=fr'$\varepsilon={epsi}$')
ax1.set_xlabel('t')
ax1.set_ylabel('y')
ax1.grid()
ax1.legend()
ax1.axes.set_ylim(-5, 5)
ax1.set_title(r'Solutions exactes pour différents $\varepsilon$')
# Approximations du cas epsilon = 0 avec un schéma d'ordre 1 et un schéma d'ordre 4
def EE(phi, tt, y0):
h = tt[1] - tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for i in range(len(tt) - 1):
uu[i+1] = uu[i] + h*phi(tt[i],uu[i])
return uu
def RK4(phi,tt,y0):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
k1 = phi( tt[n], uu[n] )
k2 = phi( tt[n]+h/2 , uu[n]+k1*h/2 )
k3 = phi( tt[n]+h/2 , uu[n]+h*k2/2 )
k4 = phi( tt[n+1] , uu[n]+h*k3 )
uu[n+1] = uu[n] + (k1+2*k2+2*k3+k4)*h/6
return uu
phi = lambda t,y : -5*(np.sin(t)-y)+np.cos(t)
y0 = 0
ax2.plot(tt, EE(phi, tt, y0), ':', label='Euler explicite')
ax2.plot(tt, RK4(phi, tt, y0), ':', label='RK4')
y_sol = sp.lambdify(t, sol.rhs.subs({beta:-5,epsilon:0}), 'numpy')
ax2.plot(tt, y_sol(tt), label='Exacte')
ax2.set_xlabel('t')
ax2.set_ylabel('y')
ax2.grid()
ax2.legend()
ax2.axes.set_ylim(-5, 5)
ax2.set_title(r"Solution exacte VS approximations avec EE et RK4 pour $\varepsilon=0$");
EXERCICE 2 : Schémas MPL à 3 pas d’ordre maximal
Écrire le schéma à \(q=3\) pas linéaire (il comportera donc 7 coefficients).
Il n’est pas nécessaire de répondre à cette question pour la suite de l’exercice.
Les 7 équations linéaires \(\xi(i)=1\) pour \(i=0,\dots,6\) constituent un système linéaire dont les 7 inconnues sont les coefficients du schéma. Calculer les coefficients de ce système consistant d’ordre 6 et montrer qu’il n’est pas zéro-stable.Sans effectuer de calculs, répondre aux deux questions suivantes :
- quel est l’ordre maximal d’un schéma à 3 pas implicite convergent ?
- quel est l’ordre maximal d’un schéma à 3 pas explicite convergent ?
Dans la suite on fixera \(b_{-1}=0\) et \(b_1 = -4/3\). De plus, \(a_1\) sera fixé en fonction de votre prénom et nom (pour cela, remplacez mon prénom et mon nom par les vôtres dans la cellule ci-dessous).
Déterminer les 4 paramètres manquants pour que le schéma soit d’ordre maximal et convergent.Tester empiriquement l’ordre de convergence du schéma ainsi obtenu sur le problème de Cauchy \(y'(t) = \sin^2(t)y(t)\) et \(y(0) = 1\) pour \(t\in[0,2]\).
Code
# ================================================
# On fixe $a_1$ pour le point 4
# ================================================
nom = "Faccanoni"
prenom = "Gloria"
S = sum([ord(c) for c in nom+prenom])
Lista = [ sp.Rational(num,12) for num in range(-24,13) ]
idx = S%len(Lista)
display(Markdown(f"**{nom} {prenom}** : $a_1 = {Lista[idx]}$"))Faccanoni Gloria : \(a_1 = 0\)
Un schéma implicite à \(q=3\) étapes s’écrit \( u_{n+1}=a_0u_n+a_1u_{n-1}+a_2u_{n-2}+h\big(b_{-1}\varphi_{n+1}+b_0\varphi_n+b_1\varphi_{n-1}+b_2\varphi_{n-2}\big) \)
Il contient les 7 paramètres \(a_0\), \(a_1\), \(a_2\), \(b_{-1}\), \(b_0\), \(b_1\), \(b_2\), il est donc théoriquement possible de construire un schéma d’ordre 6 en résolvant le système linéaire de 7 équations et 7 inconnues suivant
\(
\begin{cases}
\xi(0) = a_{0} + a_{1} + a_{2}=1 \\
\xi(1) = - a_{1} - 2 a_{2} + b_{0} + b_{1} + b_{2} + b_{-1}=1 \\
\xi(2) = a_{1} + 4 a_{2} - 2 \left(b_{1} + 2 b_{2} - b_{-1}\right)=1 \\
\xi(3) = - a_{1} - 8 a_{2} + 3 \left(b_{1} + 4 b_{2} + b_{-1}\right)=1 \\
\xi(4) = a_{1} + 16 a_{2} - 4 \left(b_{1} + 8 b_{2} - b_{-1}\right)=1 \\
\xi(5) = - a_{1} - 32 a_{2} + 5 \left(b_{1} + 16 b_{2} + b_{-1}\right)=1 \\
\xi(6) = a_{1} + 64 a_{2} - 6 \left(b_{1} + 32 b_{2} - b_{-1}\right)=1
\end{cases}
\)
Code
q = 3
def xi(i, q):
sa = sum((-j)**i * aa[j] for j in range(q))
sb = b_m1 + sum((-j)**(i-1) * bb[j] for j in range(1, q))
if i == 1:
sb += bb[0]
return (sa).factor() + (i * sb).factor()
# Ordre maximal selon Dalquist
omega = q + 2 if q%2==0 else q + 1 # i = 0, ..., omega
# Définition des symboles
a_0, a_1, a_2 = sp.symbols('a_0 a_1 a_2', real=True)
b_0, b_1, b_2 = sp.symbols('b_0 b_1 b_2', real=True)
b_m1 = sp.symbols('b_{-1}', real=True)
aa = [a_0, a_1, a_2]
bb = [b_0, b_1, b_2]
# display(Markdown(f"Pour une méthode à q={q} pas, on affiche $\\xi(i)$ pour i = 0, ..., {2*q+1}"))
# for i in range(7):
# display(Markdown(f"$\\xi({i}) = {sp.latex(xi(i, q))}$"))
display(Markdown(f"Une méthode à q={q} pas contient {2*q+1} paramètres, on peut donc résoudre le système de {2*q+1} équations suivant (garantissant la consistance d'ordre {2*q}) pour les déterminer :"))
systeme = [sp.Eq(xi(i, q),1) for i in range(7)]
display(Markdown(r"$\begin{cases}"+' '.join([sp.latex(s)+r"\\" for s in systeme]) + r"\end{cases}$"))
Une méthode à q=3 pas contient 7 paramètres, on peut donc résoudre le système de 7 équations suivant (garantissant la consistance d’ordre 6) pour les déterminer :
\(\begin{cases}a_{0} + a_{1} + a_{2} = 1\\ - a_{1} - 2 a_{2} + b_{0} + b_{1} + b_{2} + b_{-1} = 1\\ a_{1} + 4 a_{2} - 2 \left(b_{1} + 2 b_{2} - b_{-1}\right) = 1\\ - a_{1} - 8 a_{2} + 3 \left(b_{1} + 4 b_{2} + b_{-1}\right) = 1\\ a_{1} + 16 a_{2} - 4 \left(b_{1} + 8 b_{2} - b_{-1}\right) = 1\\ - a_{1} - 32 a_{2} + 5 \left(b_{1} + 16 b_{2} + b_{-1}\right) = 1\\ a_{1} + 64 a_{2} - 6 \left(b_{1} + 32 b_{2} - b_{-1}\right) = 1\\\end{cases}\)
Calculons la solution de ce système linéaire:
Code
sol = sp.solve(systeme)
display(sol){a_0: -27/11,
a_1: 27/11,
a_2: 1,
b_0: 27/11,
b_1: 27/11,
b_2: 3/11,
b_{-1}: 3/11}Cépendant on sait, d’après la première barrière de Dalquist, qu’un schéma multipas linéaire implicite à \(q\) pas avec \(q\) impair est au mieux d’ordre \(q+1\). Donc un schéma à 3 pas ne peut pas être d’ordre 6. On peut le vérifier en calculant le polynôme caractéristique du schéma : \( \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^3-a_0r^2-a_1r-a_2=(r-1)(r^2+\frac{38}{11}r+1). \)
La méthode n’est pas zéro-stable car une des racines n’est pas dans le cercle unité (voir le calcul des racines ci-dessous) :
Code
r = sp.Symbol('r')
p = q-1
rho = sp.poly(r**(p+1) - sum( aa[j].subs(sol) * r**(p-j) for j in range(p+1)), r)
display(Markdown(f"Le premier polynome caractéristique est :"))
txt = sp.latex(sp.Eq(sp.symbols(r'\varrho(r)'),rho.as_expr()),order='old')
txt += "=" + sp.latex(rho.as_expr().factor(),order='old')
display(Math(txt))
display(Markdown(f"Ses racines sont :"))
racines = sp.solve(rho,r)
display(racines)
display(Markdown(f"dont les modules valent :"))
display([abs(rac).evalf() for rac in racines])
Le premier polynome caractéristique est :
\(\displaystyle \varrho(r) = -1 - \frac{27 r}{11} + \frac{27 r^{2}}{11} + r^{3}=\frac{\left(-1 + r\right) \left(11 + 38 r + 11 r^{2}\right)}{11}\)
Ses racines sont :
[1, -19/11 - 4*sqrt(15)/11, -19/11 + 4*sqrt(15)/11]dont les modules valent :
[1.00000000000000, 3.13563030771179, 0.318915146833667]La première barrière de Dahlquist affirme qu’un schéma à \(q=3\) pas
- s’il est implicite, il est au mieux d’ordre \(q+1=4\);
- s’il est explicite, il est au mieux d’ordre \(q=3\).
Code
for ell in Lista:
choix = {b_m1 : 0 , a_1 : ell , b_1 : sp.Rational(-4,3)}
systeme = [sp.Eq(xi(i, q).subs(choix),1) for i in range(4)]
sol = sp.solve(systeme, [a_0, a_2, b_0, b_2])
display(Markdown(f"Choix : ${sp.latex(choix)}$"))
display(Markdown(r"Système : $\begin{cases}"+' '.join([sp.latex(s)+r"\\" for s in systeme]) + r"\end{cases}$"))
display(Markdown(f"Solution : ${sp.latex(sol)}$"))
display(Markdown("=============================="))
Choix : \(\left\{ a_{1} : -2, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - 2 = 1\\ - 2 a_{2} + b_{0} + b_{2} + \frac{2}{3} = 1\\ 4 a_{2} - 4 b_{2} + \frac{2}{3} = 1\\ - 8 a_{2} + 12 b_{2} - 2 = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : 2, \ a_{2} : 1, \ b_{0} : \frac{17}{12}, \ b_{2} : \frac{11}{12}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{23}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{23}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} + \frac{7}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{3}{4} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{25}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{47}{24}, \ a_{2} : \frac{23}{24}, \ b_{0} : \frac{23}{16}, \ b_{2} : \frac{43}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{11}{6}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{11}{6} = 1\\ - 2 a_{2} + b_{0} + b_{2} + \frac{1}{2} = 1\\ 4 a_{2} - 4 b_{2} + \frac{5}{6} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{13}{6} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{23}{12}, \ a_{2} : \frac{11}{12}, \ b_{0} : \frac{35}{24}, \ b_{2} : \frac{7}{8}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{7}{4}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{7}{4} = 1\\ - 2 a_{2} + b_{0} + b_{2} + \frac{5}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{11}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{9}{4} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{15}{8}, \ a_{2} : \frac{7}{8}, \ b_{0} : \frac{71}{48}, \ b_{2} : \frac{41}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{5}{3}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{5}{3} = 1\\ - 2 a_{2} + b_{0} + b_{2} + \frac{1}{3} = 1\\ 4 a_{2} - 4 b_{2} + 1 = 1\\ - 8 a_{2} + 12 b_{2} - \frac{7}{3} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{11}{6}, \ a_{2} : \frac{5}{6}, \ b_{0} : \frac{3}{2}, \ b_{2} : \frac{5}{6}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{19}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{19}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} + \frac{1}{4} = 1\\ 4 a_{2} - 4 b_{2} + \frac{13}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{29}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{43}{24}, \ a_{2} : \frac{19}{24}, \ b_{0} : \frac{73}{48}, \ b_{2} : \frac{13}{16}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{3}{2}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{3}{2} = 1\\ - 2 a_{2} + b_{0} + b_{2} + \frac{1}{6} = 1\\ 4 a_{2} - 4 b_{2} + \frac{7}{6} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{5}{2} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{7}{4}, \ a_{2} : \frac{3}{4}, \ b_{0} : \frac{37}{24}, \ b_{2} : \frac{19}{24}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{17}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{17}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} + \frac{1}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{5}{4} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{31}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{41}{24}, \ a_{2} : \frac{17}{24}, \ b_{0} : \frac{25}{16}, \ b_{2} : \frac{37}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{4}{3}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{4}{3} = 1\\ - 2 a_{2} + b_{0} + b_{2} = 1\\ 4 a_{2} - 4 b_{2} + \frac{4}{3} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{8}{3} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{5}{3}, \ a_{2} : \frac{2}{3}, \ b_{0} : \frac{19}{12}, \ b_{2} : \frac{3}{4}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{5}{4}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{5}{4} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{1}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{17}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{11}{4} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{13}{8}, \ a_{2} : \frac{5}{8}, \ b_{0} : \frac{77}{48}, \ b_{2} : \frac{35}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{7}{6}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{7}{6} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{1}{6} = 1\\ 4 a_{2} - 4 b_{2} + \frac{3}{2} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{17}{6} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{19}{12}, \ a_{2} : \frac{7}{12}, \ b_{0} : \frac{13}{8}, \ b_{2} : \frac{17}{24}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{13}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{13}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{1}{4} = 1\\ 4 a_{2} - 4 b_{2} + \frac{19}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{35}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{37}{24}, \ a_{2} : \frac{13}{24}, \ b_{0} : \frac{79}{48}, \ b_{2} : \frac{11}{16}\right\}\)
==============================
Choix : \(\left\{ a_{1} : -1, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - 1 = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{1}{3} = 1\\ 4 a_{2} - 4 b_{2} + \frac{5}{3} = 1\\ - 8 a_{2} + 12 b_{2} - 3 = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{3}{2}, \ a_{2} : \frac{1}{2}, \ b_{0} : \frac{5}{3}, \ b_{2} : \frac{2}{3}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{11}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{11}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{5}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{7}{4} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{37}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{35}{24}, \ a_{2} : \frac{11}{24}, \ b_{0} : \frac{27}{16}, \ b_{2} : \frac{31}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{5}{6}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{5}{6} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{1}{2} = 1\\ 4 a_{2} - 4 b_{2} + \frac{11}{6} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{19}{6} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{17}{12}, \ a_{2} : \frac{5}{12}, \ b_{0} : \frac{41}{24}, \ b_{2} : \frac{5}{8}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{3}{4}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{3}{4} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{7}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{23}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{13}{4} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{11}{8}, \ a_{2} : \frac{3}{8}, \ b_{0} : \frac{83}{48}, \ b_{2} : \frac{29}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{2}{3}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{2}{3} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{2}{3} = 1\\ 4 a_{2} - 4 b_{2} + 2 = 1\\ - 8 a_{2} + 12 b_{2} - \frac{10}{3} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{4}{3}, \ a_{2} : \frac{1}{3}, \ b_{0} : \frac{7}{4}, \ b_{2} : \frac{7}{12}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{7}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{7}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{3}{4} = 1\\ 4 a_{2} - 4 b_{2} + \frac{25}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{41}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{31}{24}, \ a_{2} : \frac{7}{24}, \ b_{0} : \frac{85}{48}, \ b_{2} : \frac{9}{16}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{1}{2}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{1}{2} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{5}{6} = 1\\ 4 a_{2} - 4 b_{2} + \frac{13}{6} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{7}{2} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{5}{4}, \ a_{2} : \frac{1}{4}, \ b_{0} : \frac{43}{24}, \ b_{2} : \frac{13}{24}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{5}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{5}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{11}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{9}{4} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{43}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{29}{24}, \ a_{2} : \frac{5}{24}, \ b_{0} : \frac{29}{16}, \ b_{2} : \frac{25}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{1}{3}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{1}{3} = 1\\ - 2 a_{2} + b_{0} + b_{2} - 1 = 1\\ 4 a_{2} - 4 b_{2} + \frac{7}{3} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{11}{3} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{7}{6}, \ a_{2} : \frac{1}{6}, \ b_{0} : \frac{11}{6}, \ b_{2} : \frac{1}{2}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{1}{4}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{1}{4} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{13}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{29}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{15}{4} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{9}{8}, \ a_{2} : \frac{1}{8}, \ b_{0} : \frac{89}{48}, \ b_{2} : \frac{23}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{1}{6}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{1}{6} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{7}{6} = 1\\ 4 a_{2} - 4 b_{2} + \frac{5}{2} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{23}{6} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{13}{12}, \ a_{2} : \frac{1}{12}, \ b_{0} : \frac{15}{8}, \ b_{2} : \frac{11}{24}\right\}\)
==============================
Choix : \(\left\{ a_{1} : - \frac{1}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} - \frac{1}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{5}{4} = 1\\ 4 a_{2} - 4 b_{2} + \frac{31}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{47}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{25}{24}, \ a_{2} : \frac{1}{24}, \ b_{0} : \frac{91}{48}, \ b_{2} : \frac{7}{16}\right\}\)
==============================
Choix : \(\left\{ a_{1} : 0, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{4}{3} = 1\\ 4 a_{2} - 4 b_{2} + \frac{8}{3} = 1\\ - 8 a_{2} + 12 b_{2} - 4 = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : 1, \ a_{2} : 0, \ b_{0} : \frac{23}{12}, \ b_{2} : \frac{5}{12}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{1}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{1}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{17}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{11}{4} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{49}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{23}{24}, \ a_{2} : - \frac{1}{24}, \ b_{0} : \frac{31}{16}, \ b_{2} : \frac{19}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{1}{6}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{1}{6} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{3}{2} = 1\\ 4 a_{2} - 4 b_{2} + \frac{17}{6} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{25}{6} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{11}{12}, \ a_{2} : - \frac{1}{12}, \ b_{0} : \frac{47}{24}, \ b_{2} : \frac{3}{8}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{1}{4}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{1}{4} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{19}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{35}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{17}{4} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{7}{8}, \ a_{2} : - \frac{1}{8}, \ b_{0} : \frac{95}{48}, \ b_{2} : \frac{17}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{1}{3}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{1}{3} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{5}{3} = 1\\ 4 a_{2} - 4 b_{2} + 3 = 1\\ - 8 a_{2} + 12 b_{2} - \frac{13}{3} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{5}{6}, \ a_{2} : - \frac{1}{6}, \ b_{0} : 2, \ b_{2} : \frac{1}{3}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{5}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{5}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{7}{4} = 1\\ 4 a_{2} - 4 b_{2} + \frac{37}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{53}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{19}{24}, \ a_{2} : - \frac{5}{24}, \ b_{0} : \frac{97}{48}, \ b_{2} : \frac{5}{16}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{1}{2}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{1}{2} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{11}{6} = 1\\ 4 a_{2} - 4 b_{2} + \frac{19}{6} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{9}{2} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{3}{4}, \ a_{2} : - \frac{1}{4}, \ b_{0} : \frac{49}{24}, \ b_{2} : \frac{7}{24}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{7}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{7}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{23}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{13}{4} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{55}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{17}{24}, \ a_{2} : - \frac{7}{24}, \ b_{0} : \frac{33}{16}, \ b_{2} : \frac{13}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{2}{3}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{2}{3} = 1\\ - 2 a_{2} + b_{0} + b_{2} - 2 = 1\\ 4 a_{2} - 4 b_{2} + \frac{10}{3} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{14}{3} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{2}{3}, \ a_{2} : - \frac{1}{3}, \ b_{0} : \frac{25}{12}, \ b_{2} : \frac{1}{4}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{3}{4}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{3}{4} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{25}{12} = 1\\ 4 a_{2} - 4 b_{2} + \frac{41}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{19}{4} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{5}{8}, \ a_{2} : - \frac{3}{8}, \ b_{0} : \frac{101}{48}, \ b_{2} : \frac{11}{48}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{5}{6}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{5}{6} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{13}{6} = 1\\ 4 a_{2} - 4 b_{2} + \frac{7}{2} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{29}{6} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{7}{12}, \ a_{2} : - \frac{5}{12}, \ b_{0} : \frac{17}{8}, \ b_{2} : \frac{5}{24}\right\}\)
==============================
Choix : \(\left\{ a_{1} : \frac{11}{12}, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + \frac{11}{12} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{9}{4} = 1\\ 4 a_{2} - 4 b_{2} + \frac{43}{12} = 1\\ - 8 a_{2} + 12 b_{2} - \frac{59}{12} = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{13}{24}, \ a_{2} : - \frac{11}{24}, \ b_{0} : \frac{103}{48}, \ b_{2} : \frac{3}{16}\right\}\)
==============================
Choix : \(\left\{ a_{1} : 1, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
Système : \(\begin{cases}a_{0} + a_{2} + 1 = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{7}{3} = 1\\ 4 a_{2} - 4 b_{2} + \frac{11}{3} = 1\\ - 8 a_{2} + 12 b_{2} - 5 = 1\\\end{cases}\)
Solution : \(\left\{ a_{0} : \frac{1}{2}, \ a_{2} : - \frac{1}{2}, \ b_{0} : \frac{13}{6}, \ b_{2} : \frac{1}{6}\right\}\)
==============================
Le cas correspondant à mon prénom et nom est le suivant :
Code
# voir la cellule initiale pour fixer le paramètre
choix = {b_m1 : 0 , a_1 : Lista[idx] , b_1 : sp.Rational(-4,3)}
display(Markdown(f"Pour le choix ${sp.latex(choix)}$"))
systeme = [sp.Eq(xi(i, q).subs(choix),1) for i in range(4)]
display(Markdown("le système de 4 équations devient :"))
display(Markdown(r"$\begin{cases}"+' '.join([sp.latex(s)+r"\\" for s in systeme]) + r"\end{cases}$"))
Pour le choix \(\left\{ a_{1} : 0, \ b_{1} : - \frac{4}{3}, \ b_{-1} : 0\right\}\)
le système de 4 équations devient :
\(\begin{cases}a_{0} + a_{2} = 1\\ - 2 a_{2} + b_{0} + b_{2} - \frac{4}{3} = 1\\ 4 a_{2} - 4 b_{2} + \frac{8}{3} = 1\\ - 8 a_{2} + 12 b_{2} - 4 = 1\\\end{cases}\)
Ces 4 équations permettent d’exprimer \(b_0\), \(b_2\), \(a_0\) et \(a_2\) en fonction de \(a_1\), \(b_1\) et \(b_{-1}\):
Code
sol = sp.solve(systeme, [a_0, a_2, b_0, b_2])
# Un seul dictionnaire avec toutes les valeurs:
sol.update(choix)
display(sol){a_0: 1, a_2: 0, b_0: 23/12, b_2: 5/12, b_{-1}: 0, a_1: 0, b_1: -4/3}Vérifions si le schéma ainsi obtenu est zéro-stable en calculant les racines du premier polynôme caractéristique.
Le premier polynôme caractéristique est \( \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^3-r^2=r^3-a_0r^2-a_1r-a_2. \)
La méthode est zéro-stable ssi ses racines sont de module inférieur ou égal à 1 et si elles sont simples lorsqu’elles ont module égal à 1. Ce qui est bien le cas ici :
Code
r = sp.Symbol('r')
p = q-1
rho = sp.poly(r**(p+1) - sum( aa[j].subs(sol).subs(choix) * r**(p-j) for j in range(p+1)), r)
display(Markdown(f"Le premier polynome caractéristique est :"))
display(Math(sp.latex(rho.as_expr() ,order='old')))
display(Markdown(f"Ses racines, avec leur multiplicité, sont :"))
racines = sp.roots(rho,r)
display(racines)
display([float(abs(rac)) for rac in racines.keys()])Le premier polynome caractéristique est :
\(\displaystyle - r^{2} + r^{3}\)
Ses racines, avec leur multiplicité, sont :
{1: 1, 0: 2}[1.0, 0.0]Remarque : juste en considérant la condition de consistance, on peut factoriser le polynôme caractéristique par \(r-1\), car \(a_0 = 1-a_1-a_2\) ainsi :
\(\begin{align*} \varrho(r) &=(r-1)\left(r^2+\left(a_2-\frac{1}{2}\right)r+a_2\right) \\ &=(r-1)\left(r-\frac{1-2a_2-\sqrt{4a_2^2-20a_2+1}}{4}\right)\left(r-\frac{1-2a_2+\sqrt{4a_2^2-20a_2+1}}{4}\right) \end{align*}\)
Code
# AUTOMATISATION : convertion des valeurs sympy en nombres flottants de NumPy et affectation
dico = {cle: float(valeur) for cle, valeur in sol.items()}
for cle, valeur in dico.items():
globals()[str(cle)] = valeur
##################################################
# EDO
##################################################
t0 = 0
tfinal = 2
y0 = 1
phi_np = lambda t,y: (np.sin(t))**2*y
phi_sp = lambda t,y: (sp.sin(t))**2*y
##################################################
# solution exacte
##################################################
t = sp.Symbol('t')
y = sp.Function('y')
edo = sp.Eq( sp.diff(y(t),t) , phi_sp(t,y(t)) )
solpar = sp.dsolve(edo, ics={y(t0): y0})
display(Markdown(f"**Solution exacte :** \({sp.latex(solpar)}\)"))
sol_exacte = sp.lambdify(t,solpar.rhs,'numpy')
##################################################
# schéma
##################################################
def multipas(phi, tt, sol_exacte):
h = tt[1] - tt[0]
uu = np.empty_like(tt)
# initialisation avec solution exacte pour ordre de convergence
uu[0:3] = sol_exacte(tt[0:3])
# recurrence
for i in range(2,len(tt) - 1):
uu[i+1] = a_0*uu[i] + a_1*uu[i-1] + a_2*uu[i-2] + h*(b_0*phi(tt[i],uu[i]) + b_1*phi(tt[i-1],uu[i-1]) + b_2*phi(tt[i-2],uu[i-2]))
return uu
##################################################
# ordre
##################################################
H = []
err = []
plt.figure(figsize=(16,5))
ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)
N = 10
for k in range(6):
N += 5
tt = np.linspace(t0, tfinal, N + 1)
H.append( tt[1] - tt[0] )
yy = sol_exacte(tt)
uu = multipas(phi_np, 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.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'Le schéma a ordre {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)Solution exacte : \(y{\left(t \right)} = e^{\frac{t}{2} - \frac{\sin{\left(2 t \right)}}{4}}\)
EXERCICE 3 : Étude d’un schéma RK implicite à 2 étages
Considérons le schéma RK à \(s=2\) étages dont la matrice de Butcher est la suivante :
\(
\begin{array}{c|cc}
\alpha & 0 & \alpha \\
1 & 1 & 0\\
\hline
& b_0 & b_1
\end{array}
\)
Les coefficients \(b_0\) et \(b_1\) dépendent de la première lettre de votre prénom :
\( \begin{array}{c|cc} \text{Lettre} & b_0 & b_1 \\ \hline A-F & {2}/{3} & {1}/{3} \\ G-L & {1}/{3} & {2}/{3} \\ M-R & {4}/{3} & -{1}/{3} \\ S-Z & -{1}/{3} & {4}/{3} \\ \end{array} \)
Attention : ils sont donnés comme des fractions sympy, il faudra les convertir en flottants pour les utiliser dans vos codes d’approximation numérique.
Déterminer \(\alpha\) pour que le schéma soit d’ordre maximal.
Tester empiriquement l’ordre de convergence de ce schéma sur le problème de Cauchy \(y'(t) = 4t(1-t^2)-y(t)\) et \(y(0) = 1\) pour \(t\in[0,1]\).
Étudier la A-stabilité du schéma.
Code
# ================================================
# On fixe $b_0$ et $b_1$
# ================================================
prenom = "Gloria"
if prenom[0].upper() >= 'A' and prenom[0].upper() <= 'F':
b0 = sp.Rational(2, 3)
b1 = sp.Rational(1, 3)
elif prenom[0].upper() >= 'G' and prenom[0].upper() <= 'L':
b0 = sp.Rational(1, 3)
b1 = sp.Rational(2, 3)
elif prenom[0].upper() >= 'M' and prenom[0].upper() <= 'R':
b0 = sp.Rational(4, 3)
b1 = sp.Rational(-1, 3)
elif prenom[0].upper() >= 'S' and prenom[0].upper() <= 'Z':
b0 = sp.Rational(-1, 3)
b1 = sp.Rational(4, 3)
display(Markdown(f"**{prenom}** : $b_0 = {b0}$, $b_1 = {b1}$"))Gloria : \(b_0 = 1/3\), \(b_1 = 2/3\)
C’est un schéma implicite. La barrière de Butcher affirme qu’un schéma RK à \(s=2\) étages implicite ne peut pas être d’ordre supérieur à \(2s=4\).
Calculons son ordre en fonction de \(\alpha\):
il est consistant pour tout \(\alpha\) (et tout couple \((b_0,b_1)\) donné) car
\( \begin{cases} b_0 + b_1 = 1\\ a_{0, 0} + a_{0, 1} = \alpha = c_{0}\\ a_{1, 0} + a_{1, 1} = 1 = c_{1} \end{cases} \)il est d’ordre 2 ssi \(\alpha=1-\frac{1}{2b_0}\) car
\( b_{0} c_{0} + b_{1} c_{1} = \frac{1}{2} \quad \iff \quad b_{0}\times \alpha + (1-b_0)\times 1 = \frac{1}{2} \quad \iff \quad \alpha = 1-\frac{1}{2b_0} \)il est d’ordre 3 ssi \(b_0=\frac{3}{4}\) car
\( b_{0} c_{0}^2 + b_{1} c_{1}^2 = \frac{1}{3} \quad \iff \quad b_0 \left(1-\frac{1}{2b_0}\right)^2 + (1-b_0)\times 1^2 = \frac{1}{3} \quad \iff \quad b_0 = \frac{3}{4} . \)
Conclusion : tous les schémas donnés sont d’ordre 2 si on pose \(\alpha\) comme ci-dessous :
\(
\begin{array}{c|cc|c}
\text{Lettre} & b_0 & b_1 & \alpha\\
\hline
A-F & {2}/{3} & {1}/{3} & 1/4\\
G-L & {1}/{3} & {2}/{3} & -1/2\\
M-R & {4}/{3} & -{1}/{3} & 5/8\\
S-Z & -{1}/{3} & {4}/{3} & 5/2\\
\end{array}
\)
Code
prlat = lambda *args : display(Math(''.join( sp.latex(arg) for arg in args )))
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**"))
But = matrice_Butcher(s, A, b, c)
Eqs = []
display(Markdown(f"**On a {s+1} conditions pour avoir consistance = pour être d'ordre 1**"))
Eq_1 = ordre_1(s, A, b, c, j)
Eqs.append(Eq_1)
display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
Eq_2 = ordre_2(s, A, b, c, j)
Eqs.append([Eq_2])
display(Markdown("**On doit ajouter 2 conditions pour être d'ordre 3**"))
Eq_3 = ordre_3(s, A, b, c, j)
Eqs.append(Eq_3)
# display(Markdown("**On doit ajouter 4 conditions pour être d'ordre 4**"))
# Eq_4 = ordre_4(s, A, b, c, j)
# Eqs.append(Eq_4)
return Eqs
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 But
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))
Eq_1 = [sp.Eq( sum(b).simplify(), 1 ) ]
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 ))
Eq_1.append( sp.Eq( somma, c[i] ) )
return Eq_1
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))
Eq_2 = sp.Eq( sum([b[i]*c[i] for i in range(s)]).simplify(), sp.Rational(1,2) )
return Eq_2
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))
Eq_3_1 = sp.Eq( sum([b[i]*c[i]**2 for i in range(s)]).simplify(), sp.Rational(1,3) )
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))
Eq_3_2 = sp.Eq( sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sp.Rational(1,6) )
return [Eq_3_1, Eq_3_2]
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))
Eq_4_1 = sp.Eq( sum([b[i]*c[i]**3 for i in range(s)]).simplify(), sp.Rational(1,4) )
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))
Eq_4_2 = sp.Eq( sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sp.Rational(1,8) )
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))
Eq_4_3 = sp.Eq( sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify(), sp.Rational(1,12) )
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))
Eq_4_4 = sp.Eq( 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(), sp.Rational(1,24) )
return [Eq_4_1, Eq_4_2, Eq_4_3, Eq_4_4]
##################################################
# Conditions avec s étages
##################################################
s = 2
# alpha = sp.symbols(r'\alpha')
b = [ b0, b1]
alpha= 1-1/(2*b0)
display(Math(f"\\alpha = {alpha}, b = {b}"))
c = [alpha,1]
A = [[0, alpha],
[1,0]]
Eqs = ordre_RK(s, A, b, c)
\(\displaystyle \alpha = -1/2, b = [1/3, 2/3]\)
Matrice de Butcher
\(\displaystyle \left[\begin{matrix}- \frac{1}{2} & 0 & - \frac{1}{2}\\1 & 1 & 0\\ & \frac{1}{3} & \frac{2}{3}\end{matrix}\right]\)
On a 3 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}=- \frac{1}{2}\text{ doit être égale à }- \frac{1}{2}\)
\(\displaystyle \sum_{j=1}^s a_{1j}=1\text{ doit être égale à }1\)
On doit ajouter 1 condition pour être d’ordre 2
\(\displaystyle \sum_{j=1}^s b_j c_j=\frac{1}{2}\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=\frac{3}{4}\text{ doit être égale à }\frac{1}{3}\)
\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=- \frac{1}{2}\text{ doit être égale à }\frac{1}{6}\)
Le schéma correspond à la méthode suivante :
\( \begin{cases} u_0 = y_0\\ k_0 = \varphi\Big(t_n+c_0h,u_n+h\big(a_{00}k_0+a_{01}k_1\big)\Big) \\ k_1 = \varphi\Big(t_n+c_1h,u_n+h\big(a_{10}k_0+a_{11}k_1\big)\Big) \\ u_{n+1} = u_n + h\big(b_0k_0+b_1k_1\big) \end{cases} \)
soit encore, dans notre cas particulier,
\( \begin{cases} u_0 = y_0\\ k_0 = \varphi\Big(t_n+\alpha h,u_n+\alpha h k_1\Big) \\ k_1 = \varphi\Big(t_{n+1},u_n+hk_0\Big) \\ u_{n+1} = u_n + h\big(b_0k_0+b_1k_1\big) \end{cases} \)
Code
##################################################
# EDO
##################################################
t0 = 0
tfinal = 1
y0 = 1
phi = lambda t,y: 4*t*(1-t**2)-y
##################################################
# solution exacte
##################################################
t = sp.Symbol('t')
y = sp.Function('y')
edo = sp.Eq( (y(t)).diff() , phi(t,y(t)) )
solpar = sp.dsolve(edo , ics={y(t0):y0})
display(solpar)
sol_exacte = sp.lambdify(t,solpar.rhs,'numpy')
##################################################
# schéma
##################################################
# Paramètres du schéma recuperés depuis la matrice de Butcher obtenue précédemment
cc = np.array(c)
bb = np.array(b)
AA = np.array(A)
# Schéma de Runge-Kutta implicite à s étages
def RK(phi, tt, y0):
h = tt[1] - tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for n in range(len(tt) - 1):
k_init = phi(tt[n],uu[n])
KK = fsolve( lambda kk: [ -kk[i] + phi( tt[n]+cc[i]*h , uu[n] + h*sum([kk[j]*AA[i,j] for j in range(s)]) ) for i in range(s) ], k_init*np.ones((s,1)) )
uu[n+1] = uu[n] + h*sum([KK[j]*b[j] for j in range(s)])
return uu
##################################################
# ordre
##################################################
H = []
err = []
plt.figure(figsize=(16,5))
ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)
N = 50
for k in range(6):
N += 10
tt = np.linspace(t0, tfinal, N + 1)
H.append( tt[1] - tt[0] )
yy = sol_exacte(tt)
uu = RK(phi, tt, y0)
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.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'Le schéma a ordre {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)\(\displaystyle y{\left(t \right)} = - 4 t^{3} + 12 t^{2} - 20 t + 20 - 19 e^{- t}\)
On applique le schéma au problème de Cauchy \(y'(t) = -\beta y(t)\) et \(y(0) = 1\) avec \(\beta>0\) dont la solution est \(y(t)=e^{-\beta t}\).
On obtient alors la relation de récurrence suivante :
Code
u_n = sp.symbols(r'u_n')
u_np1 = sp.symbols(r'u_{n+1}')
dt = sp.symbols(r'h', positive=True)
x = sp.symbols(r'x', positive=True)
k_0 = sp.symbols(r'k_0')
k_1 = sp.symbols(r'k_1')
kk = [k_0, k_1]
beta = sp.symbols(r'\beta', positive=True)
phi = lambda y: -beta*y
Eqs = [ sp.Eq(kk[i] , phi(u_n + dt*sum([kk[j]*A[i][j] for j in range(s)]))) for i in range(s) ]
display(Eqs)
sol = sp.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sp.Eq( u_np1 , u_n + dt*sum([KK[j]*b[j] for j in range(s)]).factor() )
display(schema)
[Eq(k_0, -\beta*(-h*k_1/2 + u_n)), Eq(k_1, -\beta*(h*k_0 + u_n))]\(\displaystyle u_{n+1} = \frac{\beta h u_{n} \left(\beta h - 2\right)}{\beta^{2} h^{2} + 2} + u_{n}\)
On peut expliciter \(u_n\) en fonction de \(u_0\) et du produit \(\beta h\) (on notera \(x\) le produit \(\beta h\)) :
Code
q = (schema.rhs/u_n).simplify().subs(beta*dt,x)
display(Math(sp.latex(sp.Eq(sp.Symbol('q(x)'),q.simplify()))))\(\displaystyle q(x) = \frac{2 \left(x^{2} - x + 1\right)}{x^{2} + 2}\)
La suite vérifie \(u_n \xrightarrow[n\to \infty]{} 0\) ssi \(\left|q(x)\right|<1\). En étudiant brièvement la fonction \(q(x)\), on trouve que
\( \begin{array}{c|cc|c} \text{Lettre} & b_0 & b_1 & \text{A-stabilité} \\ \hline A-F & {2}/{3} & {1}/{3} & \forall h \\ G-L & {1}/{3} & {2}/{3} & \beta h <2 \\ M-R & {4}/{3} & -{1}/{3} & h<\frac{2(\sqrt{7}-1)}{3}\\ S-Z & -{1}/{3} & {4}/{3} & h<\frac{\sqrt{37}-1}{9}\\ \end{array} \)
Code
sp.plot(q, (x,0,10), ylim=(-1,1),
ylabel='q(x)', xlabel='x', title='q(x) en fonction de x')
sp.solve(abs(q)<1, x)
\(\displaystyle x < 2\)