⏱ 2h30, deposer le notebook dans Moodle
Considérons le problème de Cauchy suivant
\( \begin{cases}
y'(t) = \lambda \Big(y(t) - \cos(t)\Big) - \sin(t) &\text{pour } t\in[0,10], \\
y(0) = 1+\varepsilon.
\end{cases} \)
On peut résoudre ce problème de Cauchy à la main en remarquant que l’edo se réécrit comme \( \lambda \big(y(t) - \cos(t)\big) = y'(t) +\sin(t) = \big(y(t) - \cos(t)\big)'. \) Si on pose \(z(t)=y(t) - \cos(t)\) on a le problème de Cauchy \(z'(t) = \lambda z(t)\) avec \(z(0) = y(0) - 1\). La solution est donc \(z(t) = (y(0) - 1) e^{\lambda t}\) et donc \( y(t) = \varepsilon e^{\lambda t} + \cos(t). \)
On peut vérifier nos calculs en utilisant sympy :
%reset -f
import sympy as sp
t = sp.symbols('t')
y = sp.Function('y')(t)
lamb = sp.symbols('\lambda', real=True)
epsilon = sp.symbols('\epsilon', real=True)
dydt = y.diff(t)
ode = sp.Eq(dydt, lamb*(y-sp.cos(t))-sp.sin(t))
sol = sp.dsolve(ode, y, ics={y.subs(t, 0): 1+epsilon})
sol
\(\displaystyle y{\left(t \right)} = \epsilon e^{\lambda t} + \cos{\left(t \right)}\)
import matplotlib.pyplot as plt
import numpy as np
tt = np.linspace(0, 10, 100)
fig, ax = plt.subplots(figsize=(8, 6))
for epsi in [-0.1,0,0.1]:
y_sol = sp.lambdify(t, sol.rhs.subs({lamb:10,epsilon:epsi}), 'numpy')
plt.plot(tt, y_sol(tt), label=fr'$\epsilon={epsi}$')
ax.set_xlabel('t')
ax.set_ylabel('y')
ax.grid()
ax.legend()
ax.axes.set_ylim(-3, 3);
import matplotlib.pyplot as plt
import numpy as np
tt = np.linspace(0, 0.5, 100)
fig, ax = plt.subplots(figsize=(8, 6))
for epsi in [-0.1,0,0.1]:
y_sol = sp.lambdify(t, sol.rhs.subs({lamb:-10,epsilon:epsi}), 'numpy')
plt.plot(tt, y_sol(tt), label=fr'$\epsilon={epsi}$')
ax.set_xlabel('t')
ax.set_ylabel('y')
ax.grid()
ax.legend()
ax.axes.set_ylim(1-0.2, 1+0.2);
La matrice de Butcher d’un schéma RK à \(s=4\) étages quelconque s’écrit :
\(
\begin{array}{c|cccc}
c_0 & a_{00} & a_{01} & a_{02} & a_{03}\\
c_1 & a_{10} & a_{11} & a_{12} & a_{13}\\
c_2 & a_{20} & a_{21} & a_{22} & a_{23}\\
c_3 & a_{30} & a_{31} & a_{32} & a_{33}\\
\hline
& b_0 & b_1 & b_2 & b_3
\end{array}
\)
Elle contient 24 coefficients.
Écrire la matrice de Butcher d’un schéma RK à \(s=4\) étages explicite consistant (elle ne comportera que les 9 coefficients \(c_1\), \(c_2\), \(c_3\), \(b_1\), \(b_2\), \(b_3\), \(a_{21}\), \(a_{31}\), \(a_{32}\)).
Sans effectuer de calculs, donner l’ordre maximal d’un schéma explicite à 4 étages.
Soit \(\kappa\) le nombre de lettres de votre prénom, dans mon cas je pose kappa = len("Gloria").
Poser \(c_1=\frac{1}{\kappa}\), \(c_2=1-c_1\) et \(c_3=1\). Déterminer les autres 6 paramètres pour que le schéma explicite soit d’ordre maximal.
Tester empiriquement l’ordre de convergence de ce schéma sur le problème de Cauchy suivant : \(y'(t) = 10(y(t)-\sin(t))+\cos(t)\) et \(y(0) = 0\) pour \(t\in[0,\pi/2]\). Est-ce un schéma adapté à ce type de problème ?
Déterminer sous quelle condition sur \(h\) ce schéma est A-stable.
%reset -f
import sympy
sympy.init_printing()
from IPython.display import display, Math, Markdown
prlat = lambda *args : display(Math(''.join( sympy.latex(arg) for arg in args )))
import numpy as np
import matplotlib.pyplot as plt
# from scipy.optimize import fsolve
Si le schéma est explicite, la matrice \(\mathbb A\) est triangulaire inférieure et la diagonale ne contient que des zéros, donc \(a_{00}=a_{01}=a_{02}=a_{03}=a_{11}=a_{12}=a_{13}=a_{22}=a_{23}=a_{33}=0\).
Si le schéma est consistant alors
\( \begin{cases}
\displaystyle\sum_{j=0}^3 b_{i}=1& \\
\displaystyle c_i=\sum_{j=0}^3 a_{ij}&i=0,\dots,3
\end{cases} \)
ainsi \(c_0=0\), \(a_{10}=c_1\), \(a_{20}=c_2-a_{21}\), \(a_{30}=c_3-a_{31}-a_{32}\) et enfin \(b_0=1-b_1-b_2-b_3\).
Avec ces 15 relations, il ne reste plus que 24-15=9 coefficients à déterminer, à savoir \(a_{21}\), \(a_{31}\), \(a_{32}\), \(b_1\), \(b_2\), \(b_3\), \(c_1\), \(c_2\) et \(c_3\):
\(
\begin{array}{c|cccc}
0 & 0 & 0 & 0 & 0\\
c_1 & c_1 & 0& 0 & 0\\
c_2 & c_2-a_{21} & a_{21}& 0 & 0\\
c_3 & c_3-a_{31}-a_{32} & a_{31}& a_{32} & 0\\
\hline
& 1-b_1-b_2-b_3 & b_1 & b_2 & b_3
\end{array}
\)
La barrière de Butcher affirme qu’un schéma RK à \(s=4\) étages explicite ne peut pas être d’ordre supérieur à \(s=4\).
On a déjà pris en compte les conditions de consistance. Il reste à déterminer les paramètres pour que le schéma soit d’ordre maximal. Cela revient à résoudre le système d’équations non-linéaires suivant :
\( \begin{cases} b_{1} c_{1} + b_{2} c_{2} + b_{3} c_{3} = \frac{1}{2}, \\ b_{1} c_{1}^{2} + b_{2} c_{2}^{2} + b_{3} c_{3}^{2} = \frac{1}{3}, \\ a_{21} b_{2} c_{1} + a_{31} b_{3} c_{1} + a_{32} b_{3} c_{2} = \frac{1}{6}, \\ b_{1} c_{1}^{3} + b_{2} c_{2}^{3} + b_{3} c_{3}^{3} = \frac{1}{4}, \\ a_{21} b_{2} c_{1} c_{2} + a_{31} b_{3} c_{1} c_{3} + a_{32} b_{3} c_{2} c_{3} = \frac{1}{8}, \\ a_{21} b_{2} c_{1}^{2} + a_{31} b_{3} c_{1}^{2} + a_{32} b_{3} c_{2}^{2} = \frac{1}{12}, \\ a_{21} a_{32} b_{3} c_{1} = \frac{1}{24} \end{cases} \)
On peut aussi utiliser sympy pour retrouver ce système :
def ordre_RK(s, A=None, b=None, c=None):
j = sympy.symbols('j')
if A is None:
A = sympy.MatrixSymbol('a', s, s)
else:
A = sympy.Matrix(A)
if c is None:
c = sympy.symbols(f'c_0:{s}')
if b is None:
b = sympy.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**"))
# ordre_1(s, A, b, c, j)
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.extend(Eq_3)
display(Markdown("**On doit ajouter 4 conditions pour être d'ordre 4**"))
Eq_4 = ordre_4(s, A, b, c, j)
Eqs.extend(Eq_4)
return Eqs
def matrice_Butcher(s, A, b, c):
But = sympy.Matrix(A)
But = But.col_insert(0, sympy.Matrix(c))
last = [sympy.Symbol(" ")]
last.extend(b)
But = But.row_insert(s,sympy.Matrix(last).T)
display(Math(sympy.latex(sympy.Matrix(But))))
return But
def ordre_1(s, A, b, c, j):
texte = "\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 = sympy.summation(A[i,j],(j,0,s-1)).simplify()
texte = r'\sum_{j=1}^s a_{'+str(i)+r'j}=' + sympy.latex( somma )
texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
display(Math( texte ))
return None
def ordre_2(s, A, b, c, j):
texte = '\sum_{j=1}^s b_j c_j=' +sympy.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 = sympy.Eq( sum([b[i]*c[i] for i in range(s)]).simplify(), sympy.Rational(1,2) )
return Eq_2
def ordre_3(s, A, b, c, j):
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))
Eq_3_1 = sympy.Eq( sum([b[i]*c[i]**2 for i in range(s)]).simplify(), sympy.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 + sympy.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{6}"
display(Math(texte))
Eq_3_2 = sympy.Eq( sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,6) )
return [Eq_3_1, Eq_3_2]
def ordre_4(s, A, b, c, j):
texte = '\sum_{j=1}^s b_j c_j^3='
texte += sympy.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 = sympy.Eq( sum([b[i]*c[i]**3 for i in range(s)]).simplify(), sympy.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 + sympy.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{8}"
display(Math(texte))
Eq_4_2 = sympy.Eq( sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sympy.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 + sympy.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{12}"
display(Math(texte))
Eq_4_3 = sympy.Eq( sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify(), sympy.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 + sympy.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{24}"
display(Math(texte))
Eq_4_4 = sympy.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(), sympy.Rational(1,24) )
return [Eq_4_1, Eq_4_2, Eq_4_3, Eq_4_4]
##################################################
# Conditions avec s étages
##################################################
s = 4
# Construction de la matrice de Butcher en supposant
# que le schéma est Explicite
# et qu'il est consistant (i.e. au moins d'ordre 1)
b_1, b_2, b_3 = sympy.symbols('b_1 b_2 b_3', real=True)
c_1, c_2, c_3 = sympy.symbols('c_1 c_2 c_3', real=True)
a_21, a_31, a_32 = sympy.symbols('a_21 a_31 a_32', real=True)
b = [1-b_1-b_2-b_3, b_1, b_2, b_3]
c = [0, c_1, c_2, c_3]
A = [[0, 0, 0, 0],
[c_1, 0, 0, 0],
[c_2-a_21, a_21, 0, 0],
[c_3-a_31-a_32, a_31, a_32, 0]]
Eqs = ordre_RK(s, A, b, c)
# for eq in Eqs:
# display(eq)Matrice de Butcher
\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0 & 0\\c_{1} & c_{1} & 0 & 0 & 0\\c_{2} & - a_{21} + c_{2} & a_{21} & 0 & 0\\c_{3} & - a_{31} - a_{32} + c_{3} & a_{31} & a_{32} & 0\\ & - b_{1} - b_{2} - b_{3} + 1 & b_{1} & b_{2} & b_{3}\end{matrix}\right]\)
On doit ajouter 1 condition pour être d’ordre 2
\(\displaystyle \sum_{j=1}^s b_j c_j=b_{1} c_{1} + b_{2} c_{2} + b_{3} c_{3}\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=b_{1} c_{1}^{2} + b_{2} c_{2}^{2} + b_{3} c_{3}^{2}\text{ doit être égale à }\frac{1}{3}\)
\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=a_{21} b_{2} c_{1} + a_{31} b_{3} c_{1} + a_{32} b_{3} c_{2}\text{ doit être égale à }\frac{1}{6}\)
On doit ajouter 4 conditions pour être d’ordre 4
\(\displaystyle \sum_{j=1}^s b_j c_j^3=b_{1} c_{1}^{3} + b_{2} c_{2}^{3} + b_{3} c_{3}^{3}\text{ doit être égale à }\frac{1}{4}\)
\(\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=a_{21} b_{2} c_{1} c_{2} + a_{31} b_{3} c_{1} c_{3} + a_{32} b_{3} c_{2} c_{3}\text{ doit être égale à }\frac{1}{8}\)
\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=a_{21} b_{2} c_{1}^{2} + a_{31} b_{3} c_{1}^{2} + a_{32} b_{3} c_{2}^{2}\text{ doit être égale à }\frac{1}{12}\)
\(\displaystyle \sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k=a_{21} a_{32} b_{3} c_{1}\text{ doit être égale à }\frac{1}{24}\)
On résout ce système NON LINÉAIRE pour trouver les coefficients dans le cas particulier où \(c_1\), \(c_2\) et \(c_3\) sont donnés. Voici le démarche de résolution :
on trouve \(b_1\), \(b_2\) et \(b_3\) en résolvant le sous-sytème (linéaire) :
\(
\begin{cases}
b_{1} c_{1} + b_{2} c_{2} + b_{3} c_{3} = \frac{1}{2}, \\
b_{1} c_{1}^{2} + b_{2} c_{2}^{2} + b_{3} c_{3}^{2} = \frac{1}{3}, \\
b_{1} c_{1}^{3} + b_{2} c_{2}^{3} + b_{3} c_{3}^{3} = \frac{1}{4}, \\
\end{cases}
\)
on trouve \(a_{21}\), \(a_{31}\) et \(a_{32}\) en résolvant le sous-sytème (linéaire) :
\(
\begin{cases}
a_{21} b_{2} c_{1} + a_{31} b_{3} c_{1} + a_{32} b_{3} c_{2} = \frac{1}{6}, \\
a_{21} b_{2} c_{1} c_{2} + a_{31} b_{3} c_{1} c_{3} + a_{32} b_{3} c_{2} c_{3} = \frac{1}{8}, \\
a_{21} b_{2} c_{1}^{2} + a_{31} b_{3} c_{1}^{2} + a_{32} b_{3} c_{2}^{2} = \frac{1}{12}, \\
\end{cases}
\)
on trouve le dernier paramètre \(a_{32}\) en résolvant l’équation :
\(
a_{21} a_{32} b_{3} c_{1} = \frac{1}{24}.
\)
Remarque Dans l’énoncé j’ai donné \(c_3=1\) car cela simplifie les calculs à la main, mais en réalité il suffit de fixer \(c_1\) et \(c_2\) et on retrouve \(c_3=1\) en résolvant le système :
# 7 équations (non linéaires) et 9 inconnues
# On fixe 2 paramètres et on résoud le système non linéaire
kappa = len("Gloria")
choix = {c_1 : sympy.Rational(1,kappa), c_2 : sympy.Rational(kappa-1,kappa)}
display(Markdown("On fixe $c_1$ et $c_2$"))
display(choix)
sol = sympy.solve( [eq.subs(choix) for eq in Eqs] )[0]
display(Markdown("En imposant d'être d'ordre 4, on trouve les 7 autres paramètres :"))
display(sol)
display(Markdown("Autrement dit, la matrice de Butcher d'un schéma explicite consistant à 4 étages"))
But = matrice_Butcher(s, A, b, c)
display(Markdown("devient, avec les valeurs de $c_1$ et $c_2$ données en imposant l'ordre 4"))
But = But.subs(choix).subs(sol)
display(But) On fixe \(c_1\) et \(c_2\)
\(\displaystyle \left\{ c_{1} : \frac{1}{6}, \ c_{2} : \frac{5}{6}\right\}\)
En imposant d’être d’ordre 4, on trouve les 7 autres paramètres :
\(\displaystyle \left\{ a_{21} : \frac{5}{2}, \ a_{31} : 10, \ a_{32} : -1, \ b_{1} : \frac{3}{5}, \ b_{2} : \frac{3}{5}, \ b_{3} : - \frac{1}{10}, \ c_{3} : 1\right\}\)
Autrement dit, la matrice de Butcher d’un schéma explicite consistant à 4 étages
\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0 & 0\\c_{1} & c_{1} & 0 & 0 & 0\\c_{2} & - a_{21} + c_{2} & a_{21} & 0 & 0\\c_{3} & - a_{31} - a_{32} + c_{3} & a_{31} & a_{32} & 0\\ & - b_{1} - b_{2} - b_{3} + 1 & b_{1} & b_{2} & b_{3}\end{matrix}\right]\)
devient, avec les valeurs de \(c_1\) et \(c_2\) données en imposant l’ordre 4
\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0 & 0\\\frac{1}{6} & \frac{1}{6} & 0 & 0 & 0\\\frac{5}{6} & - \frac{5}{3} & \frac{5}{2} & 0 & 0\\1 & -8 & 10 & -1 & 0\\ & - \frac{1}{10} & \frac{3}{5} & \frac{3}{5} & - \frac{1}{10}\end{matrix}\right]\)
##################################################
# EDO
##################################################
t0 = 0
y0 = 0
tfinal = np.pi/2
# Prothero-Robinson problem (stiff ou numériquement mal posé selon lambda)
phi_sym = lambda t,y: sympy.cos(t)-10*(y-sympy.sin(t))
phi = lambda t,y: np.cos(t)-10*(y-np.sin(t))
##################################################
# solution exacte
##################################################
t = sympy.Symbol('t')
y = sympy.Function('y')
edo = sympy.Eq( (y(t)).diff() , phi_sym(t,y(t)) )
solpar = sympy.dsolve(edo , ics={y(t0):y0})
display(solpar)
sol_exacte = sympy.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(But[:-1,0]).astype(float)
bb = np.array(list(But[-1,1:])).astype(float)
AA = np.array(But[:-1,1:]).astype(float)
# Schéma de Runge-Kutta explicite à 4 étages
def RK(phi, tt, y0):
h = tt[1] - tt[0]
uu = np.zeros_like(tt)
# initialisation
uu[0] = y0
# recurrence
for n in range(len(tt) - 1):
k0 = phi(tt[n] , uu[n] )
k1 = phi(tt[n] + cc[1,0]*h, uu[n] + h* AA[1,0]*k0)
k2 = phi(tt[n] + cc[2,0]*h, uu[n] + h*(AA[2,0]*k0 + AA[2,1]*k1))
k3 = phi(tt[n] + cc[3,0]*h, uu[n] + h*(AA[3,0]*k0 + AA[3,1]*k1 + AA[3,2]*k2))
uu[n+1] = uu[n] + h*(bb[0]*k0 + bb[1]*k1 + bb[2]*k2 + bb[3]*k3)
return uu
##################################################
# ordre
##################################################
H = []
err = []
plt.figure(figsize=(16,5))
ax1 = plt.subplot(1,2,1)
ax1.set_xlabel('$t$')
ax1.set_ylabel('$y$')
ax1.grid(True)
ax2 = plt.subplot(1,2,2)
ax2.set_xlabel(r'$\ln(h)$')
ax2.set_ylabel(r'$\ln(e)$')
ax2.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.plot( np.log(H), np.log(err), 'r-o')
ax2.set_title(f'Le schéma a ordre {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}');\(\displaystyle y{\left(t \right)} = \sin{\left(t \right)}\)
c = list(But[:-1,0])
b = list(But[-1,1:])
A = But[:-1,1:]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 :
u_n = sympy.symbols(r'u_n')
u_np1 = sympy.symbols(r'u_{n+1}')
dt = sympy.symbols(r'h', positive=True)
x = sympy.symbols(r'x', positive=True)
beta = sympy.symbols(r'\beta', positive=True)
phi_A = lambda y : -beta*y
k0 = phi_A(u_n)
k1 = phi_A(u_n + dt* A[1,0]*k0)
k2 = phi_A(u_n + dt*(A[2,0]*k0 + A[2,1]*k1))
k3 = phi_A(u_n + dt*(A[3,0]*k0 + A[3,1]*k1 + A[3,2]*k2))
for i,k in enumerate([k0,k1,k2,k3]):
display(Math(sympy.latex(sympy.Eq(sympy.Symbol(f'k_{i}'), k.factor(beta*dt)))))
u_np1 = u_n + dt*(b[0]*k0 + b[1]*k1 + b[2]*k2 + b[3]*k3)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('u_{n+1}'),u_np1.simplify()))))
\(\displaystyle k_{0} = - \beta u_{n}\)
\(\displaystyle k_{1} = - \frac{\beta u_{n} \left(- \beta h + 6\right)}{6}\)
\(\displaystyle k_{2} = - \frac{\beta u_{n} \left(5 \beta^{2} h^{2} - 10 \beta h + 12\right)}{12}\)
\(\displaystyle k_{3} = - \frac{\beta u_{n} \left(5 \beta^{3} h^{3} + 10 \beta^{2} h^{2} - 12 \beta h + 12\right)}{12}\)
\(\displaystyle u_{n+1} = \frac{u_{n} \left(\beta^{4} h^{4} - 4 \beta^{3} h^{3} + 12 \beta^{2} h^{2} - 24 \beta h + 24\right)}{24}\)
# Variante "semi-automatique"
# k_0 = sympy.symbols(r'k_0')
# k_1 = sympy.symbols(r'k_1')
# k_2 = sympy.symbols(r'k_2')
# k_3 = sympy.symbols(r'k_3')
# kk = [k_0, k_1, k_2, k_3]
# phi_A = lambda y: -beta*y
# Eqs = [ sympy.Eq( kk[i] , phi_A(u_n + dt*sum([kk[j]*A[i,j] for j in range(s)])) ) for i in range(s) ]
# for eq in Eqs:
# display(eq)
# sol = sympy.solve( Eqs, kk ) # c'est un dictionnaire
# KK = [sol[kk[i]] for i in range(s)]
# schema = sympy.Eq( u_np1 , u_n + dt*sum([KK[j]*b[j] for j in range(s)]) )
# # display(schema)
# display(Math(sympy.latex(sympy.Eq(sympy.Symbol('u_{n+1}'),schema.rhs.simplify()))))
On a donc une suite du type \(u_{n+1}=q(\beta h) u_n\). On peut expliciter \(u_n\) en fonction de \(u_0\) et du produit \(\beta h\) (on notera \(x\) le produit \(\beta h\)) car \(u_n=q(x)^nu_0\) avec \(q(x)=u_{n+1}/u_n\) :
q = (u_np1/u_n).simplify().subs(beta*dt,x)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('q(x)'),q.simplify()))))\(\displaystyle q(x) = \frac{x^{4}}{24} - \frac{x^{3}}{6} + \frac{x^{2}}{2} - x + 1\)
La suite vérifie \(\lim_{n\to \infty} u_n = 0\) ssi \(\left|q(x)\right|<1\). En étudiant brièvement la fonction \(q(x)\), on trouve que \(q(x)<1\) ssi \(x<x_0\) donné ci-dessous:.
sympy.plot(q, (x,0,4), ylim=(-1,1),
ylabel='q(x)', xlabel='x', title='q(x) en fonction de x')
sol = sympy.solveset(q-1, x,domain=sympy.S.Reals)
sol0, x0 = sol.args
display(Markdown(f"$\\beta h<{x0.evalf()}$"))
\(\beta h<2.78529356340528\)
Écrire le schéma multipas à \(q=2\) pas (comportant 5 coefficients).
Sans effectuer de calculs, répondre aux deux questions suivantes :
Parmi ces schémas, calculer les coefficients du seul schéma consistant d’ordre 4. Montrer que le schéma ainsi trouvé est zéro-stable.
Tester empiriquement l’ordre de convergence du schéma ainsi obtenu sur le problème de Cauchy \(y'(t) = 4t(t^2-1)+y(t)\) et \(y(0) = 1\) pour \(t\in[0,1]\).
Considérons les schémas explicites à 2 pas (comportant 4 coefficients).
%reset -f
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
from IPython.display import display, Math, MarkdownUn schéma implicite à \(q=2\) étapes s’écrit
\( u_{n+1}=a_0u_n+a_1u_{n-1} +h\big(b_{-1}\varphi_{n+1}+b_0\varphi_n+b_1\varphi_{n-1}\big) \)
La première barrière de Dahlquist affirme qu’un schéma à \(q=2\) pas
Il contient les 5 paramètres \(a_0\), \(a_1\), \(b_{-1}\), \(b_0\), \(b_1\), il est donc théoriquement possible de construire un schéma d’ordre 4 en résolvant le système linéaire de \(q+1=5\) équations suivant
\( \begin{cases} \xi(0) = a_{0} + a_{1}=1 \\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1}=1 \\ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right)=1 \\ \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right)=1 \\ \xi(4) = a_{1} - 4 \left(b_{1} - b_{-1}\right)=1 \end{cases} \)
On trouve \(a_0=0\), \(a_1=1\), \(b_{-1}=\frac{1}{3}\), \(b_0=\frac{4}{3}\) et \(b_1=\frac{1}{3}\).
Ci-dessous la résolution avec sympy :
q = 2
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()
# Définition des symboles
a_0, a_1 = sp.symbols('a_0 a_1 ', real=True)
b_0, b_1 = sp.symbols('b_0 b_1 ', real=True)
b_m1 = sp.symbols('b_{-1}', real=True)
aa = [a_0, a_1]
bb = [b_0, b_1]
systeme = [sp.Eq(xi(i, q),1) for i in range(2*q+1)]
display(Markdown(r"$\begin{cases}"+' '.join([sp.latex(s)+r"\\" for s in systeme]) + r"\end{cases}$"))
sol = sp.solve(systeme)
display(sol)\(\begin{cases}a_{0} + a_{1} = 1\\ - a_{1} + b_{0} + b_{1} + b_{-1} = 1\\ a_{1} - 2 \left(b_{1} - b_{-1}\right) = 1\\ - a_{1} + 3 \left(b_{1} + b_{-1}\right) = 1\\ a_{1} - 4 \left(b_{1} - b_{-1}\right) = 1\\\end{cases}\)
\(\displaystyle \left\{ a_{0} : 0, \ a_{1} : 1, \ b_{0} : \frac{4}{3}, \ b_{1} : \frac{1}{3}, \ b_{-1} : \frac{1}{3}\right\}\)
Le schéma s’écrit donc :
n = sp.Symbol('n')
txt = r"u_{n+1} = "
txt += sp.latex( sum( aa[j].subs(sol) * sp.Symbol(f"u_{{{n-j}}}") for j in range(q)) )
txt += "+ h \left("
txt += sp.latex( sum( bb[j].subs(sol) * sp.Symbol(f"\\varphi_{{{n-j}}}") for j in range(q)) )
txt += r"+"
txt += sp.latex( b_m1.subs(sol) * sp.Symbol(r"\varphi_{n+1}") )
txt += r"\right)"
display(Math( txt ))\(\displaystyle u_{n+1} = u_{n - 1}+ h \left(\frac{\varphi_{n - 1}}{3} + \frac{4 \varphi_{n}}{3}+\frac{\varphi_{n+1}}{3}\right)\)
Vérifions si le schéma ainsi obtenu est zéro-stable en calculant les racines du premier polynôme caractéristique qui est défini par
\( \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-1=(r-1)(r+1) \)
La méthode est zéro-stable car ses racines sont de module égal à 1 mais simples.
Remarque : on reconnaît le schéma de Milne Simpson à deux pas.
Voici le code pour vérifier la zéro-stabilité avec sympy:
r = sp.Symbol('r')
p = q-1
rho_gen = sp.poly(r**(p+1) - sum( aa[j] * r**(p-j) for j in range(p+1)), r)
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 :"))
display(Math(sp.latex(rho_gen.as_expr() ,order='old') + "=" + sp.latex(rho.as_expr() ,order='old')))
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 - a_{1} - r a_{0} + r^{2}=-1 + r^{2}\)
Ses racines sont :
\(\displaystyle \left[ -1, \ 1\right]\)
dont les modules valent :
\(\displaystyle \left[ 1.0, \ 1.0\right]\)
# AUTOMATISATION : conversion 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():
if str(cle) == "b_{-1}":
cle = "b_m1"
globals()[cle] = valeur
else :
globals()[str(cle)] = valeur
# TEST de l'affectation
for var in ['a_0', 'a_1', 'b_0', 'b_1', 'b_m1']:
print(f"{var} = {globals()[var]}")
a_0 = 0.0
a_1 = 1.0
b_0 = 1.3333333333333333
b_1 = 0.3333333333333333
b_m1 = 0.3333333333333333##################################################
# EDO
##################################################
t0 = 0
tfinal = 1
y0 = 1
phi = lambda t,y: 4*t*(t**2-1)+y
##################################################
# solution exacte
##################################################
t = sp.Symbol('t')
y = sp.Function('y')
edo = sp.Eq( sp.diff(y(t),t) , phi(t,y(t)) )
solpar = sp.dsolve(edo, ics={y(t0): y0})
display(solpar)
sol_exacte = sp.lambdify(t,solpar.rhs,'numpy')
##################################################
# schéma
##################################################
def multipas(phi, tt, sol_exacte):
h = tt[1] - tt[0]
uu = np.zeros_like(tt)
# initialisation avec solution exacte pour ordre de convergence
uu[0:q] = sol_exacte(tt[0:q])
# recurrence
for i in range(q-1,len(tt) - 1):
eq = lambda x : -x + a_0*uu[i] + a_1*uu[i-1] + h*(b_0*phi(tt[i],uu[i]) + b_1*phi(tt[i-1],uu[i-1]) + b_m1*phi(tt[i+1],x) )
tmp = fsolve ( eq , uu[i] )
uu[i+1] = tmp[0]
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)
ax2 = plt.subplot(1,2,2)
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
ax2.grid(True)
N = 10
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 = 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.legend()
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.title(f'Le schéma a ordre {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}');\(\displaystyle y{\left(t \right)} = - 4 t^{3} - 12 t^{2} - 20 t + 21 e^{t} - 20\)
Il contient les 4 paramètres \(a_0\), \(a_1\), \(b_0\), \(b_1\), il est donc théoriquement possible de construire un schéma d’ordre 3 en résolvant le système linéaire de \(4\) équations suivant
\( \begin{cases} \xi(0) = a_{0} + a_{1}=1 \\ \xi(1) = - a_{1} + b_{0} + b_{1} =1 \\ \xi(2) = a_{1} - 2 b_{1} =1 \\ \xi(3) = - a_{1} + 3 b_{1} =1 \\ \end{cases} \)
On trouve \(a_0=-4\), \(a_1=5\), \(b_0=4\) et \(b_1=2\).
Ci-dessous la résolution avec sympy :
# Définition des symboles
a_0, a_1 = sp.symbols('a_0 a_1 ', real=True)
b_0, b_1 = sp.symbols('b_0 b_1 ', real=True)
aa = [a_0, a_1]
bb = [b_0, b_1]
b_m1 = 0
q = 2
systeme = [sp.Eq(xi(i, q),1) for i in range(4)]
display(Markdown(r"$\begin{cases}"+' '.join([sp.latex(s)+r"\\" for s in systeme]) + r"\end{cases}$"))
sol = sp.solve(systeme)
display(sol)\(\begin{cases}a_{0} + a_{1} = 1\\ - a_{1} + b_{0} + b_{1} = 1\\ a_{1} - 2 b_{1} = 1\\ - a_{1} + 3 b_{1} = 1\\\end{cases}\)
\(\displaystyle \left\{ a_{0} : -4, \ a_{1} : 5, \ b_{0} : 4, \ b_{1} : 2\right\}\)
Le schéma s’écrit donc :
n = sp.Symbol('n')
txt = r"u_{n+1} = "
txt += sp.latex( sum( aa[j].subs(sol) * sp.Symbol(f"u_{{{n-j}}}") for j in range(q)) )
txt += r" + h \left("
txt += sp.latex( b_m1 * sp.Symbol(r"\varphi_{n+1}") + sum( bb[j].subs(sol) * sp.Symbol(f"\\varphi_{{{n-j}}}") for j in range(q)) )
txt += r"\right)"
display(Math(txt))\(\displaystyle u_{n+1} = 5 u_{n - 1} - 4 u_{n} + h \left(2 \varphi_{n - 1} + 4 \varphi_{n}\right)\)
Cépendant on sait, d’après la première barrière de Dalquist, qu’un schéma multipas linéaire explicite à \(q\) pas est au mieux d’ordre \(q\). Donc un schéma à 2 pas ne peut pas être d’ordre 3. 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^2+4r-5=(r-1)(r+5 .) \)
La méthode n’est pas zéro-stable car une des racines n’est pas dans le cercle unitaire.
Voici le code pour vérifier la zéro-stabilité avec sympy:
r = sp.Symbol('r')
p = q-1
rho_gen = sp.poly(r**(p+1) - sum( aa[j] * r**(p-j) for j in range(p+1)), r)
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 :"))
display(Math(sp.latex(rho_gen.as_expr() ,order='old') + "=" + sp.latex(rho.as_expr() ,order='old')))
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 - a_{1} - r a_{0} + r^{2}=-5 + 4 r + r^{2}\)
Ses racines sont :
\(\displaystyle \left[ -5, \ 1\right]\)
dont les modules valent :
\(\displaystyle \left[ 5.0, \ 1.0\right]\)
Pour trouver un schéma explicite d’ordre 2, on résout le système linéaire de 3 équations suivant qui contient donc un paramètre libre (ici \(a_0\)) :
\( \begin{cases} \xi(0) = a_{0} + a_{1}=1 \\ \xi(1) = - a_{1} + b_{0} + b_{1} =1 \\ \xi(2) = a_{1} - 2 b_{1} =1 \\ \end{cases} \)
On trouve \(a_1=1-a_0\), \(b_0=2-\frac{a_0}{2}\) et \(b_1=-\frac{a_0}{2}\).
Ci-dessous la résolution avec sympy :
systeme = [sp.Eq(xi(i, q),1) for i in range(3)]
display(Markdown(r"$\begin{cases}"+' '.join([sp.latex(s)+r"\\" for s in systeme]) + r"\end{cases}$"))
sol = sp.solve(systeme, (a_1, b_0, b_1))
display(sol)\(\begin{cases}a_{0} + a_{1} = 1\\ - a_{1} + b_{0} + b_{1} = 1\\ a_{1} - 2 b_{1} = 1\\\end{cases}\)
\(\displaystyle \left\{ a_{1} : 1 - a_{0}, \ b_{0} : 2 - \frac{a_{0}}{2}, \ b_{1} : - \frac{a_{0}}{2}\right\}\)
Imposons maintenant la zéro-stabilite. Le premier polynôme caractéristique est
\( \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}= r^{2} - r a_{0} - a_{1} = r^{2} - r a_{0} + (a_{0}-1) \)
Les racines sont \(1\) et \(a_0-1\). La méthode est zéro-stable ssi \(-1\le a_0-1<1\), c’est-à-dire \(\boxed{0\leq a_0 < 2}\).
Voici le code pour vérifier la zéro-stabilité avec sympy:
r = sp.Symbol('r')
p = q-1
rho_gen = sp.poly(r**(p+1) - sum( aa[j] * r**(p-j) for j in range(p+1)), r)
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 :"))
display(Math(sp.latex(rho_gen.as_expr() ,order='old') + "=" + sp.latex(rho.as_expr() ,order='old')))
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 - a_{1} - r a_{0} + r^{2}=-1 + a_{0} - r a_{0} + r^{2}\)
Ses racines sont :
\(\displaystyle \left[ 1, \ a_{0} - 1\right]\)
dont les modules valent :
\(\displaystyle \left[ 1.0, \ \left|{a_{0} - 1}\right|\right]\)
Le schéma s’écrit donc :
n = sp.Symbol('n')
txt = r"u_{n+1} = "
txt += sp.latex( sum( aa[j].subs(sol) * sp.Symbol(f"u_{{{n-j}}}") for j in range(q)) )
txt += r" + h \left("
txt += sp.latex( b_m1 * sp.Symbol(r"\varphi_{n+1}") + sum( bb[j].subs(sol) * sp.Symbol(f"\\varphi_{{{n-j}}}") for j in range(q)) )
txt += r"\right)"
display(Math(txt))\(\displaystyle u_{n+1} = a_{0} u_{n} + u_{n - 1} \left(1 - a_{0}\right) + h \left(- \frac{\varphi_{n - 1} a_{0}}{2} + \varphi_{n} \left(2 - \frac{a_{0}}{2}\right)\right)\)