⏱ 2h, déposer le notebook dans Moodle (si 1/3 temps supplémentaire, 2h40)
📚 Autorisés : notebooks de cours, notes personnelles.
Compléter ici :
On considère le schéma de Runge-Kutta explicite à 3 étages défini par le tableau de Butcher suivant :
\( \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & 0 & 0 \\ 1 & 0 & 1 & 0 \\ \hline & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{array} \)
Correction
C’est une méthode de Runge-Kutta explicite à \(s=3\) étages. Elle s’écrit comme
\( \begin{cases} k_1 = \varphi(t_n, u_n) \\ k_2 = \varphi(t_n+\frac{h}{2}, u_n+\frac{h}{2}k_1) \\ k_3 = \varphi(t_n+h, u_n+h k_2) \\ u_{n+1} = u_n + \frac{h}{6}(k_1+4k_2+k_3) \end{cases} \)
# IMPLEMENTATION
# =============================================================================
import numpy as np
import matplotlib.pyplot as plt
def RK(phi,tt,y0):
h = tt[1]-tt[0]
uu = np.zeros(len(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]+h/2*k1)
k3 = phi(tt[n]+h,uu[n]+h*k2)
uu[n+1] = uu[n] + h/6*(k1+4*k2+k3)
return uuLa barrière de Butcher donne une borne supérieure pour l’ordre de convergence : étant une méthode explicite et \(s<5\), l’ordre vérifie \(\omega\leq s\) donc \(\omega\leq 3\).
On peut calculer explicitement l’ordre de convergence avec sympy et on trouve qu’il n’est que d’ordre 2 :
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 )))
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**"))
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 = 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 = r"\sum_{j=1}^s b_j =" + f"{sum(b).simplify()}"
texte += r"\text{ doit être égale à }1"
display(Math(texte))
Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
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 ))
Eq_1.append(sympy.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=' +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 = r'\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 = r'\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]
s = 3
b = [1/sympy.S(6), 2/sympy.S(3), 1/sympy.S(6)]
c = [ 0, 1/sympy.S(2), 1]
A = [[0, 0, 0], [1/sympy.S(2), 0, 0], [0, sympy.S(1), 0]]
Eqs = ordre_RK(s, A=A, b=b, c=c)Matrice de Butcher
\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\\frac{1}{2} & \frac{1}{2} & 0 & 0\\1 & 0 & 1 & 0\\ & \frac{1}{6} & \frac{2}{3} & \frac{1}{6}\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}{2}\text{ doit être égale à }\frac{1}{2}\)
\(\displaystyle \sum_{j=1}^s a_{2j}=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{1}{3}\text{ doit être égale à }\frac{1}{3}\)
\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{1}{12}\text{ doit être égale à }\frac{1}{6}\)
# Test de l'ordre
# =============================================================================
import sympy as sp
t = sp.symbols('t')
y = sp.Function('y')(t)
t0, y0 = 1, 1
tfin = 2
phi_sp = lambda t,y : -y/t-y**2*sp.log(t)
sol_sp = sp.dsolve(y.diff(t)-phi_sp(t,y), y, ics={y.subs(t,t0):y0})
display(sol_sp)
# On convertit les expressions symboliques en fonctions numpy
sol_np = sp.lambdify(t, sol_sp.rhs, 'numpy')
phi = sp.lambdify([t,y], phi_sp(t,y), 'numpy')
# On résout le problème avec différentes valeurs de N
H, err = [], []
for N in range(10,100,10):
tt = np.linspace(t0,tfin,N)
uu = RK(phi,tt,y0)
yy = sol_np(tt)
err.append(np.linalg.norm(uu-yy,np.inf))
H.append(tt[1]-tt[0])
# On calcule la pente de la droite de meilleur ajustement en échelle log-log
pente = np.polyfit(np.log(H),np.log(err),1)[0]
plt.figure()
plt.loglog(H,err,'o-',label=f'pente = {pente:.2f}')
plt.xlabel('h')
plt.ylabel('erreur')
plt.legend() ;\(\displaystyle y{\left(t \right)} = \frac{2}{t \left(\log{\left(t \right)}^{2} + 2\right)}\)
Étudions la A-stabilité de la méthode en appliquant le schéma au problème ayant \(\varphi(t,y)=-\beta y\) avec \(\beta>0\). On a alors
\( \begin{cases} k_1 = -\beta u_n \\ k_2 = -\beta( u_n+\frac{h}{2}k_1) = -\beta u_n+\frac{\beta^2 h}{2}u_n\\ k_3 = -\beta(u_n+h k_2) = -\beta u_n + \beta^2 h u_n - \frac{\beta^3 h^2}{2}u_n\\ u_{n+1} = u_n + \frac{h}{6}(k_1+4k_2+k_3) = \left( 1 -(\beta h) +\frac{(\beta h)^2}{2} - \frac{(\beta h)^3}{12} \right) u_n \end{cases} \)
Notons \(x=\beta h>0\) et \(q(x)\) le coefficient d’amplification. On a alors
\( q(x) = \frac{u_{n+1}}{u_n} = 1 - x + \frac{x^2}{2} - \frac{x^3}{12} \)
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 \(q\) est une fonction monotone décroissante ayant un point d’inflexion en \(x=2\) et telle que \(q(0)=1\) et \(q(x)\xrightarrow[x\to+\infty]{}-\infty\). On cherche alors les solutions de \(q(x)+1=0\).
Il n’est pas possible de résoudre analytiquement cette équation, on utilisera donc une méthode numérique pour déterminer la valeur de \(x\) pour laquelle le schéma est A-stable.
De même, le schéma génère une solution monotone ssi \(q(x)>0\) . On cherche alors numériquement les solutions de \(q(x)=0\).
u_n = sympy.symbols(r'u_n')
u_np1 = sympy.symbols(r'u_{n+1}', escape=True)
dt = sympy.symbols(r'h', positive=True)
x = sympy.symbols(r'x', positive=True)
k_0 = sympy.symbols(r'k_0')
k_1 = sympy.symbols(r'k_1')
k_2 = sympy.symbols(r'k_2')
kk = [k_0, k_1, k_2]
beta = sympy.symbols(r'\beta', positive=True)
phi = lambda y: -beta*y
Eqs = [ sympy.Eq( kk[i] , phi(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)
q = (schema.rhs/u_n).simplify().subs(beta*dt,x)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('q(x)'),q.simplify()))))
# Tracer q(x) avec matplotlib
x_vals = np.linspace(0, 5, 400)
q_func = sympy.lambdify(x, q, 'numpy')
y_vals = q_func(x_vals)
plt.plot(x_vals, y_vals, label='q(x)')
plt.axhline(y=1, color='r', linestyle='--', label='y=1')
plt.axhline(y=-1, color='r', linestyle='--', label='y=-1')
plt.ylim(-1.5, 1.5)
plt.xlabel('$x$')
plt.ylabel('$q(x)$')
# Conditions de stabilité et de monotonie
from scipy.optimize import fsolve
x_A = fsolve(lambda x: q_func(x)+1, 4)[0] # q(x) = -1
x_M = fsolve(lambda x: q_func(x) , 2)[0] # q(x) = 0
display(Markdown(f"La condition de A-stabilité est satisfaite pour $\\beta h < {x_A:.6f}$"))
display(Markdown(f"La condition de monotonie est satisfaite pour $\\beta h < {x_M:.6f}$"))
plt.axvline(x=x_A, color='g', linestyle=':', label=f'$x_A = {x_A:.6f}$')
plt.axvline(x=x_M, color='b', linestyle=':', label=f'$x_M = {x_M:.6f}$');
\(\displaystyle k_{0} = - \beta u_{n}\)
\(\displaystyle k_{1} = - \beta \left(\frac{h k_{0}}{2} + u_{n}\right)\)
\(\displaystyle k_{2} = - \beta \left(h k_{1} + u_{n}\right)\)
\(\displaystyle u_{n+1} = h \left(- \frac{\beta^{3} h^{2} u_{n}}{12} + \frac{\beta^{2} h u_{n}}{2} - \beta u_{n}\right) + u_{n}\)
\(\displaystyle q(x) = - \frac{x \left(x^{2} - 6 x + 12\right)}{12} + 1\)
/usr/lib/python3/dist-packages/scipy/optimize/_minpack_py.py:177: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
warnings.warn(msg, RuntimeWarning)La condition de A-stabilité est satisfaite pour \(\beta h < 4.519842\)
La condition de monotonie est satisfaite pour \(\beta h < 2.001600\)
\( \begin{cases} y'(t) = t^2y(t), \\ y(0)=1, \end{cases} \quad t\in[0,1]. \)
Correction
Un schéma multi-pas linéaire explicite à \(q=3\) trois pas s’écrit
\( u_{n+1} = a_0 u_n + a_1 u_{n-1} + a_2 u_{n-2} + h \Big( b_0 \varphi_n + b_1 \varphi_{n-1} + b_2 \varphi_{n-2} \Big). \)
:one: Calcule des coefficients \(a_0\), \(a_1\), \(a_2\). Le premier polynomial caractéristique est
\( \varrho(r) = r^3 - a_0 r^2 - a_1 r - a_2 \)
qui a 3 racines. On sait que \(r=0\) est une racine double et que le schéma est convergent, donc \(r=1\) est une racine simple de \(\varrho\). Le polynôme caractéristique est alors
\( \varrho(r) = (r-1) r^2 = r^3 - r^2. \)
Par identification des coefficients, on trouve \(a_0=1\) et \(a_1=a_2=0\).
:two: Calcul des coefficients \(b_0\), \(b_1\), \(b_2\).
Les valeurs de \(b_0\) et \(b_2\) sont donnés dans l’énoncé :
\( u_{n+1} = u_n + \frac{h}{12}\Big( 5 \varphi_n + 12 b_1 \varphi_{n-1} + 23 \varphi_{n-2} \Big). \)
Il reste à déterminer \(b_1\). Le schéma étant convergent, donc consistant et zéro-stable. La condition de stabilité est donnée par les deux relations \(\xi(0)=1\) et \(\xi(1)=1\). La condition \(\xi(0)=1\) est déjà satisfaite, car cela équivaut à dire que \(r=1\) est une racine du polynôme caractéristique. La condition \(\xi(1)=1\) donne
\( -a_1+2a_2+b_0+b_1+b_2 = 1 \)
ce qui permet de calculer \(b_1\) :
\( b_1 = 1 - b_0 - b_2 = 1 - \frac{5}{12} - \frac{23}{12} = 1-\frac{28}{12} = -\frac{16}{12} =-\frac{4}{3}. \)
Le schéma cherché est donc
\( \boxed{ u_{n+1} = u_n + h \Big( \frac{5}{12} \varphi_n - \frac{4}{3} \varphi_{n-1} + \frac{23}{12} \varphi_{n-2} \Big). } \)
:three: Calcul de l’ordre de convergence. Il est d’ordre 1 car :
\( \xi(2) = a_{1} + 4 a_{2} - 2 \left(b_{1} + 2 b_{2} - b_{-1}\right) = 0 + 0 - 2 \left(-\frac{4}{3} + 2 \frac{23}{12} \right) = - 5 \neq 1 \)
:four: Implémentation du schéma et vérification expérimentale de l’ordre de convergence.
import numpy as np
import matplotlib.pyplot as plt
def MP(phi,tt,sol_exacte):
h = tt[1] - tt[0]
q = 3
uu = np.zeros_like(tt)
uu[:q] = sol_exacte(tt[:q])
for n in range(len(tt)-1):
k0 = phi(tt[n], uu[n])
k1 = phi(tt[n-1], uu[n-1])
k2 = phi(tt[n-2], uu[n-2])
uu[n+1] = uu[n] + h*(5*k0 - 16*k1 +23*k2)/12
return uu
# PB de Cauchy
t0, y0 = 0, 1
phi = lambda t,y: t**2*y
import sympy as sp
t = sp.symbols('t')
y = sp.Function('y')
eq = sp.Eq(y(t).diff(t), t**2*y(t))
sol = sp.dsolve(eq, y(t), ics={y(0): 1})
display(sol)
sol_exacte = sp.lambdify(t, sol.rhs, 'numpy')
H, err = [], []
for N in range(100,1000,100):
tt = np.linspace(t0, 1, N)
uu = MP(phi,tt,sol_exacte)
yy = sol_exacte(tt)
H.append(tt[1]-tt[0])
err.append(np.linalg.norm(uu-yy, np.inf))
ordre = np.polyfit(np.log(H), np.log(err), 1)[0]
plt.loglog(H,err, label=f"ordre = {ordre:.2f}")
plt.xlabel('h')
plt.ylabel('Erreur')
plt.legend();\(\displaystyle y{\left(t \right)} = e^{\frac{t^{3}}{3}}\)
