2026 Contrôle Terminale session 1 – M64 – Data Science & Scientific Computing 6 : Schémas numériques pour les problèmes de Cauchy

2026 Contrôle Terminale session 1

Auteur·rice

Gloria Faccanoni

Date de publication

25 avril 2026

⏱ 2h, déposer le notebook dans Moodle (si 1/3 temps supplémentaire, 2h40)

📚 Autorisés : notebooks de cours, notes personnelles.

⚠️ le contrôle est composé de 3 exercices indépendants.

Compléter ici :

  • NOM :
  • Prénom :
  • N° étudiant :
Code 
import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import fsolve

import sympy as sp

from IPython.display import display, Markdown

👣 Schémas multi-pas linéaires

Soit un schéma multi-pas linéaire convergent, défini par ses deux polynômes caractéristiques :

  • Premier polynôme caractéristique : \(\varrho(r)\)
  • Second polynôme caractéristique : \(\sigma(r)\)

Expliquer pourquoi, dans le cas d’un schéma multi-pas linéaire convergent, il est impossible d’avoir \(\varrho(r) = \sigma(r)\). Justifier votre réponse.

Réponse :

Un schéma convergent est consistant et zéro-stable.

  • Un schéma est consistant ssi \(\xi(0)=\xi(1)=1\).
  • Un schéma est zéro-stable ssi toutes les racines de son polynôme caractéristique sont dans le cercle unité et, si une racine est sur le cercle unité, elle est simple.

On sait que

