Propriétés du problème de Cauchy : bien posés (mathématiquement et numériquement), bien conditionnés – M64 – Data Science & Scientific Computing 6 : Schémas numériques pour les problèmes de Cauchy

Propriétés du problème de Cauchy : bien posés (mathématiquement et numériquement), bien conditionnés

Auteur·rice

Gloria Faccanoni

Date de publication

16 février 2026

Code 
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 12})

from scipy.optimize import fsolve

import sympy as sym
sym.init_printing()

from IPython.display import Markdown

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é.

Que se passe-t-il lorsqu’on utilise un schéma pour un 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é.

Dans ce cas, 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.

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

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

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

solgenLIST = sym.dsolve(edo,y(t), ics={y(0):0} )
display(Markdown(f"Solution(s) : ${sym.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]\)

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 = sym.Symbol('t')
y = sym.Function('y')


display(Markdown("**Problème de Cauchy**"))
edo = sym.Eq( sym.diff(y(t),t) , 2*sym.sqrt(abs(y(t)))  )
display(edo)

t0, y0 = 0, 0
display(sym.Eq(y(t0),y0))

# Méthode directe
display(Markdown("**Méthode directe**"))
sol = sym.dsolve(edo,y(t), ics={y(t0):y0})
display(sol)

# Méthode pas à pas
display(Markdown("**Méthode pas à pas**"))
solGEN = sym.dsolve(edo,y(t))
display(solGEN)
consts = sym.solve( [solGEN.rhs.subs(t,t0)-y0 ], dict=True)[0]
solPAR = solGEN.subs(consts).simplify()
display(solPAR)

Problème de Cauchy

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

\(\displaystyle y{\left(0 \right)} = 0\)

Méthode directe

\(\displaystyle \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\)

Méthode pas à pas

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

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

Traçons quelques courbes :

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

# f = lambda x,b : 0 if (x<=b) else (x-b)**2
# yy1 = [f(x,1)    for x in xx]
# yy75= [f(x,0.75) for x in xx]
# yy0 = [f(x,0)    for x in 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();

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) # [y(t,0.1) for t in tt]
yy0 = y(tt,0)   # [y(t,0) for t in tt]

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,1));
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) # [y(t,0.1) for t in tt]
yy0 = y(tt,0)   # [y(t,0) for t in tt]

