Propriétés du problème de Cauchy – M64 – Data Science & Scientific Computing 6 : Schémas numériques pour les problèmes de Cauchy

Propriétés du problème de Cauchy

Auteur·rice

Gloria Faccanoni

Date de publication

26 mars 2026

Code 
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve
import sympy as sp
from IPython.display import Markdown, display, Math
plt.rcParams.update({'font.size': 12})

En général, il ne suffit pas qu’un schéma numérique soit convergent pour qu’il donne de bons résultats sur n’importe quelle équation différentielle. Dès qu’on tente d’approcher une solution, de nombreux facteurs peuvent en fait nous en éloigner : une étude mathématique sophistiquée doit alors être envisagée. En particulier, il faut que le problème de Cauchy soit

  • mathématiquement bien posé : sa solution existe et est unique,
  • numériquement bien posé : sa solution n’est pas trop perturbée par une erreur initiale ou des perturbations du second membre,
  • bien conditionné (= non raide) : les méthodes numériques usuelles peuvent calculer une solution approchée raisonnable en un nombre raisonnable d’itérations.

Ces définitions sont celles adoptés aux pages 235-236 dans le livre de Jean-Pierre DEMAILLY, Analyse Numérique et équations différentielles.

Problème mathématiquement bien posé : existence et unicité de la solution

Définition de problème mathématiquement bien posé

On dit qu’un problème de Cauchy est mathématiquement bien posé s’il existe une et une seule solution.

Si ce n’est pas le cas, on dit qu’il est mathématiquement mal posé.

Exemple de problème mathématiquement mal posé

Le problème de Cauchy