\( \begin{cases} \xi(0)=1\\ \xi(1)=1 \end{cases} \iff \begin{cases} \varrho(1)=0\\ \varrho'(1)=\sigma(1) \end{cases} \)

Si \(\varrho=\sigma\), alors \(\varrho(1)=\sigma(1)\) et on a les trois conditions suivantes :

\( \begin{cases} \varrho(1)=0\\ \varrho'(1)=\sigma(1)\\ \varrho(1)=\sigma(1) \end{cases} \leadsto \varrho'(1)=\sigma(1)=\varrho(1)=0 \)

Autrement dit, 1 est une racine double du polynôme caractéristique, donc le schéma n’est pas zéro-stable.

👣 Schéma multi-pas linéaire implicite

  1. Déterminer les coefficients d’un schéma multi-pas linéaire implicite à deux pas sous les conditions suivantes :

    • \(r = \dfrac{1}{6}\) est une racine du premier polynôme caractéristique ;
    • l’ordre de convergence du schéma est supérieur ou égal à 3.

  2. On considère le problème de Cauchy \(\begin{cases} y'(t) = \dfrac{y^2(t)}{t}, & t \in [1, t_{\text{max}}], \\ y(1) = 4. & \end{cases}\) Déterminer la valeur maximale \(t_{\max}\) pour laquelle la solution est définie.
    (Indication : déterminer la solution exacte.)

  3. Implémenter le schéma déterminé à la question 1 pour résoudre le problème de Cauchy. Vérifier expérimentalement l’ordre de convergence du schéma sur l’intervalle \([1, \frac{1 + t_{\text{max}}}{2}]\).

Correction

Un schéma multi-pas linéaire implicite à \(q=2\) pas s’écrit comme

\( u_{n+1} = a_0 u_n + a_1 u_{n-1} + h \Big( b_0 \varphi_n + b_1 \varphi_{n-1} \Big) + h b_{-1}\varphi_{n+1}. \)

Le premier polynôme caractéristique a ordre 2 et est donné par

\( \varrho(r) = r^2 - a_0 r - a_1. \)

Le schéma est convergent, ainsi \(r=1\) est une racine. On sait aussi que \(r=\frac{1}{6}\) est l’autre racine, donc

\( \varrho(r) = (r-1)\left(r-\frac{1}{6}\right) = r^2 - \frac{7}{6} r + \frac{1}{6}. \)

Par identification des coefficients, on trouve \(a_0=\frac{7}{6}\) et \(a_1=-\frac{1}{6}\).

Il ne reste qu’à trouver les coefficients \(b_{-1}\), \(b_0\) et \(b_1\). Le schéma est convergent à l’ordre au moins 3 ssi \(\xi(0)=\xi(1)=\xi(2)=\xi(3)=1\). La condition \(\xi(0)=1\) est déjà satisfaite car \(r=1\) est une racine du polynôme caractéristique. On a alors les trois équations

\( \begin{cases} -a_1+b_0+b_1 + b_{-1} = 1\\ a_{1} - 2 b_{1} + 2 b_{-1} = 1\\ -a_1 + 3 b_{1} + 3 b_{-1} = 1\\ \end{cases} \quad \text{soit} \quad \begin{cases} \frac{1}{6}+b_0+b_1 + b_{-1} = 1\\ -\frac{1}{6} - 2 b_{1} + 2 b_{-1} = 1\\ \frac{1}{6} + 3 b_{1} + 3 b_{-1} = 1 \end{cases} \quad \text{ainsi} \quad \begin{cases} b_0 = \frac{40}{72}\\ b_1 = -\frac{11}{72}\\ b_{-1} = \frac{31}{72} \end{cases} \)

Conclusion : un seul schéma satisfait à toutes les conditions, il s’agit du schéma

\( u_{n+1} = \frac{7}{6} u_n - \frac{1}{6} u_{n-1} + \frac{h}{72} \left( 40 \varphi_n - 11 \varphi_{n-1} + 31 \varphi_{n+1} \right). \)

Nous pouvons vérifier nos calculs avec sympy :

Code 
p = 1     # p+1 coeffs du schema d'indices 0 à p

# ==========================================================
# Les équations
# ==========================================================
def xi(i, q):
    aa = sp.symbols(f'a_0:{q}')
    bb = sp.symbols(f'b_0:{q}')
    bm1 = sp.Symbol('b_{-1}')
    p = len(aa)
    sa = sum((-j)**i * aa[j] for j in range(p))
    sb = bm1 + sum((-j)**(i-1) * bb[j] for j in range(1, p))
    # Cas particulier : si j=0 et i=0, il calcule (-j)**i =1 et on a bien a_0 dans la somme
    # mais si j=0 et i=1 il ne calcule pas b_0 et il faut l'ajouter à la main
    if i == 1:
        sb += bb[0]
    return (sa).factor() + (i * sb).factor()

q = p + 1
display(Markdown(f"Pour une méthode à $q={p + 1}$ pas, "))
systeme = [sp.Eq(sp.Symbol(f"\\xi({i})") , sp.Eq(xi(i, q),1)) for i in range(4)]
display(Markdown(r"\(\begin{cases}"+",\\\\ ".join([sp.latex(eq) for eq in systeme])+r"\end{cases}\)"))


# ================================
# Calcul des coefficients 
# ================================

a0 = sp.symbols('a_0')
a1 = sp.symbols('a_1')
b0 = sp.symbols('b_0')
b1 = sp.symbols('b_1')
bm1 = sp.symbols('b_{-1}')
r = sp.symbols('r')

# Calculons a_0 et a_1
Eq0 = sp.Eq( (r-1)*(r-1/sp.S(6)) , r**2 - a0*r-a1 )
display(Eq0)
display(Markdown('**On trouve ainsi $a_0$ et $a_1$**'))
aa = sp.solve(Eq0, (a0,a1))
for k,v in aa.items():
    display(sp.Eq(k,v))
a0 = aa[a0]
a1 = aa[a1]

# Calculons b_0, b_1 et b_{-1}
display(Markdown(r'**On utilise $\xi(1)=\xi(2)=\xi(3)=1$ pour déterminer $b_0$, $b_1$ et $b_{-1}$**'))
sol = sp.solve( [ -a1+b0+b1+bm1 - 1,
                  a1-2*(b1-bm1) - 1,
                  -a1+3*(b1+bm1) -1
                ],
                (b0, b1, bm1)
              )
for k,v in sol.items():
    display(sp.Eq(k,v))

b0 = sol[b0]
b1 = sol[b1]
bm1 = sol[bm1]

Pour une méthode à \(q=2\) pas,

\(\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\end{cases}\)

\(\displaystyle \left(r - 1\right) \left(r - \frac{1}{6}\right) = - a_{0} r - a_{1} + r^{2}\)

On trouve ainsi \(a_0\) et \(a_1\)

\(\displaystyle a_{0} = \frac{7}{6}\)

\(\displaystyle a_{1} = - \frac{1}{6}\)

On utilise \(\xi(1)=\xi(2)=\xi(3)=1\) pour déterminer \(b_0\), \(b_1\) et \(b_{-1}\)

\(\displaystyle b_{0} = \frac{5}{9}\)

\(\displaystyle b_{1} = - \frac{11}{72}\)

\(\displaystyle b_{-1} = \frac{31}{72}\)

Calculons la solution exacte du problème de Cauchy. L’EDO est à variables séparables : \(\frac{dy}{y^2} = \frac{dt}{t}\). En intégrant de chaque côté on trouve \(-\frac{1}{y} = \ln|t| + C\). Avec \(y(1) = 4\) on trouve la constante d’intégration : \(-\frac{1}{4} = \ln(1) + C \leadsto C = -\frac{1}{4}\). Donc \(y(t) = \frac{1}{\frac{1}{4} - \ln(t)}\).

La solution exacte du problème de Cauchy est donnée par \(y(t) = \frac{4}{1-4\log(t)}\) qui n’est pas définie en \(t=\sqrt[4]{e}\). Ainsi, \(t_{\text{max}} = \sqrt[4]{e}\).

Code 
# ==================================================
# TEST
# ==================================================

phi = lambda t,y : y**2/t
t0, y0 = 1, 4


# Solution exacte
# ---------------
t = sp.symbols('t')
y = sp.Function('y')(t)
eq = sp.Eq(y.diff(t), phi(t, y))
display(eq)
ics = {y.subs(t, t0): y0 }
display(sp.Eq(y.subs(t, t0), y0))
sol = sp.dsolve(eq, y, ics=ics)
display(Markdown(f"La solution exacte est"))
display(sol)


display(Markdown(fr"On en déduit que $t < \sqrt[4]{{e}} \approx {np.exp(0.25)}$"))

# Transfomation en fonction numpy
# -------------------------------

sol = sp.lambdify(t, sol.rhs, 'numpy')
tfin = (t0+(np.exp(0.25)))/2

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = \frac{y^{2}{\left(t \right)}}{t}\)

\(\displaystyle y{\left(1 \right)} = 4\)

La solution exacte est

\(\displaystyle y{\left(t \right)} = - \frac{1}{\log{\left(t \right)} - \frac{1}{4}}\)

On en déduit que \(t < \sqrt[4]{e} \approx 1.2840254166877414\)

Code 
a0, a1 = float(a0), float(a1)
b0, b1, bm1 = float(b0), float(b1), float(bm1)


def MPL(phi, tt, sol_exacte) :
    uu = np.zeros_like(tt)
    h = tt[1]-tt[0]
    q = 2
    uu[:q] = sol_exacte(tt[:q])
    for n in range(q-1, len(tt)-1) :
        k0 = phi(tt[n], uu[n])
        k1 = phi(tt[n-1], uu[n-1])
        f = lambda x : -x + ( a0*uu[n] + a1*uu[n-1] ) + h*( b0*k0 + b1*k1 + bm1*phi(tt[n+1],x) ) 
        uu[n+1] = fsolve(f, uu[n])[0]
    return uu



# Solution numérique : vérification de l'implémentation
# -----------------------------------------------------
# tt = np.linspace(t0, tfin, 100)
# uu = MPL(phi, tt, sol)
# plt.plot(tt, uu, label='MPL')
# plt.plot(tt, sol(tt), label='Solution exacte')
# plt.legend();
# plt.figure()

# Calcul de l'ordre de convergence
# --------------------------------

H, err = [], []
for N in range(50,500,50) :
    tt = np.linspace(t0, tfin, N)
    uu = MPL(phi, tt, sol)
    err.append(np.linalg.norm(uu-sol(tt), np.inf))
    H.append(tt[1]-tt[0])

pente = np.polyfit(np.log(H), np.log(err), 1)[0]
plt.loglog(H, err, 'o-')
plt.title(f"Ordre de convergence = {pente = :g}")
plt.xlabel('h')
plt.ylabel('Erreur');

🎨 Étude de l’A-stabilité d’un schéma RK explicite à 3 étages

On considère le schéma de Runge-Kutta explicite à 3 étages défini par le tableau de Butcher suivant :

\( \def\arraystretch{2} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 \\ 2/3 & 0 & 2/3 & 0 \\ \hline & 1/4 & 0 & 3/4 \end{array} \)

On considère le problème de Cauchy \(y'(t) = -\beta y(t)\) et \(y(0) = 1\) avec \(\beta>0\).

  1. Déterminer les valeurs de \(h\) pour lesquelles le schéma est A-stable. On notera \(h_A(\beta)\) la plus grande valeur \(h\) assurant la A-stabilité.
  2. Déterminer les valeurs de \(h\) garantissant que le schéma produit une solution monotone. On notera \(h_M(\beta)\) la plus grande valeur de \(h\) assurant la monotonie.
  3. Implémenter le schéma et vérifier expérimentalement son A-stabilité en prenant \(h=h_M(\beta)/2\).

Le schéma est donné par

\( \begin{aligned} k_1 &= \varphi\Big(t_n, u_n\Big) \\ k_2 &= \varphi\Big(t_n + \frac{1}{3}h, u_n + \frac{1}{3}hk_1\Big) \\ k_3 &= \varphi\Big(t_n + \frac{2}{3}h, u_n + \frac{2}{3}hk_2 \Big) \\ u_{n+1} &= u_n + \frac{1}{4}h\Big(k_1 + 3k_3\Big) \end{aligned} \)

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}\). Avec \(\varphi(t, y) = -\beta y\), on obtient alors la relation de récurrence suivante :

