⏱ 2h, déposer le notebook dans Moodle (si 1/3 temps supplémentaire, 2h40)
📚 Autorisés : notebooks de cours, notes personnelles.
⚠️ Le contrôle continu est composé de 3 exercices indépendants. Deux versions de chaque exercice sont proposées ou bien l’énoncé dépend de votre numéro étudiant.
Compléter ici :
On considère le système d’équations différentielles suivant
\( \begin{cases} x'(t) = 2x(t)+3y^2(t),\\ y'(t) = 3x^2(t)-2y(t),\\ x(0) = -\frac{1}{2},\\ y(0) = \frac{1}{2}. \end{cases} \quad t \in [0, 5]. \)
Vérifier analytiquement que \(I(t) = x^3(t)-2x(t)y(t)-y^3(t)\) est un invariant.
Calculer une solution approchée avec la méthode de Heun avec \(30\) sous-intervalles. Tracer ensuite
Pour estimer l’ordre de convergence de la méthode de Heun, on ne peut pas utiliser la solution exacte car elle n’est pas connue. On peut cependant utiliser l’invariant \(I(t)\) dont on connaît la valeur exacte. Vérifier sur l’invariant que l’ordre de convergence est au moins 2. Avez-vous une explication pour le résultat obtenu ?
Correction
import numpy as np
import matplotlib.pyplot as plt
# 2D system
# x' = 2x +3y^2
# y' = 3x^2 - 2y
# x(0) = -1/2
# y(0) = 1/2
# Invariant: x^3 - y^3 - 2xy
# ==============================
t0 = 0
tfin = 5
x0 = -1/2
y0 = 1/2
phi = lambda t, zz : np.array([ 2*zz[0]+3*zz[1]**2,
3*zz[0]**2 - 2*zz[1] ])
Invariant = lambda zz : zz[0]**3 - zz[1]**3 -2*zz[0]*zz[1]
I_init = Invariant([x0, y0])On vérifie que \(I\) est bien un invariant pour le système donné. On peut calculer \(I'\) à la main ou utiliser sympy.
Notons \(\mathcal{I}(x,y) = x^3 - y^3 - 2xy\) ainsi \(I(t)=\mathcal{I}(x(t),y(t))\).
\(\begin{aligned} I'(t) &= \frac{\partial \mathcal{I}}{\partial x} \cdot x'(t) + \frac{\partial \mathcal{I}}{\partial y} \cdot y'(t)\\ &= \Big(3x^2(t)-2y(t)\Big)x'(t) - y'(t) \Big(3y^2(t)+2x(t)\Big)\\ &= y'(t)x'(t)-y'(t)x'(t)\\ &= 0. \end{aligned}\)
On a donc bien \(I'(t) = 0\) ainsi \(I(t)\) est invariant et vaut \(I(0) = \frac{1}{4}\) pour tout \(t\).
import sympy as sp
t = sp.symbols('t')
x = sp.Function('x')(t)
y = sp.Function('y')(t)
eq1 = sp.Eq(sp.diff(x, t), phi(t, [x, y])[0])
eq2 = sp.Eq(sp.diff(y, t), phi(t, [x, y])[1])
I = Invariant([x, y])
display(sp.Eq(sp.Symbol('I'), I))
dI = sp.diff(I, t)
dI = dI.subs(sp.diff(x, t), eq1.rhs)
dI = dI.subs(sp.diff(y, t), eq2.rhs)
dI = dI.simplify()
display(sp.Eq(sp.Symbol("I'(t)"), dI))
display(sp.Eq(sp.Symbol("I(0)"), I_init))\(\displaystyle I = x^{3}{\left(t \right)} - 2 x{\left(t \right)} y{\left(t \right)} - y^{3}{\left(t \right)}\)
\(\displaystyle I'(t) = 0\)
\(\displaystyle I(0) = 0.25\)
Même si ce n’est pas explicitement demandé dans l’énoncé, c’est toujours une bonne idée d’afficher le portrait de phase pour ce genre de système. Ci-dessous on affiche, de plus, la courbe de niveau de l’invariant \(I(t)\) correspondant à \(I(t) = I(0)\).
# Plotting the vector field
# ==========================================
Y1,Y2 = np.meshgrid(np.linspace(-1.5,1.5,41),np.linspace(-1.5,1.5,41))
tt = np.linspace(t0,tfin,100)
V1,V2 = phi(tt,[Y1,Y2])
r1 = np.sqrt(1+V1**2)
r2 = np.sqrt(1+V2**2)
M = np.sqrt(V1**2 + V2**2)
plt.figure(figsize=(10,7))
plt.quiver(Y1, Y2, V1/r1, V2/r2, M, cmap=plt.cm.viridis, scale_units='xy',scale=10.)
plt.contour(Y1,Y2,Invariant([Y1,Y2]),levels=[I_init],colors='r',linewidths=2)
plt.plot(x0,y0,'bo')
plt.grid()
plt.xlabel(r'$x$')
plt.ylabel(r'$y$');
# Schéma de Heun
# =============================================================================
def Heun_system(phi,tt, yy0):
nb_eq = len(yy0) # nombre d'équations
N = len(tt)-1 # N+1 points, donc N intervalles
uu = np.zeros((N+1,nb_eq)) # colonnes = variables, lignes = temps
uu[0] = yy0
h = tt[1]-tt[0]
for n in range(N):
K1 = phi(tt[n],uu[n])
K2 = phi(tt[n+1],uu[n]+h*K1)
uu[n+1] = uu[n] + h/2 * ( K1 + K2 )
return uu# TEST 1
# =======
N = 30
t_values = np.linspace(t0, tfin, N+1)
yy0 = np.array([x0, y0])
uu = Heun_system(phi, t_values, yy0)
# from scipy.integrate import odeint
# uu = odeint(phi, yy0, t_values, tfirst=True)
II = Invariant(uu.T)
fig, axx = plt.subplots(1,3, figsize=(10,5))
axx[0].plot(t_values, uu[:,0], 'o-',label='x')
axx[0].plot(t_values, uu[:,1], 'd-',label='y')
axx[0].set_xlabel('t')
axx[0].legend()
axx[0].grid()
axx[1].plot(uu[:,0], uu[:,1], 'o-')
axx[1].set_xlabel('x')
axx[1].set_ylabel('y')
axx[1].grid()
axx[2].plot(t_values, II, 'o-')
axx[2].plot(t_values, I_init*np.ones(N+1), 'r--')
axx[2].set_xlabel('t')
axx[2].set_ylabel('Invariant')
axx[2].grid()
plt.tight_layout()
# TEST 2 : ordre de convergence
# ==============================
H, err = [], []
for k in range(7):
N += 2**k
t_values = np.linspace(t0, tfin, N+1)
uu = Heun_system(phi, t_values, yy0)
II = Invariant(uu.T)
err.append(np.linalg.norm(II-I_init))
H.append(t_values[1]-t_values[0])
if err[-1] <1.e-14 :
break
plt.loglog(H, err, 'o-', label='erreur')
a,b = np.polyfit(np.log(H), np.log(err), 1)
plt.loglog(H, np.exp(a*np.log(H)+b), label=f'pente {a:.2f}')
plt.legend()
plt.show()
L’ordre calculé sur l’invariant ne peut pas être inférieur à l’ordre théorique sur la solution. En revanche, il peut être meilleur grâce à la structure particulière de l’invariant et aux compensations possibles entre les erreurs sur chaque composante de la solution.
On considère le système d’équations différentielles suivant
\( \begin{cases} x'(t) = -(2x(t)+3y^2(t)),\\ y'(t) = 2y(t)-3x^2(t),\\ x(0) = -\frac{1}{2},\\ y(0) = \frac{1}{2}. \end{cases} \quad t \in [0, 5]. \)
Vérifier analytiquement que \(I(t) = y^3(t)-x^3(t)+2x(t)y(t)\) est un invariant.
Calculer une solution approchée avec la méthode d’Euler Modifiée avec \(30\) sous-intervalles. Tracer ensuite
Pour estimer l’ordre de convergence de la méthode d’Euler Modifiée, on ne peut pas utiliser la solution exacte car elle n’est pas connue. On peut cependant utiliser l’invariant \(I(t)\) dont on connaît la valeur exacte. Vérifier sur l’invariant que l’ordre de convergence est au moins 2. Avez-vous une explication pour le résultat obtenu ?
Correction
import numpy as np
import matplotlib.pyplot as plt
# 2D system
# x' = -2x -3y^2
# y' = -3x^2 + 2y
# x(0) = -1/2
# y(0) = 1/2
# Invariant: y^3 - x^3 + 2xy
# ==============================
t0 = 0
tfin = 5
x0 = -1/2
y0 = 1/2
phi = lambda t, zz : np.array([ -2*zz[0] -3*zz[1]**2,
-3*zz[0]**2 + 2*zz[1] ])
Invariant = lambda zz : zz[1]**3 - zz[0]**3 +2*zz[0]*zz[1]
I_init = Invariant([x0, y0])On vérifie que \(I\) est bien un invariant pour le système donné. On peut calculer \(I'\) à la main ou utiliser sympy.
Notons \(\mathcal{I}(x,y) = y^3 - x^3 + 2xy\) ainsi \(I(t)=\mathcal{I}(x(t),y(t))\).
\(\begin{aligned} I'(t) &= \frac{\partial \mathcal{I}}{\partial x} \cdot x'(t) + \frac{\partial \mathcal{I}}{\partial y} \cdot y'(t)\\ &= \Big(3y^2(t)+2x(t)\Big)y'(t) - x'(t) \Big(3x^2(t)-2y(t)\Big)\\ &= -x'(t)y'(t)+x'(t)y'(t)\\ &= 0. \end{aligned}\)
On a donc bien \(I'(t) = 0\) ainsi \(I(t)\) est invariant et vaut \(I(0) = -\frac{1}{4}\) pour tout \(t\).
import sympy as sp
t = sp.symbols('t')
x = sp.Function('x')(t)
y = sp.Function('y')(t)
eq1 = sp.Eq(sp.diff(x, t), phi(t, [x, y])[0])
eq2 = sp.Eq(sp.diff(y, t), phi(t, [x, y])[1])
I = Invariant([x, y])
display(sp.Eq(sp.Symbol('I'), I))
dI = sp.diff(I, t)
dI = dI.subs(sp.diff(x, t), eq1.rhs)
dI = dI.subs(sp.diff(y, t), eq2.rhs)
dI = dI.simplify()
display(sp.Eq(sp.Symbol("I'(t)"), dI))
display(sp.Eq(sp.Symbol("I(0)"), I_init))\(\displaystyle I = - x^{3}{\left(t \right)} + 2 x{\left(t \right)} y{\left(t \right)} + y^{3}{\left(t \right)}\)
\(\displaystyle I'(t) = 0\)
\(\displaystyle I(0) = -0.25\)
Même si ce n’est pas explicitement demandé dans l’énoncé, c’est toujours une bonne idée d’afficher le portrait de phase pour ce genre de système. Ci-dessous on affiche, de plus, la courbe de niveau de l’invariant \(I(t)\) correspondant à \(I(t) = I(0)\).
# Plotting the vector field
# ==========================================
Y1,Y2 = np.meshgrid(np.linspace(-1.5,1.5,41),np.linspace(-1.5,1.5,41))
tt = np.linspace(t0,tfin,100)
V1,V2 = phi(tt,[Y1,Y2])
r1 = np.sqrt(1+V1**2)
r2 = np.sqrt(1+V2**2)
M = np.sqrt(V1**2 + V2**2)
plt.figure(figsize=(10,7))
plt.quiver(Y1, Y2, V1/r1, V2/r2, M, cmap=plt.cm.viridis, scale_units='xy',scale=10.)
plt.contour(Y1,Y2,Invariant([Y1,Y2]),levels=[I_init],colors='r',linewidths=2)
plt.plot(x0,y0,'bo')
plt.grid()
plt.xlabel(r'$x$')
plt.ylabel(r'$y$');
# Schéma d'Euler Modifié
# =============================================================================
def EM_system(phi,tt, yy0):
nb_eq = len(yy0) # nombre d'équations
N = len(tt)-1 # N+1 points, donc N intervalles
uu = np.zeros((N+1,nb_eq)) # colonnes = variables, lignes = temps
uu[0] = yy0
h = tt[1]-tt[0]
for n in range(N):
K1 = phi(tt[n],uu[n])
uu[n+1] = uu[n] + h*phi(tt[n]+h/2,uu[n]+h/2*K1)
return uu# TEST 1
#=======
N = 30
t_values = np.linspace(t0, tfin, N+1)
yy0 = np.array([x0, y0])
uu = EM_system(phi, t_values, yy0)
# from scipy.integrate import odeint
# uu = odeint(phi, yy0, t_values, tfirst=True)
II = Invariant(uu.T)
fig, axx = plt.subplots(1,3, figsize=(10,5))
axx[0].plot(t_values, uu[:,0], 'o-',label='x')
axx[0].plot(t_values, uu[:,1], 'd-',label='y')
axx[0].set_xlabel('t')
axx[0].legend()
axx[0].grid()
axx[1].plot(uu[:,0], uu[:,1], 'o-')
axx[1].set_xlabel('x')
axx[1].set_ylabel('y')
axx[1].grid()
axx[2].plot(t_values, II, 'o-')
axx[2].plot(t_values, I_init*np.ones(N+1), 'r--')
axx[2].set_xlabel('t')
axx[2].set_ylabel('Invariant')
axx[2].grid()
plt.tight_layout()
# TEST 2 : ordre de convergence
# ==============================
H, err = [], []
for k in range(7):
N += 2**k
t_values = np.linspace(t0, tfin, N+1)
uu = EM_system(phi, t_values, yy0)
II = Invariant(uu.T)
err.append(np.linalg.norm(II-I_init))
H.append(t_values[1]-t_values[0])
if err[-1] <1.e-14 :
break
plt.loglog(H, err, 'o-', label='erreur')
a,b = np.polyfit(np.log(H), np.log(err), 1)
plt.loglog(H, np.exp(a*np.log(H)+b), label=f'pente {a:.2f}')
plt.legend()
plt.show()
L’ordre calculé sur l’invariant ne peut pas être inférieur à l’ordre théorique sur la solution. En revanche, il peut être meilleur grâce à la structure particulière de l’invariant et aux compensations possibles entre les erreurs sur chaque composante de la solution.
La méthode multi-pas linéaire suivante est-elle convergente ? Justifier.
\( u_{n+1} = -\frac{1}{4} u_{n} + \frac{1}{2} u_{n-1} + \frac{3}{4} u_{n-2} + \frac{h}{8} (9\varphi_{n} + 5\varphi_{n-1}) \)
Correction
Une méthode multi-pas linéaire est convergente ssi elle est consistante et zéro-stable.
Cette méthode est bien zéro-stable, car toutes les racines de son polynôme caractéristique sont dans le cercle unité et seule la racine \(1\) est sur la frontière du cercle unité.
Cependant, elle n’est pas consistante, car \(\xi(1)\neq1\).
Vérifions-le avec symmpy :
from sympy import S, symbols, Symbol, latex, factor, solve, Eq, poly, Rational, roots
from IPython.display import display, Math, Markdown, Latex
# ========================
q = 3
# ========================
r = Symbol('r')
# ========================
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()
# ========================
aa = [-1/S(4), 1/S(2), 3/S(4)]
bb = [9/S(8), -5/S(8), 0]
b_m1 = 0
display(Markdown(f"C'est une méthode à ${q = }$ pas"))
display(Markdown(f"Pour étudier la zéro-stabilité, on écrit le premier polynôme caractéristique"))
rho = poly(r**(q) - sum( aa[j] * r**(q-1-j) for j in range(q)), r)
texte = r"$\varrho(r)=" + latex(rho.as_expr(),order='old') + r"$"
display(Math(texte))
display(Markdown(f"Il a 3 racines :"))
racines = roots(rho,r)
for k,v in racines.items():
display(Markdown(f" - ${latex(k)}$ de module ${latex(abs(k))} \\approx {float(abs(k))}$"))
display(Markdown(f"Pour étudier la consistance, on calcule $\\xi(i)$"))
for i in range(2):
display(Markdown(f"$\\xi({i}) = {latex(xi(i, q))}$"))C’est une méthode à \(q = 3\) pas
Pour étudier la zéro-stabilité, on écrit le premier polynôme caractéristique
\(\displaystyle \varrho(r)=- \frac{3}{4} - \frac{r}{2} + \frac{r^{2}}{4} + r^{3}\)
Il a 3 racines :
Pour étudier la consistance, on calcule \(\xi(i)\)
\(\xi(0) = 1\)
\(\xi(1) = - \frac{3}{2}\)
La méthode multi-pas linéaire suivante est-elle convergente ? Justifier.
\( u_{n+1} = -\frac{1}{4} u_{n} + \frac{1}{2} u_{n-1} +\frac{3}{4} u_{n-2} + \frac{h}{8} (19\varphi_{n} + 5\varphi_{n-2}) \)
Correction
Une méthode multi-pas est convergente ssi elle est consistante et zéro-stable.
Elle est zéro-stable, car toutes les racines de son polynôme caractéristique sont dans le cercle unité et seul \(1\) est sur la frontière du cercle unité.
Elle est aussi consistante d’ordre \(3\), car \(\xi(i)=1\) pour \(i=1,2,3\). Étant une méthode explicite à 3 pas, c’est l’ordre maximal possible d’après la première barrière de Dalquist.
Vérifions ces calculs avec sympy :
from sympy import S, symbols, Symbol, latex, factor, solve, Eq, poly, Rational, roots
from IPython.display import display, Math, Markdown, Latex
# ========================
q = 3
# ========================
r = Symbol('r')
# ========================
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()
# ========================
aa = [-1/S(4), 1/S(2), 3/S(4)]
bb = [19/S(8), 0, 5/S(8)]
b_m1 = 0
display(Markdown(f"C'est une méthode à ${q = }$ pas"))
display(Markdown(f"Pour étudier la zéro-stabilité, on écrit le premier polynôme caractéristique"))
rho = poly(r**(q) - sum( aa[j] * r**(q-1-j) for j in range(q)), r)
texte = r"$\varrho(r)=" + latex(rho.as_expr(),order='old') + r"$"
display(Math(texte))
display(Markdown(f"dont les racines sont :"))
racines = roots(rho,r)
for k,v in racines.items():
display(Markdown(f" - ${latex(k)}$ de module ${latex(abs(k))} \\approx {float(abs(k))}$"))
display(Markdown(f"Pour étudier la consistance (et l'ordre de consistance), on calcule $\\xi(i)$"))
for i in range(4):
display(Markdown(f"$\\xi({i}) = {latex(xi(i, q))}$"))C’est une méthode à \(q = 3\) pas
Pour étudier la zéro-stabilité, on écrit le premier polynôme caractéristique
\(\displaystyle \varrho(r)=- \frac{3}{4} - \frac{r}{2} + \frac{r^{2}}{4} + r^{3}\)
dont les racines sont :
Pour étudier la consistance (et l’ordre de consistance), on calcule \(\xi(i)\)
\(\xi(0) = 1\)
\(\xi(1) = 1\)
\(\xi(2) = 1\)
\(\xi(3) = 1\)
On considère le schéma de Runge-Kutta explicite à 2 étages défini par le tableau de Butcher suivant :
\( \begin{array}{c|cc} 1 & 1 & 0 \\ 1/3 & 1/3 & 0 \\ \hline & 1/4 & 3/4 \end{array} \)
Après avoir calculé l’ordre de convergence théorique, implémenter le schéma et vérifier expérimentalement son ordre de convergence sur un exemple que vous choisirez.
Correction
On calcule l’ordre de convergence théorique du schéma de Runge-Kutta à 2 étages donné.
C’est un schéma semi-implicite. D’après la barrière de Butcher, l’ordre est au plus 4.
En vérifiant les conditions de convergence, on trouve que l’ordre est 2.
# =============================================================================
# Fonctions pour calculer les conditions de consistance pour un schéma RK.
#
# La fonction ordre_RK calcule les conditions pour les ordres 1, 2, 3 et 4 (même si théoriquement cela n'est pas possible)
# - affiche la matrice de Butcher
# - affiche les conditions de consistance obtenues pour chaque ordre par la fonction ordre_i
# - renvoie une liste d'équations pour une éventuelle résolution
#
# Chaque fonction ordre_i calcule les conditions pour être d'ordre i
# - peut afficher les conditions de consistance [actuellement en commentaire]
# - renvoie la liste d'équations
#
# =============================================================================
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, type="implicite", A=None, b=None, c=None):
j = sympy.symbols('j')
if A is None:
if type == "implicite":
A = sympy.MatrixSymbol('a', s, s)
elif type=="semi-implicite":
A = sympy.Matrix(s, s, lambda i,j: sympy.Symbol('a{}{}'.format(i, j)) if i >= j else 0)
elif type=="explicite":
A = sympy.Matrix(s, s, lambda i,j: sympy.Symbol('a{}{}'.format(i, j)) if i > j else 0)
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(*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(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(*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)
display(*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]
# =============================================================================
# Calcul des conditions de consistance pour un schéma RK
# On doit indiquer le nombre d'étages s et le type de schéma
# parmi implicite, semi-implicite ou explicite
# =============================================================================
s = 2
type = "semi-implicite" # implicite, semi-implicite, explicite
A = [[1, 0], [1/sympy.S(3), 0]]
c = [1, 1/sympy.S(3)]
b = [1/sympy.S(4),3/sympy.S(4)]
display(Markdown(f"### Conditions pour un schéma {type} avec s = {s} étages "))
Eqs = ordre_RK(s=s, type=type)
display(Markdown(f"### Pour notre exemple on trouve "))
Eqs = ordre_RK(s=s, type=type, A=A, b=b, c=c)Matrice de Butcher
\(\displaystyle \left[\begin{matrix}c_{0} & a_{00} & 0\\c_{1} & a_{10} & a_{11}\\ & b_{0} & b_{1}\end{matrix}\right]\)
On a 3 conditions pour avoir consistance (= pour être d’ordre 1)
\(\displaystyle b_{0} + b_{1} = 1\)
\(\displaystyle a_{00} = c_{0}\)
\(\displaystyle a_{10} + a_{11} = c_{1}\)
On doit ajouter 1 condition pour être d’ordre 2
\(\displaystyle b_{0} c_{0} + b_{1} c_{1} = \frac{1}{2}\)
On doit ajouter 2 conditions pour être d’ordre 3
\(\displaystyle b_{0} c_{0}^{2} + b_{1} c_{1}^{2} = \frac{1}{3}\)
\(\displaystyle a_{00} b_{0} c_{0} + a_{10} b_{1} c_{0} + a_{11} b_{1} c_{1} = \frac{1}{6}\)
On doit ajouter 4 conditions pour être d’ordre 4
\(\displaystyle b_{0} c_{0}^{3} + b_{1} c_{1}^{3} = \frac{1}{4}\)
\(\displaystyle a_{00} b_{0} c_{0}^{2} + a_{10} b_{1} c_{0} c_{1} + a_{11} b_{1} c_{1}^{2} = \frac{1}{8}\)
\(\displaystyle a_{00} b_{0} c_{0}^{2} + a_{10} b_{1} c_{0}^{2} + a_{11} b_{1} c_{1}^{2} = \frac{1}{12}\)
\(\displaystyle a_{00}^{2} b_{0} c_{0} + a_{00} a_{10} b_{1} c_{0} + a_{10} a_{11} b_{1} c_{0} + a_{11}^{2} b_{1} c_{1} = \frac{1}{24}\)
Matrice de Butcher
\(\displaystyle \left[\begin{matrix}1 & 1 & 0\\\frac{1}{3} & \frac{1}{3} & 0\\ & \frac{1}{4} & \frac{3}{4}\end{matrix}\right]\)
On a 3 conditions pour avoir consistance (= pour être d’ordre 1)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
On doit ajouter 1 condition pour être d’ordre 2
\(\displaystyle \text{True}\)
On doit ajouter 2 conditions pour être d’ordre 3
\(\displaystyle \text{True}\)
\(\displaystyle \text{False}\)
On doit ajouter 4 conditions pour être d’ordre 4
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
On implemente le schéma et on vérifie expérimentalement son ordre de convergence sur le problème de Cauchy \(y'(t) = t^2e^{-t}-(3t^2+1)y(t)\) avec \(y(0) = 1\).
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
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])
k1 = fsolve(lambda k: k - phi(tt[n+1], uu[n]+h*k), k_init)[0]
k2 = phi(tt[n]+h/3, uu[n]+h*k1/3)
uu[n+1] = uu[n] + h*(k1+3*k2)/4
return uu
t0 = 0
y0 = 1
phi = lambda t, y: -(3*t**2+1)*y+t**2*np.exp(-t)
sol_exacte = lambda t: (2+np.exp(t**3))/3*np.exp(-t**3-t)
H, err = [], []
for N in range(50,1000,50):
tt = np.linspace(0, 1, N)
yy = sol_exacte(tt)
uu = RK(phi, tt, y0)
H.append(tt[1]-tt[0])
err.append(np.linalg.norm(yy-uu)) # norme L2
pente = np.polyfit(np.log(H), np.log(err), 1)[0]
plt.loglog(H, err, 'o-', label=f'pente={pente:.2f}')
plt.legend();
On considère le schéma de Runge-Kutta explicite à 2 étages défini par le tableau de Butcher suivant :
\( \begin{array}{c|cc} 0 & 0 & 0 \\ 2/3 & 1/3 & 1/3 \\ \hline & 1/4 & 3/4 \end{array} \)
Après avoir calculé l’ordre de convergence théorique, implémenter le schéma et vérifier expérimentalement son ordre de convergence sur un exemple que vous choisirez.
Correction
On calcule l’ordre de convergence théorique du schéma de Runge-Kutta à 2 étages donné.
C’est un schéma semi-implicite. D’après la barrière de Butcher, l’ordre est au plus 4.
En vérifiant les conditions de convergence, on trouve que l’ordre est 3.
# =============================================================================
# Fonctions pour calculer les conditions de consistance pour un schéma RK.
#
# La fonction ordre_RK calcule les conditions pour les ordres 1, 2, 3 et 4 (même si théoriquement cela n'est pas possible)
# - affiche la matrice de Butcher
# - affiche les conditions de consistance obtenues pour chaque ordre par la fonction ordre_i
# - renvoie une liste d'équations pour une éventuelle résolution
#
# Chaque fonction ordre_i calcule les conditions pour être d'ordre i
# - peut afficher les conditions de consistance [actuellement en commentaire]
# - renvoie la liste d'équations
#
# =============================================================================
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, type="implicite", A=None, b=None, c=None):
j = sympy.symbols('j')
if A is None:
if type == "implicite":
A = sympy.MatrixSymbol('a', s, s)
elif type=="semi-implicite":
A = sympy.Matrix(s, s, lambda i,j: sympy.Symbol('a{}{}'.format(i, j)) if i >= j else 0)
elif type=="explicite":
A = sympy.Matrix(s, s, lambda i,j: sympy.Symbol('a{}{}'.format(i, j)) if i > j else 0)
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(*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(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(*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)
display(*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]
# =============================================================================
# Calcul des conditions de consistance pour un schéma RK
# On doit indiquer le nombre d'étages s et le type de schéma implicite, semi-implicite ou explicite
# =============================================================================
s = 2
type = "semi-implicite" # implicite, semi-implicite, explicite
A = [[0, 0], [1/sympy.S(3), 1/sympy.S(3)]]
b = [1/sympy.S(4), 3/sympy.S(4)]
c = [0, 2/sympy.S(3)]
display(Markdown(f"### Conditions pour un schéma {type} avec s = {s} étages "))
Eqs = ordre_RK(s=s, type=type, A=A, b=b, c=c)Matrice de Butcher
\(\displaystyle \left[\begin{matrix}0 & 0 & 0\\\frac{2}{3} & \frac{1}{3} & \frac{1}{3}\\ & \frac{1}{4} & \frac{3}{4}\end{matrix}\right]\)
On a 3 conditions pour avoir consistance (= pour être d’ordre 1)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
On doit ajouter 1 condition pour être d’ordre 2
\(\displaystyle \text{True}\)
On doit ajouter 2 conditions pour être d’ordre 3
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
On doit ajouter 4 conditions pour être d’ordre 4
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
On implemente le schéma et on vérifie expérimentalement son ordre de convergence sur le problème de Cauchy \(y'(t) = t^2e^{-t}-(3t^2+1)y(t)\) avec \(y(0) = 1\).
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
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):
k1 = phi(tt[n], uu[n])
k2 = fsolve(lambda k: k - phi(tt[n]+h*2/3, uu[n]+h*(k1+k)/3), phi(tt[n], uu[n]))[0]
uu[n+1] = uu[n] + h*(k1+3*k2)/4
return uu
t0 = 0
y0 = 1
phi = lambda t, y: -(3*t**2+1)*y+t**2*np.exp(-t)
sol_exacte = lambda t: (2+np.exp(t**3))/3*np.exp(-t**3-t)
H, err = [], []
for N in range(50,1000,50):
tt = np.linspace(0, 1, N)
yy = sol_exacte(tt)
uu = RK(phi, tt, y0)
H.append(tt[1]-tt[0])
err.append(np.linalg.norm(yy-uu)) # norme L2
pente = np.polyfit(np.log(H), np.log(err), 1)[0]
plt.loglog(H, err, 'o-', label=f'pente={pente:.2f}')
plt.legend();