plt.figure(figsize=(10,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) # [y(t,0.01) for t in tt]
yy0 = y(tt,0)    # [y(t,0) for t in tt]

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 = sym.Symbol('t')
y = sym.Function(r'y_{\alpha}')

edo = sym.Eq( sym.diff(y(t),t) , 3*y(t)-3*t  )
#display(edo)

solgen = sym.dsolve(edo,y(t))
#display(solgen)

alpha = sym.Symbol('alpha')
t0, y0 = 0, alpha
consts = sym.solve( [solgen.rhs.subs(t,t0)-y0 ], dict=True)[0]
solpar = solgen.subs(consts).simplify()
display(solpar)

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

\(\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 
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

# Définition des symboles
t = sym.Symbol('t')
epsilon = sym.Symbol(r'\varepsilon')

# Définition de la fonction y(t)
y = sym.Function(r'y_{\varepsilon}')(t)

# Équation différentielle
edo = sym.Eq(sym.diff(y, t), 5*y - 5)
display(edo)

# Condition initiale
display( sym.Eq(y.subs(t, 0), 1 + epsilon))

# Résolution
sol = sym.dsolve(edo, y, ics={y.subs(t, 0): 1 + epsilon})
display(sol)

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

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

\(\displaystyle \frac{d}{d t} y_{\varepsilon}{\left(t \right)} = 5 y_{\varepsilon}{\left(t \right)} - 5\)

\(\displaystyle y_{\varepsilon}{\left(0 \right)} = \varepsilon + 1\)

\(\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 = sym.Symbol('t')
y = sym.Function('y')

edo = sym.Eq( sym.diff(y(t),t) ,  3*y(t)/t-5/t**3  )
display(edo)

solgen = sym.dsolve(edo,y(t))
display(solgen)

a = sym.Symbol('a')
t0, y0 = 1, a
consts = sym.solve( [solgen.rhs.subs(t,t0)-y0 ], dict=True)[0]
solpar = solgen.subs(consts)
display(solpar)

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = \frac{3 y{\left(t \right)}}{t} - \frac{5}{t^{3}}\)

\(\displaystyle y{\left(t \right)} = \frac{C_{1} t^{5} + 1}{t^{2}}\)

\(\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=(10,7))
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     = 1
tfinal = 5
y0     = 1


def EE(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros(len(tt))
    uu[0] = y0
    for i in range(len(tt)-1):
        uu[i+1] = uu[i]+h*phi(tt[i],uu[i])
    return uu

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

def EM(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros(len(tt))
    uu[0] = y0
    for i in range(len(tt)-1):
        utilde = uu[i]+h/2*phi(tt[i],uu[i])
        uu[i+1] = uu[i]+h*phi(tt[i]+h/2,utilde)
    return uu

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

def HE(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros(len(tt))
    uu[0] = y0
    for i in range(len(tt)-1):
        utilde = uu[i]+h*phi(tt[i],uu[i])
        uu[i+1] = uu[i] + h/2 * (phi(tt[i],uu[i])+phi(tt[i+1],utilde))
    return uu


NN = [20,100,300,5000]

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

plt.subplot(2,3,1)
for N in NN:
    tt = np.linspace(t0,tfinal,N)
    uu = EE(phi,tt,y0)
    plt.plot(tt,uu,'.-')
plt.plot(tt, sol_exacte(tt,y0),'c-',lw=2)
plt.legend(['N=20','N=100','N=300','N=5000','Exacte'])
plt.title('$y_0=1$ - Exacte vs EE')
plt.axis([t0,tfinal,-c,c]);

plt.subplot(2,3,3)
for N in NN:
    tt = np.linspace(t0,tfinal,N)
    uu = EI(phi,tt,y0)
    plt.plot(tt,uu,'.-')
plt.plot(tt,[sol_exacte(t,y0) for t in tt],'c-',lw=2)
plt.legend(['N=20','N=100','N=300','N=5000','Exacte'])
plt.title('$y_0=1$ - Exacte vs EI')
plt.axis([t0,tfinal,-c,c]);

plt.subplot(2,3,4)
for N in NN:
    tt = np.linspace(t0,tfinal,N)
    uu = EM(phi,tt,y0)
    plt.plot(tt,uu,'.-')
plt.plot(tt,[sol_exacte(t,y0) for t in tt],'c-',lw=2)
plt.legend(['N=20','N=100','N=300','N=5000','Exacte'])
plt.title('$y_0=1$ - Exacte vs EM')
plt.axis([t0,tfinal,-c,c]);

plt.subplot(2,3,5)
for N in NN:
    tt = np.linspace(t0,tfinal,N)
    uu = HE(phi,tt,y0)
    plt.plot(tt,uu,'.-')
plt.plot(tt,[sol_exacte(t,y0) for t in tt],'c-',lw=2)
plt.legend(['N=20','N=100','N=300','N=5000','Exacte'])
plt.title('$y_0=1$ - Exacte vs HE')
plt.axis([t0,tfinal,-c,c]);

plt.subplot(2,3,6)
for N in NN:
    tt = np.linspace(t0,tfinal,N)
    uu = CN(phi,tt,y0)
    plt.plot(tt,uu,'.-')
plt.plot(tt,[sol_exacte(t,y0) for t in tt],'c-',lw=2)
plt.legend(['N=20','N=100','N=300','N=5000','Exacte'])
plt.title('$y_0=1$ - Exacte vs CN')
plt.axis([t0,tfinal,-c,c]);

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

Définition de problème bien conditionné

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 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 le pas \(h\) devient petit car la solution devient de plus en plus raide (stiff) et seul un pas très petit permet d’approcher la solution exacte indépendamment de la méthode numérique utilisée.

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'$b={b}$')
    plt.axis([0,1,0,1])
    plt.grid();

Intuitivement : la solution exacte \(y(t)=e^{-\beta t}\) décroit très rapidement lorsque \(\beta\) est grand ; 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.

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\). On voit donc que plus \(\beta\) est grand, plus le pas \(h\) doit être petit pour que la solution numérique puisse approcher la solution exacte. Cela signifie que le problème devient de plus en plus raide (stiff).

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)

plt.figure(1, figsize=(15, 5))
plt.subplot(1,2,1)
g = 10
b = 10
Y0 = [g/b, g/b+1.e-8]
for y0 in Y0:
    yy = exacte(tt,b,g,y0) #  [exacte(t,b,g,y0) for t in tt]
    plt.plot(tt,yy,label=('$y_0=$'+str(y0)));
plt.title("$g=b=$"+str(g))
plt.axis([0,1,g/b-1.e-8,g/b+1.e-8])
plt.legend()
plt.grid()

plt.subplot(1,2,2)
g = 40
b = 40
Y0 = [g/b, g/b+1.e-8]
for y0 in Y0:
    yy = exacte(tt,b,g,y0) #  [exacte(t,b,g,y0) for t in tt]
    plt.plot(tt,yy,label=(r'$y_0=$'+str(y0)));
plt.title("$g=b=$"+str(g))
plt.axis([0,1,g/b-1.e-8,g/b+1.e-8])
plt.legend()
plt.grid();

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

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

plt.figure(figsize=(18, 7))
plt.subplot(1,2,1)
y0 = Y0[0]
for N in NN:
    tt = np.linspace(t0,tfinal,N+1)
    yy = EE(phi,tt,y0,b,g)
    plt.plot(tt,yy,label=('$N=$'+str(N)+' noeuds'));
plt.title("Euler explicite, $g=b=$"+str(g)+ ", $y_0$="+str(y0))
plt.axis([0,1,g/b-1.e-8,g/b+1.e-8])
plt.legend()
plt.grid()

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

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

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

plt.figure(figsize=(18,7))
plt.subplot(1,2,1)
y0 = Y0[0]
for N in NN:
    tt = np.linspace(t0,tfinal,N+1)
    yy = EI(phi,tt,y0,b,g)
    plt.plot(tt,yy,label=('$N=$'+str(N)+' noeuds'));
plt.title("Euler implicite, $g=b=$"+str(g)+ ", $y_0$="+str(y0))
plt.axis([0,1,g/b-1.e-8,g/b+1.e-8])
plt.legend()
plt.grid()

plt.subplot(1,2,2)
y0 = Y0[1]
for N in NN:
    tt = np.linspace(t0,tfinal,N+1)
    yy = EI(phi,tt,y0,b,g)
    plt.plot(tt,yy,'-*',label=('$N=$'+str(N)+' noeuds'));
plt.title("Euler implicite, $g=b=$"+str(g)+ ", $y_0$="+str(y0))
plt.axis([0,1,g/b-1.e-8,g/b+1.e-8])
plt.legend()
plt.grid()

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 = sym.Symbol('t')
y = sym.Function('y')

edo = sym.Eq( sym.diff(y(t),t) ,  -150*y(t)+30  )
display(edo)
eps = sym.Symbol(r'\varepsilon')
t0, y0 = 0, sym.Rational(1,5)+eps
display(sym.Eq(y(t0),y0))

sol = sym.dsolve(edo,y(t), ics={y(t0):y0})
display(sol)

func = sym.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()
plt.grid()

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = 30 - 150 y{\left(t \right)}\)

\(\displaystyle y{\left(0 \right)} = \varepsilon + \frac{1}{5}\)

\(\displaystyle 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 = sym.Symbol('t')
y = sym.Function('y')

edo = sym.Eq( sym.diff(y(t),t) ,  -1000*(y(t)-sym.exp(-t))-sym.exp(-t)  )
display(edo)
eps = sym.Symbol(r'\varepsilon')
t0, y0 = 0, 0+eps
display(sym.Eq(y(t0),y0))

sol = sym.dsolve(edo,y(t), ics={y(t0):y0})
display(sol.expand())

func = sym.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()
plt.grid()

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = - 1000 y{\left(t \right)} + 999 e^{- t}\)

\(\displaystyle y{\left(0 \right)} = \varepsilon\)

\(\displaystyle y{\left(t \right)} = \varepsilon e^{- 1000 t} + e^{- t} - e^{- 1000 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).