\( \begin{aligned} k_1 &= -\beta u_n \\ k_2 &= -\beta\Big(u_n + \frac{1}{3}hk_1\Big) = -\beta u_n - \frac{\beta h}{3}k_1 = k_1 - \frac{\beta h}{3}k_1 = \Big(1 - \frac{1}{3}\beta h\Big)k_1 \\ k_3 &= -\beta\Big(u_n + \frac{2}{3}hk_2\Big) = -\beta u_n - \frac{2}{3}\beta hk_2 = k_1 - \frac{2}{3}\beta h\Big(1 - \frac{1}{3}\beta h\Big)k_1 = \Big(\frac{2}{9}(\beta h)^2 - \frac{2}{3}(\beta h) + 1\Big) k_1 \\ u_{n+1} &= u_n + \frac{1}{4}h\Big(k_1 + 3k_3\Big) = \left(-\frac{(\beta h)^3}{6} + \frac{(\beta h)^2}{2} -(\beta h)+1\right)u_n \end{aligned} \)

On notera \(x = \beta h>0\) et \(q(x)\) le facteur d’amplification du schéma :

\( q(x) = \frac{u_{n+1}}{u_n} = -\frac{x^3}{6} + \frac{x^2}{2} -x+1 \)

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 avec \(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 fsolve pour déterminer une valeur approchée 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\).