\( \begin{cases} y'(t) = \sqrt[3]{y(t)}, & \text{pour } t>0,\\ y(0) = 0. \end{cases} \)

ne satisfait pas les hypothèses du théorème de Cauchy-Lipschitz, car la fonction \(\varphi(t,y)=\sqrt[3]{y}\) n’est pas lipschitzienne par rapport à la deuxième variable en \(y=0\) (on a \(\partial_y \varphi(t,y)=\frac{1}{3\sqrt[3]{y^2}} \xrightarrow[y\to0]{}+\infty\)).
On ne peut donc pas garantir l’unicité de la solution et, en effet, il est simple de vérifier que, pour tout \(t\ge0\), les trois fonctions

  • \(y_1(t)=0\) pour tout \(t\)
  • \(y_{2,3}(t)=\pm\sqrt{\frac{8}{27}t^3}=\pm\frac{2\sqrt{2}}{3\sqrt{3}}\sqrt{t^3}=\pm\frac{2\sqrt{6}}{9}\sqrt{t^3}\)

sont solution du problème de Cauchy donné.

Code 
# Remarque : si on demande à sympy de calculer toutes les solutions, il ne trouve pas la solution constante = 0

t = sp.Symbol('t')
y = sp.Function('y')

edo = sp.Eq( sp.diff(y(t),t) , y(t)**(1/sp.S(3))  )
t0, y0 = 0, 0
display(Markdown(f"Equation différentielle : ${sp.latex(edo)}$"))
display(Markdown(f"Condition initiale : $y({t0}) = {y0}$"))

solgenLIST = sp.dsolve(edo,y(t), ics={y(0):0} )
display(Markdown(f"Solution(s) : ${sp.latex(solgenLIST)}$"))

Equation différentielle : \(\frac{d}{d t} y{\left(t \right)} = \sqrt[3]{y{\left(t \right)}}\)

Condition initiale : \(y(0) = 0\)

Solution(s) : \(\left[ y{\left(t \right)} = - \frac{2 \sqrt{6} t^{\frac{3}{2}}}{9}, \ y{\left(t \right)} = \frac{2 \sqrt{6} t^{\frac{3}{2}}}{9}\right]\)

Que se passe-t-il lorsqu’on utilise un schéma pour un problème mathématiquement MAL posé ?

Lorsqu’on utilise une méthode numérique, on ne peut pas savoir quelle solution ce schéma approche et différents schémas peuvent approcher différentes solutions.

Dans l’exemple précédent, si on utilise la méthode d’Euler explicite, on trouve que \(u_n=0\) pour tout \(n\), ce qui correspond à la solution \(y_1\).

Correction

Sur \(\mathbb{R}^+\) la fonction \(\varphi(t,y)=2\sqrt{y}\) n’est pas lipschitzienne par rapport à \(y\) au voisinage de \((t,0)\) (car \(\partial_y\varphi=1/\sqrt{y}\to\infty\) lorsque \(y\to0\)), donc le théorème d’existence et unicité locale n’est pas valable au voisinage de \((0,0)\). Même raisonnement sur \(\mathbb{R}^-\).

L’EDO est à variables séparables, on peut donc calculer explicitement toutes les solutions du problème de Cauchy :

  • elle admet une seule solution constante, la fonction \(y(t)=0\) pour tout \(t\in\mathbb{R}^+\),
  • et des solutions de la forme \(y(t)=(t+c)^{2}\) pour tout \(t\ge c\).

En imposant la CI on trouve que, pour tout \(b\in\mathbb{R}^+\), les fonctions

\( y_b(t)= \begin{cases} 0 , &\text{si }0\le t\le b, \\ (t-b)^{2} , &\text{si }t\ge b, \end{cases} \)

sont de classe \(\mathcal{C}^1(\mathbb{R}^+)\) et sont solution du problème de Cauchy donné.

Code 
# Remarque : si on demande à sympy de calculer les solutions, il donne une réponse "enigmatique"
# ===============================================================================================

t = sp.Symbol('t')
y = sp.Function('y')

edo = sp.Eq( sp.diff(y(t),t) , 2*sp.sqrt(abs(y(t)))  )
display(Markdown(f"Equation différentielle : ${sp.latex(edo)}$"))

t0, y0 = 0, 0
display(Markdown(f"Condition initiale : $y({t0}) = {y0}$"))

sol = sp.dsolve(edo,y(t), ics={y(t0):y0})
display(Markdown(f"Solution : ${sp.latex(sol)}$"))

Equation différentielle : \(\frac{d}{d t} y{\left(t \right)} = 2 \sqrt{\left|{y{\left(t \right)}}\right|}\)

Condition initiale : \(y(0) = 0\)

Solution : \(\int\limits^{y{\left(t \right)}} \frac{1}{\sqrt{\left|{y}\right|}}\, dy = 2 t + \int\limits^{0} \frac{1}{\sqrt{\left|{y}\right|}}\, dy\)

Traçons quelques courbes :

Code 
xx = np.linspace(0,2,101)
yyr = np.zeros_like(xx)

f = lambda x, b : np.where(x<=b, 0, (x-b)**2)
yy1 = f(xx,1)
yy75= f(xx,0.75)
yy0 = f(xx,0)

plt.figure(figsize=(8,5))
plt.plot(xx,yyr  ,'r-', label=r'$y(t)=0$ $\forall$ $t$')
plt.plot(xx,yy1  ,'b:', label=r'$y_{b=1}(t)$')
plt.plot(xx,yy75 ,'g*', label=r'$y_{b=0.75}(t)$')
plt.plot(xx,yy0  ,'yd', label=r'$y_{b=0}(t)$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');

Code 
# INTERACTIVE PLOT
# =================
# xx = np.linspace(0, 2, 101)
# fig = plt.figure(figsize=(8, 5))
# plt.xlabel('t')
# plt.ylabel('y(t)')
# plt.grid(True)
# plt.legend()
# plt.axis([0, 2, 0, 1.5])

# def plot_solution(b,fig):
#     yy = f(xx, b)
#     plt.plot(xx, yy, 'b-', linewidth=2, label=f'$y_{{b={b:.2f}}}(t)$')
#     plt.title(f'Solution pour b = {b:.2f}')
#     display(fig)
    
# from ipywidgets import interact, FloatSlider
# @interact(b=FloatSlider(min=0, max=1.8, step=0.1, value=0, description='b'))

La méthode d’Euler explicite construit la suite

\(\begin{cases} u_0=y_0=0,\\ u_{n+1}=u_n+h \varphi(t_n,u_n)=u_n+hu_n^{1/2},& n=0,1,2,\dots N_h-1 \end{cases}\)

par conséquent \(u_n=0\) pour tout \(n\): la méthode d’Euler explicite approche la solution constante \(y(t)=0\) pour tout \(t\in\mathbb{R}^+\).

La méthode d’Euler implicite construit la suite

\(\begin{cases} u_0=y_0=0,\\ u_{n+1}=u_n+h \varphi(t_{n+1},u_{n+1})=u_n+hu_{n+1}^{1/2},& n=0,1,2,\dots N_h-1. \end{cases}\)

par conséquent \(u_0=0\) mais \(u_1\) dépend de la méthode de résolution de l’équation implicite \(x=0+h\sqrt{x}\). Bien-sûr \(x=0\) est une solution mais \(x=h^{2}\) est aussi solution. Si le schéma choisit \(u_1=h^{2}\), alors \(u_n>0\) pour tout \(n\in\mathbb{N}^*\).

Notons que le problème de Cauchy avec une CI \(y(0)=y_0>0\) admet une et une seule solution, la fonction \(y(t)=\frac{1}{4}(t+2\sqrt{y_0})^{2}\). Dans ce cas, les deux schémas approchent la même unique solution.

Problème numériquement bien posé : dépendance continue de la solution par rapport aux données.

Une fois calculée la solution numérique \(\{u_n\}_{n=0}^{N}\) d’un problème de Cauchy mathématiquement bien posé, il est légitime de chercher à savoir dans quelle mesure l’erreur \(|y(t_n)-u_n|\) est petite. Cela dépend bien sur du schéma choisi, mais aussi du problème en question. En effet, étant donné que les erreurs d’arrondi amènent toujours à résoudre un problème perturbé, il est important de savoir si la solution d’un problème perturbé est proche de la solution du problème non perturbé.

Définition de problème numériquement bien posé

On dit qu’un problème de Cauchy est numériquement bien posé si la solution d’un problème faiblement perturbé (second membre ou condition initiale) possède une solution proche de celle du problème original.

Si ce n’est pas le cas, on dit qu’il est numériquement mal posé.

Question pratique : Généralement on utilise un schéma pour approcher la solution d’un problème dont on ne connait pas la solution exacte. Comment peut-on conjecturer que le problème est numériquement mal posé ?

Idée : on teste un même schéma avec différentes discrétisations (très différentes entre elles) ou on teste différents schémas. Si les solutions obtenues sont très différentes, le problème est probablement numériquement mal posé.

Exemple de problème BIEN posé numériquement

Soit \(y_0\) donné et considérons une petite perturbation \(\varepsilon\). Calculons \(\lim_{t\to+\infty} y_{\varepsilon}\) où \(y_{\varepsilon}\) est la solution exacte du problème de Cauchy

\(\begin{cases} y_{\varepsilon}'(t) = -y_{\varepsilon}(t), &\forall t>0,\\ y_{\varepsilon}(0) = y_0+\varepsilon. \end{cases}\)

La solution exacte est \(y(t)=y_0e^{-t}+\varepsilon e^{-t}\) : l’effet de la perturbation \(\varepsilon\) s’atténue lorsque \(t\to+\infty\) puisque \(\varepsilon e^{-t}\xrightarrow[t\to+\infty]{}0\).
Cela suggère que si une erreur est faite dans une étape d’une méthode d’itération, l’effet de cette erreur s’atténue au cours du temps: le problème est numériquement bien posé.

Code 
y0 = 1

y = lambda t,eps : (y0+eps)*np.exp(-t)

tt  = np.linspace(0,5,101)
yye = y(tt,0.1) 
yy0 = y(tt,0)   

plt.figure(figsize=(8,5))
plt.plot(tt,yye, label="$y(0)=1.1$")
plt.plot(tt,yy0, label="$y(0)=1$")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');
plt.axis([0,5,0,1.1]);

Exemple de problème BIEN posé numériquement

Considérons le problème de Cauchy \(\begin{cases} y'(t) = y(t)\sin(t), &\forall t>0,\\ y(0) = -1+\varepsilon. \end{cases}\)

La solution est \(y(t)=-e^{1-\cos(t)}+\varepsilon e^{1-\cos(t)}\) : l’effet de la perturbation \(\varepsilon\) reste borné lorsque \(t\to+\infty\) puisque \(|\varepsilon e^{1-\cos(t)}|\le \varepsilon e^2\). Cela suggère que si une erreur est faite dans une étape d’une méthode d’itération, l’effet de cette erreur reste borné au cours du temps: le problème est numériquement bien posé.

Code 
y0 = -1

y = lambda t,eps : (y0+eps)*np.exp(1-np.cos(t))

tt = np.linspace(0,40,401)
yye = y(tt,0.1) 
yy0 = y(tt,0)   

plt.figure(figsize=(8,5)) 
plt.plot(tt,yye,label="$y(0)=-0.9$")
plt.plot(tt,yy0,label="$y(0)=-1$")
plt.legend(bbox_to_anchor=(1,1));

Exemple de problème MAL posé numériquement

Considérons le problème de Cauchy

\(\begin{cases} y'(t) = y(t), &\forall t>0,\\ y(0) = y_0+\varepsilon. \end{cases}\)

La solution est \(y(t)=y_0e^{t}+\varepsilon e^{t}\) : l’effet de la perturbation \(\varepsilon\) s’amplifie lorsque \(t\to+\infty\) puisque \(\varepsilon e^{t}\xrightarrow[t\to+\infty]{}+\infty\). Cela suggère que si une erreur est faite dans une étape d’une méthode d’itération, l’effet de cette erreur s’amplifie au cours du temps: le problème est numériquement mal posé.

Code 
y0 = 0

y = lambda t,eps : (y0+eps)*np.exp(t)

tt  = np.linspace(0,5,101)
yye = y(tt,0.01) 
yy0 = y(tt,0)    

plt.figure(figsize=(8,5)) 
plt.plot(tt, yye, lw='2', label="$y(0)=0.01$")
plt.plot(tt, yy0, lw='5', label="$y(0)=0$")
plt.legend(bbox_to_anchor=(1,1))
plt.axis([0,5,0,1.1]);

Exemple de problème MAL posé numériquement

Soit \(\alpha\) un nombre quelconque. Étudions la solution exacte \(y_{\alpha}\) du problème de Cauchy

\( \begin{cases} y_{\alpha}'(t) = 3y_{\alpha}(t)-3t, &\forall t\in\mathbb{R},\\ y_{\alpha}(0) = \alpha. \end{cases} \)

Comparons ensuite les solutions \(y_{1/3}\) et \(y_{0.333333}\) en \(t=10\).

La solution, définie sur \(\mathbb{R}\), est donnée par

\( y_{\alpha}(t)=\left(\alpha-\frac{1}{3}\right)e^{3t}+t+\frac{1}{3}. \)

Évaluons \(y_{\alpha}\) en \(t=10\):

  • si \(\alpha=\frac{1}{3}\) alors \(y_{\alpha}(10)=\frac{31}{3}\),
  • si \(\alpha=0.333333\) alors \(y_{\alpha}(10)=(0.333333-1/3)e^{30}+10+1/3=-e^{30}/3000000+31/3\approx31/3+10^7/3\).

Conséquence : si nous faisons le calcul avec l’approximation \(\alpha=0.333333\) au lieu de \(1/3\), nous avons \(y(10)=31/3-e^{30}/3000000\) ce qui représente une différence avec la précédente valeur de \(e^{30}/3000000\approx10^7/3\).

Cet exemple nous apprend qu’une petite erreur sur la condition initiale (erreur relative d’ordre \(10^{-6}\)) peut provoquer une très grande erreur sur \(y(10)\) (erreur relative d’ordre \(10^{6}\)). Ainsi, si le calculateur mis à notre disposition ne calcule qu’avec \(6\) chiffres significatifs (en virgule flottante), alors \(\alpha=1/3\) devient \(\alpha=0.333333\) et il est inutile d’essayer d’inventer une méthode numérique pour calculer \(y(10)\). En effet, la seule erreur sur la condition initiale provoque déjà une erreur inadmissible sur la solution. Nous sommes en présence ici d’un problème numériquement mal posé

Vérifions-le avec sympy :

Code 
t = sp.Symbol('t')
y = sp.Function(r'y_{\alpha}')
alpha = sp.Symbol('alpha')

edo = sp.Eq( sp.diff(y(t),t) , 3*y(t)-3*t  )
t0, y0 = 0, alpha
solpar = sp.dsolve(edo,y(t),ics={y(t0):y0})
display(Markdown("Solution : $\\displaystyle "+sp.latex(solpar)+"$" ))

sol1 = solpar.subs(t,10).subs(alpha,sp.Rational(1,3))
sol2 = solpar.subs(t,10).subs(alpha,0.333333)
display(Markdown("Solution avec $\\alpha=\\frac{1}{3}$ : $\\displaystyle "+sp.latex(sol1)+"$" ))
display(Markdown("Solution avec $\\alpha=0.333333$ : $\\displaystyle "+sp.latex(sol2)+"$" ))
display(Markdown("Différence : $\\displaystyle "+sp.latex((sol1.rhs-sol2.rhs).evalf())+"$"))

Solution : \(\displaystyle y_{\alpha}{\left(t \right)} = t + \left(\alpha - \frac{1}{3}\right) e^{3 t} + \frac{1}{3}\)

Solution avec \(\alpha=\frac{1}{3}\) : \(\displaystyle y_{\alpha}{\left(10 \right)} = \frac{31}{3}\)

Solution avec \(\alpha=0.333333\) : \(\displaystyle y_{\alpha}{\left(10 \right)} = \frac{31}{3} - 3.33333333324415 \cdot 10^{-7} e^{30}\)

Différence : \(\displaystyle 3562158.19374618\)

Exercices

La solution exacte est donnée par \(y_\varepsilon(t)=\varepsilon e^{5t}+1\) comme on peut le vérifier avec sympy :

Code 
t = sp.Symbol('t')
epsilon = sp.Symbol(r'\varepsilon')
y = sp.Function(r'y_{\varepsilon}')(t)

edo = sp.Eq(sp.diff(y, t), 5*y - 5)
sol = sp.dsolve(edo, y, ics={y.subs(t, 0): 1 + epsilon})
display(Markdown("Solution : $\\displaystyle "+sp.latex(sol)+"$" ))

# Calcul des solutions pour epsilon=10^(-17) et epsilon=0
sol1 = sol.subs({t: 50, epsilon: sp.Rational(1, 10**17)})
sol0 = sol.subs({t: 50, epsilon: 0})
display(Markdown("Solution avec $\\varepsilon=10^{-17}$ : $\\displaystyle "+sp.latex(sol1)+"$" ))
display(Markdown("Solution avec $\\varepsilon=0$ : $\\displaystyle "+sp.latex(sol0)+"$" ))

# Affichage de la différence entre les solutions
display(Markdown("Différence : $\\displaystyle "+sp.latex((sol1.rhs - sol0.rhs).evalf())+"$"))

Solution : \(\displaystyle y_{\varepsilon}{\left(t \right)} = \varepsilon e^{5 t} + 1\)

Solution avec \(\varepsilon=10^{-17}\) : \(\displaystyle y_{\varepsilon}{\left(50 \right)} = 1 + \frac{e^{250}}{100000000000000000}\)

Solution avec \(\varepsilon=0\) : \(\displaystyle y_{\varepsilon}{\left(50 \right)} = 1\)

Différence : \(\displaystyle 3.74645461450267 \cdot 10^{91}\)

Une très petite erreur sur la condition initiale (de l’ordre de la précision machine) provoque une très grande erreur sur \(y(50)\). Ainsi, il est inutile d’essayer d’inventer une méthode numérique pour calculer \(y(50)\). En effet, la seule erreur sur la condition initiale provoque déjà une erreur inadmissible sur la solution. Nous sommes en présence ici d’un problème numériquement mal posé

Correction

Calculons les solutions de l’edo pour \(t>0\). On pose \(u(t)=\frac{y(t)}{t}\) ainsi \(y(t)=tu(t)\) et \(y'(t)=u(t)+tu'(t)\). On obtient l’edo linéaire d’ordre 1 suivante : \( tu'(t)-2u(t)=-5t^{-3}. \) On a \(a(t)=t\), \(b(t)=-2\) et \(g(t)=-5t^{-3}\) donc
+ \(A(t)=\int \frac{b(t)}{a(t)}\mathrm{d}t=\int -\frac{2}{t}\mathrm{d}t= -2\ln(t)=\ln(t^{-2})\), + \(B(t)=\int \frac{g(t)}{a(t)}e^{A(t)}\mathrm{d}t=\int -5t^{-4} e^{A(t)}\mathrm{d}t=\int -5t^{-6} \mathrm{d}t = t^{-5}\).

On conclut que \(u(t)=ce^{-A(t)}+B(t)e^{-A(t)}=ct^2+t^{-3}\) et \(\boxed{y(t)=tu(t)=ct^3+t^{-2}}\).

En imposant la condition \(a=y(1)\) on trouve l’unique solution du problème de Cauchy donné: \( y(t)=(a-1)t^3+\dfrac{1}{t^2}. \)

Vérifions ce calcul avec sympy:

Code 
t = sp.Symbol('t')
y = sp.Function('y')
a = sp.Symbol('a')

edo = sp.Eq( sp.diff(y(t),t) ,  3*y(t)/t-5/t**3  )
t0, y0 = 1, a
solpar = sp.dsolve(edo,y(t),ics={y(t0):y0})
display(Markdown("Solution : $\\displaystyle "+sp.latex(solpar)+"$" ))

Solution : \(\displaystyle y{\left(t \right)} = \frac{t^{5} \left(a - 1\right) + 1}{t^{2}}\)

Traçons les solutions exactes pour \(a=1\), \(a=1\pm\frac{1}{1000}\) et \(a=1\pm\frac{1}{100}\) sur l’intervalle \([1;5]\) :

Code 
sol_exacte = lambda t,y0 : (y0-1)*t**3+1/t**2

# INITIALISATION
t0     = 1
tfinal = 5
c = 1

tt = np.linspace(t0,tfinal,101)
Y0 = [c+1.e-2,c+1.e-3,c,c-1.e-3,c-1.e-2]

plt.figure(figsize=(8,5))
for y0 in Y0:
    yy = sol_exacte(tt,y0) # [sol_exacte(t,y0) for t in tt] 
    plt.plot(tt, yy, label=f'$y_0={y0}$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid()
plt.title('Solutions exactes pour differents valeurs de $y_0$')
plt.axis([t0,tfinal,-c,c]);

On constate qu’une petite perturbation de la condition initiale entraîne une grosse perturbation de la solution, en effet \( \lim_{t\to+\infty}y(t)= \begin{cases} +\infty &\text{si } a>1,\\ 0^+ &\text{si } a=1,\\ -\infty &\text{si } a<1. \end{cases} \)
Si \(a=1\), la solution exacte est \(y(t)=t^{-2}\) mais le problème est numériquement mal posé car toute (petite) erreur de calcul a le même effet qu’une perturbation de la condition initiale: on “réveille” le terme dormant \(t^{3}\).

Lorsque l’on essaye d’approcher la solution avec les méthodes classiques à un pas, bien sûr lorsque \(N\to+\infty\) on converge vers la solution exacte, cependant en pratique \(N\) est fini et on voit qu’inévitablement nous allons approcher une autre solution que celle correspondant à \(a=1\).

Code 
phi = lambda t, y: 3*y/t - 5/t**3
t0, tfinal, y0 = 1, 5, 1
NN = [20, 100, 300, 5000] # Nombre de points pour les différentes courbes

# --- Définition des "Steps" de chaque schéma ---
def step_EE(phi, t, u, h): return u + h * phi(t, u)
def step_EI(phi, t, u, h): return fsolve(lambda x: u + h * phi(t + h, x) - x, u)[0]
def step_EM(phi, t, u, h): return u + h * phi(t + h/2, u + h/2 * phi(t, u))
def step_HE(phi, t, u, h): 
    u_pred = u + h * phi(t, u)
    return u + h/2 * (phi(t, u) + phi(t + h, u_pred))
def step_CN(phi, t, u, h): 
    return fsolve(lambda x: u + h/2 * (phi(t, u) + phi(t + h, x)) - x, u)[0]

# --- Solveur Générique ---
def solve_ode(method_step, phi, tt, y0):
    uu = np.zeros_like(tt)
    uu[0] = y0
    h = tt[1] - tt[0]
    for i in range(len(tt) - 1):
        uu[i+1] = method_step(phi, tt[i], uu[i], h)
    return uu

# --- Liste des méthodes à tester ---
# (Nom, Fonction de step, Position dans la figure)
methods = [
    ("EE", step_EE, 1),
    ("EI", step_EI, 3),
    ("EM", step_EM, 4),
    ("HE", step_HE, 5),
    ("CN", step_CN, 6)
]

# --- Boucle d'affichage ---
plt.figure(figsize=(15, 10))

for name, step_func, pos in methods:
    plt.subplot(2, 3, pos)
    for N in NN:
        tt = np.linspace(t0, tfinal, N)
        uu = solve_ode(step_func, phi, tt, y0)
        plt.plot(tt, uu, '.-', label=f'N={N}')
    
    # Solution exacte
    t_fine = np.linspace(t0, tfinal, 1000)
    plt.plot(t_fine, sol_exacte(t_fine, y0), 'c-', lw=2, label='Exacte')
    
    plt.title(f'$y_0={y0}$ - Exacte vs {name}')
    plt.legend()
    plt.axis([t0, tfinal, -c, c])

plt.tight_layout()
plt.show()

Problème de Cauchy bien posé mathématiquement et numériquement

Problème bien conditionné (= non raide) : solution approchée en un nombre raisonnable d’itérations.

Définition de problème bien conditionné

Une équation différentielle est dite “stiff” lorsqu’elle contient des composantes qui évoluent à des vitesses très différentes. Par exemple, considérons une équation différentielle d’ordre 2 (un système des deux équations différentielles). Si l’une varie lentement et l’autre chute très rapidement vers un équilibre, alors on dit que le problème est stiff (ou raide). Ou encore on s’interesse à un problème de Cauchy dont la solution admet une composante rapide qui chute très rapidement vers un équilibre et une composante lente qui évolue lentement.

Nous verrons que certains schémas forcent à choisir un pas de temps \(h\) minuscule, calé sur la vitesse de la composante rapide, même si celle-ci a déjà atteint son équilibre et ne “sert plus à rien” pour la dynamique globale.

On dit qu’un problème de Cauchy est bien conditionné si les méthodes numériques usuelles peuvent donner sa solution en un nombre raisonnable d’opérations.

Dans le cas contraire on parle de problème raide.

Exemple de problème MAL conditionné (= raide = stiff)

Considérons le problème de Cauchy avec \(\beta>0\) :

\( \begin{cases} y'(t) = -\beta y(t), &\forall t\in[0;1],\\ y(0) = 1. \end{cases} \)

Ce problème est

  • mathématiquement bien posé pour tout \(\beta\in\mathbb{R}\) : il existe une et une seule solution, elle est définie sur \(\mathbb{R}\) et est donnée par \(y(t)=e^{-\beta t}\);
  • numériquement bien posé pour tout \(\beta>0\) : l’erreur sur la solution sera au pire de l’ordre de l’erreur sur la condition initiale ;
  • mal conditionné ou stiff si \(\beta \gg 1\) : plus \(\beta\) est grand, plus la solution est raide et seul un pas très petit permet d’approcher la solution exacte dans la partie où la solution chute très rapidement vers son équilibre \(0\). Certaines méthodes forcent à choisir un pas de temps \(h\) calé sur la vitesse de la composante rapide.

Nous allons tracer la solution exacte pour \(\beta=1\), \(\beta=10\), \(\beta=50\) et \(\beta=100\) sur l’intervalle \([0;1]\) pour illustrer ce phénomène.
Nous allons aussi étudier analytiquement le comportement du schéma d’Euler explicite pour ce problème.

Voici ci-dessous le tracé de la solution exacte pour \(\beta\) allant de \(1\) à \(100\) :

Code 
exacte = lambda t,b : np.exp(-b*t) 

t0 = 0
tfinal = 1
tt = np.linspace(t0,tfinal,501)

plt.figure(1, figsize=(18, 4))
for i,b in enumerate([1,10,50,100]):
    plt.subplot(1,4,i+1)
    yy = exacte(tt,b) #[exacte(t,b) for t in tt]
    plt.plot(tt, yy)
    plt.title(f'$\\beta={b}$')
    plt.axis([0,1,0,1])
    plt.grid();

La méthode d’Euler explicite pour cette EDO s’écrit

\( u_{n+1}=(1-\beta h)u_n=\dots=(1-\beta h)^{n+1}u_0. \)

La suite obtenue est une suite géométrique de raison \(q=1-\beta h\). On sait qu’une telle suite

  • diverge si \(|q|>1\) ou \(q=-1\),
  • est stationnaire si \(q=1\),
  • converge vers \(0\) si \(|q|<1\) (de façon monotone vers \(0\) si \(0<q<1\)).

ainsi la solution numérique converge vers \(0\) ssi \(0<\beta h<2\) : plus \(\beta\) est grand, plus le pas \(h\) doit être petit pour que la solution numérique puisse approcher la solution exacte.

Intuitivement : la solution exacte \(y(t)=e^{-\beta t}\) décroit très rapidement lorsque \(\beta\) est grand (le problème devient de plus en plus raide) ; pour que la solution numérique puisse suivre cette décroissance rapide, il faut qu’il y ait des points dans cette zone de décroissance rapide, c’est-à-dire que le pas \(h\) soit petit.

On verra plus tard que, si on n’est pas intéressé par cette zone de décroit rapide, mais plutôt par le comportement asymptotique de la solution, on peut utiliser des méthodes numériques spécifiques qui permettent de prendre des pas plus grands. Cela nous amène à la notion d’A-stabilité que nous allons définir dans les prochains chapitres.

Exercices

Considérons l’équation différentielle ordinaire \(y'(t) = ay(t) + b\), où \(a \neq 0\) et \(b\) sont des constantes réelles et considérons une donnée initiale \(y(0) = y_0\neq0\).

On cherche à réécrire cette équation sous la forme d’une équation \(v'(t) = qv(t)\), où \(v\) est une nouvelle fonction à déterminer.

Posons \(v(t) = y(t) + \frac{b}{a}\). Alors \(v(0) = y_0 + \frac{b}{a}\) et

\(\begin{aligned} v'(t) &= y'(t)\\ &= ay(t) + b\\ &=a\left(v(t)-\frac{b}{a}\right)+b\\ &=av(t). \end{aligned}\)

La solution de cette équation est \(v(t) = v(0) e^{at}\).

On conclut alors que \(y(t) = v(t) + \frac{b}{a} = v(0) e^{at} + \frac{b}{a} = \left( y_0 + \frac{b}{a} \right) e^{at} - \frac{b}{a}\), autrement dit

\( \left( y(t)+\frac{b}{a} \right) = \left( y_0+\frac{b}{a} \right) e^{at}. \)

Une suite \((u_n)\) est dite arithmético-géométrique si elle satisfait une relation de la forme \(u_{n+1}=au_n+b\) avec \(a\neq1\) et \(b\) des réels constants.

On cherche à réécrire cette suite sous la forme d’une suite géométrique \((v_n)\) telle que \(v_{n+1}=qv_n\) avec \(q\) une constante réelle.

On pose \(v_n=u_n-\frac{b}{1-a}\). Alors

\(\begin{aligned} v_{n+1} &= u_{n+1}-\frac{b}{1-a}\\ &=au_n+b-\frac{b}{1-a}\\ &=a\left(v_n+\frac{b}{1-a}\right)+b-\frac{b}{1-a}\\ &=av_n+\frac{ab}{1-a}+b-\frac{b}{1-a}\\ &=av_n+\frac{ab+b(1-a)-b}{1-a}\\ &=av_n. \end{aligned}\)

La suite \((v_n)\) est bien géométrique de raison \(q=a\) donc \(v_n = a^n v_0\).

On conclut alors que \(u_n = a^n v_0 + \frac{b}{1-a}= a^n \left(u_0-\frac{b}{1-a}\right) + \frac{b}{1-a}\), autrement dit

\( \left( u_n-\frac{b}{1-a} \right) = \left( u_0-\frac{b}{1-a} \right) a^n . \)

Correction : solution exacte

Il s’agit de l’edo linéaire d’ordre 1 \(y'(t)=-by(t)+g\), dont la solution est

\( y(t)-\frac{g}{b}=\left( y_0-\frac{g}{b} \right)e^{-bt} \)

et

\( \lim_{t\to+\infty}y(t)=\frac{g}{b} \text{ pour tout } y_0\in\mathbb{R} \)

le problème est numériquement bien posé, c’est-à-dire que l’erreur sur la solution sera au pire de l’ordre de l’erreur sur la condition initiale.

Cependant, lorsque \(b\) devient de plus en plus grand, la solution devient de plus en plus raide:

Code 
exacte = lambda t, b, g, y0 : (y0-g/b)*np.exp(-b*t)+g/b 

t0 = 0.
tfinal = 1.0
tt = np.linspace(t0,tfinal,10001)

def plot_exacts(tt, b, g):
    Y0 = [g/b, g/b + 1.e-8]
    for y0 in Y0:
        yy = exacte(tt, b, g, y0)
        plt.plot(tt, yy, label=f'$y_0={y0}$')
    plt.title(f"$g=b={g}$")
    plt.axis([0, 1, g/b - 1.e-8, g/b + 1.e-8])
    plt.legend()
    plt.grid()

plt.figure(1, figsize=(15, 5))

plt.subplot(1, 2, 1)
plot_exacts(tt, b=10, g=10)

plt.subplot(1, 2, 2)
plot_exacts(tt, b=40, g=40)

Correction : Méthode d’Euler explicite

La méthode d’Euler explicite donne la suite définie par récurrence \(u_{n+1}=Au_n+B\) avec \(A=1-bh\) et \(B=gh\), soit encore

\( \left( u_n-\frac{B}{1-A} \right) = \left( u_0-\frac{B}{1-A} \right) A^n \quad\leadsto\quad u_n-\frac{g}{b}=(1-hb)^n\left( y_0 -\frac{g}{b}\right). \)

Comparons la solution exacte et la solution approchée :

\(\begin{aligned} y_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)e^{-bt_n}\\ u_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)(1-hb)^n \end{aligned}\)

Pour que \(u_n\simeq y_n\) il faut que \((1-hb)^n\simeq e^{-bt_n} \in]0,1]\), donc il faut que \(0\le(1-hb)<1\), soit encore \(h<\frac{1}{b}\). Comme \(h=\frac{T-t_0}{N}=\frac{1}{N}\), cela implique \(N\ge b\). Dans cet exemple, pour obtenir une bonne convergence il est nécessaire de prendre un pas \(h\) de pus en plus petit que \(b\) est grand, c’est-à-dire un coût en calcul plus élevé. On parle de problème raide ou mal conditionné.

Code 
phi = lambda t, y, b, g : -b*y+g

def EE(phi,tt,y0,b,g):
    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],b,g)
    return uu

t0 = 0.
tfinal = 1.0

def plot_euler_exp(phi, t0, tfinal, b, g, y0, NN):
    for N in NN:
        tt = np.linspace(t0, tfinal, N+1)
        yy = EE(phi, tt, y0, b, g)
        plt.plot(tt, yy, label=f'$N={N}$ noeuds')
    plt.title(f"Euler explicite, $g=b={g}$, $y_0={y0}$")
    plt.axis([0, 1, g/b - 1.e-8, g/b + 1.e-8])
    plt.legend()
    plt.grid()

g = 40
b = 40
Y0 = [g/b, g/b + 1.e-8]
NN = [19, 23, 45, 100]

plt.figure(figsize=(15, 5))

plt.subplot(1, 2, 1)
plot_euler_exp(phi, t0, tfinal, b, g, y0=Y0[0], NN=NN)

plt.subplot(1, 2, 2)
plot_euler_exp(phi, t0, tfinal, b, g, y0=Y0[1], NN=NN)

Correction : Méthode d’Euler implicite

La méthode d’Euler implicite donne la suite définie par récurrence \(u_{n+1}=Au_n+B\) avec \(A=\frac{1}{1+bh}\) et \(B=\frac{gh}{1+bh}\), soit encore

\( \left( u_n-\frac{B}{1-A} \right) = \left( u_0-\frac{B}{1-A} \right) A^n \quad\leadsto\quad u_n-\frac{g}{b}=\frac{1}{(1+hb)^n}\left( y_0 -\frac{g}{b}\right). \)

Comparons la solution exacte et la solution approchée :

\(\begin{aligned} y_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)e^{-bt_n}\\ u_n&=\frac{g}{b}+\left( y_0 -\frac{g}{b}\right)\left(\frac{1}{1+hb}\right)^n \end{aligned}\)

Pour que \(u_n\simeq y_n\) il faut que \(\frac{1}{(1+hb)^n}\simeq e^{-bt_n} \in]0,1]\), donc il faut que \(0\le\frac{1}{1+hb}<1\) c’est-à-dire pour n’importe quelle valeur de \(h\). Dans cet exemple, pour obtenir la convergence il n’est pas nécessaire de limiter le pas \(h\). Cependant, pour une edo quelconque, à chaque pas de temps on doit résoudre une équation non-linéaire.

Code 
phi = lambda t,y,b,g : -b*y+g

def EI(phi,tt,y0,b,g):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        uu[i+1] = g/b + (y0-g/b)/(1+h*b)**i
        uu[i+1] = fsolve(lambda x: x-uu[i]-h*phi(tt[i+1],x,b,g), uu[i])[0]
    return uu

t0 = 0.
tfinal = 1

def plot_euler_imp(phi, t0, tfinal, b, g, y0, NN, marker='-'):
    for N in NN:
        tt = np.linspace(t0, tfinal, N+1)
        yy = EI(phi, tt, y0, b, g)
        plt.plot(tt, yy, marker, label=f'$N={N}$ noeuds')
    plt.title(f"Euler implicite, $g=b={g}$, $y_0={y0}$")
    plt.axis([0, 1, g/b - 1.e-8, g/b + 1.e-8])
    plt.legend()
    plt.grid()

g = 40
b = 40
Y0 = [g/b, g/b + 1.e-8]
NN = [19, 23, 45, 100]

plt.figure(figsize=(15, 5))

plt.subplot(1, 2, 1)
plot_euler_imp(phi, t0, tfinal, b, g, y0=Y0[0], NN=NN, marker='-')

plt.subplot(1, 2, 2)
plot_euler_imp(phi, t0, tfinal, b, g, y0=Y0[1], NN=NN, marker='-*')

Le problème de Cauchy donné admet une et une seule solution, la fonction constante \(y(t)=\frac{1}{5}\), il est donc bien posé mathématiquement.

Considérons le problème de Cauchy perturbé

\( \begin{cases} y_\varepsilon'(t)=-150y_\varepsilon(t)+30\\ y_\varepsilon(0)=\frac{1}{5}+\varepsilon. \end{cases} \)

La solution exacte \(y_\varepsilon\) est

\( \left( y_\varepsilon(t)+\frac{30}{-150} \right) = \left( y_\varepsilon(0)+\frac{30}{-150} \right) e^{-150t} \quad\leadsto\quad \left( y_\varepsilon(t)-\frac{1}{5} \right) = \left( y_\varepsilon(0)-\frac{1}{5} \right) e^{-150t} = \varepsilon e^{-150t}. \)

Si \(\varepsilon\to0\), alors \(y_{\varepsilon}(t)\to \frac{1}{5}\) : le problème est numériquement bien posé, car l’erreur sur la solution sera inférieure à l’erreur sur la condition initiale.

Code 
# Calcul de la solution exacte
# ============================

t = sp.Symbol('t')
y = sp.Function('y')

eps = sp.Symbol(r'\varepsilon')
t0, y0 = 0, sp.Rational(1,5)+eps
edo = sp.Eq( sp.diff(y(t),t) ,  -150*y(t)+30  )
sol = sp.dsolve(edo,y(t), ics={y(t0):y0})
display(Markdown(f"Solution : ${sp.latex(sol)}$"))

func = sp.lambdify((t,eps),sol.rhs,'numpy')
tt = np.linspace(0,1,101)
yy_eps = func(tt,0.01)
yy_0 = func(tt,0)
plt.figure(figsize=(8,5))
plt.plot(tt,yy_eps,label=r'$\varepsilon=0.01$')
plt.plot(tt,yy_0,label=r'$\varepsilon=0$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid()

Solution : \(y{\left(t \right)} = \varepsilon e^{- 150 t} + \frac{1}{5}\)

La méthode d’Euler implicite donne la suite définie par récurrence \(u_{n+1} = u_n+h(-150u_{n}+30) = au_n+b\) avec \(a=1-150h\) et \(b=30h\), soit encore

\(\begin{aligned} \left( u_n-\frac{b}{1-a} \right) = \left( u_0-\frac{b}{1-a} \right) a^n &\quad\leadsto\quad \left( u_n-\frac{30h}{150h} \right) = \left( u_0-\frac{30h}{150h} \right) (1-150h)^n \\ &\quad\leadsto\quad \left( u_n-\frac{1}{5} \right) = \left( u_0-\frac{1}{5} \right) (1-150h)^n. \end{aligned}\)

Comparons la solution exacte et la solution approchée :

\(\begin{aligned} y_n-\frac{1}{5} &= \varepsilon e^{-150t}\\ u_n-\frac{1}{5} &= \varepsilon (1-150h)^n \end{aligned}\)

Pour que \(u_n\simeq y_n\) il faut que \((1-150h)^n\simeq e^{-150(nh)} \in]0,1]\), donc il faut que \(|1-150h|<1\), c’est-à-dire \(h<\frac{1}{75}\).

Code 
# Calcul de la solution exacte
# ============================

t = sp.Symbol('t')
y = sp.Function('y')

eps = sp.Symbol(r'\varepsilon')
t0, y0 = 0, 0+eps
edo = sp.Eq( sp.diff(y(t),t) ,  -1000*(y(t)-sp.exp(-t))-sp.exp(-t)  )
sol = sp.dsolve(edo,y(t), ics={y(t0):y0})
display(Markdown(f"Solution : ${sp.latex(sol)}$"))

func = sp.lambdify((t,eps),sol.rhs,'numpy')
tt = np.linspace(0,1,501)
yy_1 = func(tt,1)
yy_eps = func(tt,0.01)
yy_0 = func(tt,0)
plt.figure(figsize=(8,5))
plt.plot(tt,yy_1,label=r'$\varepsilon=1$')
plt.plot(tt,yy_eps,'.',label=r'$\varepsilon=0.01$')
plt.plot(tt,yy_0,':',label=r'$\varepsilon=0$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid()

Solution : \(y{\left(t \right)} = \left(\left(\varepsilon - 1\right) e^{- 999 t} + 1\right) e^{- t}\)

Le problème est

  • bien posé mathématiquement : il admet une et une seule solution, la fonction \(y(t)=e^{-t}\), définie sur \(\mathbb{R}\),
  • bien posé numériquement : l’erreur \(\varepsilon\) sur la donnée initiale est atténuée au cours du temps,
  • raide : la solution exacte \(y(t)=e^{-t}-e^{-1000t}\) croit très rapidement (\(-e^{-1000t}\)) sur un très petit intervalle, puis elle décroit comme \(e^{-t}\). Notons que, si \(\varepsilon=1\), la solution est \(y(t)=e^{-t}\).
AstucePropriétés du problème de Cauchy : que retenir ?
  • Un problème de Cauchy est mathématiquement bien posé s’il existe une et une seule solution. Si ce n’est pas le cas, inutile de chercher à approcher la solution numériquement.
  • Un problème de Cauchy est numériquement bien posé si la solution d’un problème faiblement perturbé (second membre ou condition initiale) possède une solution proche de celle du problème original. Si ce n’est pas le cas, on peut chercher à approcher la solution numériquement, mais c’est délicat, car une petite erreur peut provoquer une très grande erreur sur la solution. Il faudra alors utiliser des méthodes numériques spécifiques (d’ordre élevé).
  • Un problème de Cauchy est bien conditionné si les méthodes numériques usuelles peuvent donner sa solution en un nombre raisonnable d’opérations. Dans le cas contraire on parle de problème raide (stiff). Si un problème est raide, il faudra utiliser des méthodes numériques spécifiques (généralement des méthodes A-stable car on verra que c’est grosso modo le même problème que pour les problèmes numériquement mal posés).