Code 
q_func = lambda x: -x**3/6+x**2/2-x+1

x = np.linspace(0,4,101)
plt.plot(x, q_func(x))
plt.hlines(1,0,4, color='green')
plt.hlines(-1,0,4, color='green')
plt.hlines(0,0,4, color='red')
plt.ylim(-1.5,1.5)
plt.xlim(0,4)
plt.grid()
plt.xlabel('x')
plt.ylabel('q(x)')
plt.title('q(x) en fonction de x')

h_A = fsolve(lambda x: q_func(x)+1, 2)[0]
h_M = fsolve(lambda x: q_func(x), 1)[0]

plt.scatter(h_A, -1, color='green')
plt.annotate(f'$h_A = {h_A:.6f}$', xy=(h_A, -1), xytext=(h_A+0.2, -0.5), arrowprops=dict(arrowstyle='->', lw=1.5), fontsize=12, color='green')
plt.scatter(h_M, 0, color='red')
plt.annotate(f'$h_M = {h_M:.6f}$', xy=(h_M, 0), xytext=(h_M+0.2, 0.5), arrowprops=dict(arrowstyle='->', lw=1.5), fontsize=12, color='red')

display(Markdown(f"La condition de A-stabilité est satisfaite pour $\\beta h < {h_A:.6f}$"))
display(Markdown(f"La condition de monotonie est satisfaite pour $\\beta h < {h_M:.6f}$"))

La condition de A-stabilité est satisfaite pour \(\beta h < 2.512745\)

La condition de monotonie est satisfaite pour \(\beta h < 1.596072\)

On peut vérifier nos calculs avec sympy :

Code 
import sympy as sp

b = [1/sp.S(4), 0, 3/sp.S(4)]
c = [0, 1/sp.S(3), 2/sp.S(3)]
A = [[0, 0, 0], [1/sp.S(3), 0, 0], [0, 2/sp.S(3), 0]]

s = len(c)

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_1 = sp.symbols(r'k_1')
k_2 = sp.symbols(r'k_2')
k_3 = sp.symbols(r'k_3')
kk = [k_1, k_2, k_3]

beta = sp.symbols(r'\beta', positive=True)
phi = lambda y: -beta*y

display(Markdown("**Les équations**"))
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) ]
for eq in Eqs:
    display(eq)

display(Markdown("**Ce qui donne**"))
sol = sp.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
for i,k in enumerate(KK):
    display(sp.Eq(kk[i]/u_n, k.factor()/u_n))

display(Markdown("**Enfin**"))
schema = sp.Eq( u_np1/u_n , (u_n + dt*sum([KK[j]*b[j] for j in range(s)])).factor()/u_n )
display(schema)


display(Markdown("**Le facteur d'amplification**"))
q = (schema.rhs).simplify().subs(beta*dt,x)
q_func = sp.lambdify(x, q, modules=['numpy'])
display(sp.Eq(sp.Symbol('q(x)'), q.factor()))

# sp.plot(q, (x,0,10), ylim=(-1.5,1.5),
#             ylabel='$q(x)$', xlabel='$x$')

# # sp.solve(abs(q)<1, x)
display(Markdown(f"La **condition de A-stabilité** est satisfaite pour $\\beta h < {sp.nsolve(abs(q)-1, x, 2):.6f}$"))
display(Markdown(f"La **condition de monotonie** est satisfaite pour $\\beta h < {sp.nsolve(abs(q), x, 2):.6f}$"))

Les équations

\(\displaystyle k_{1} = - \beta u_{n}\)

\(\displaystyle k_{2} = - \beta \left(\frac{h k_{1}}{3} + u_{n}\right)\)

\(\displaystyle k_{3} = - \beta \left(\frac{2 h k_{2}}{3} + u_{n}\right)\)

Ce qui donne

\(\displaystyle \frac{k_{1}}{u_{n}} = - \beta\)

\(\displaystyle \frac{k_{2}}{u_{n}} = \frac{\beta \left(\beta h - 3\right)}{3}\)

\(\displaystyle \frac{k_{3}}{u_{n}} = - \frac{\beta \left(2 \beta^{2} h^{2} - 6 \beta h + 9\right)}{9}\)

Enfin

\(\displaystyle \frac{u_{n+1}}{u_{n}} = - \frac{\beta^{3} h^{3}}{6} + \frac{\beta^{2} h^{2}}{2} - \beta h + 1\)

Le facteur d’amplification

\(\displaystyle q(x) = - \frac{x^{3} - 3 x^{2} + 6 x - 6}{6}\)

La condition de A-stabilité est satisfaite pour \(\beta h < 2.512745\)

La condition de monotonie est satisfaite pour \(\beta h < 1.596072\)

Vérifions empiriquement la stabilité du schéma pour \(\beta = 1\).

Code 
def RK(phi, tt, yy):
    h = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = yy[0]
    for n in range(len(tt)-1):
        k1 = phi(tt[n], uu[n])
        k2 = phi(tt[n] + h/3, uu[n] + h/3*k1)
        k3 = phi(tt[n] + h*2/3, uu[n] + 2/3*h*k2)
        uu[n+1] = uu[n] + h*(k1 + 3*k3)/4
    return uu


sol_exacte = lambda t: np.exp(-t)
phi = lambda t, y: -y

tt_M = np.arange(0, 25, h_M/2)
yy_M = sol_exacte(tt_M)
uu_M = RK(phi, tt_M, yy_M)

plt.plot(tt_M, uu_M, '^:', label='RK', lw=2)
plt.plot(tt_M, yy_M, 'o--', label='Solution exacte', lw=3)
plt.legend()
plt.show()