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})Propriétés du problème de Cauchy
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
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\).
Exercice : problème MAL posé mathématiquement
- Montrer que le problème de Cauchy suivant admet une infinité de solutions de classe \(\mathcal{C}^1(\mathbb{R})\):
\( \begin{cases} y'(t)=2\sqrt{|y(t)|}, & t>0\\ y(0)=0. \end{cases}\)
- Parmi ces solutions, quelle solution approche-t-on si on utilise la méthode d’Euler explicite? et la méthode d’Euler implicite?
- Que se passe-t-il si la donnée initiale est \(y(0)=y_0>0\)?
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é.
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.
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\)
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
Exercice : problème numériquement MAL posé à la précision machine
On considère un problème de Cauchy dépendant d’un paramètre \(\varepsilon\) et on note \(y_\varepsilon(t)\) sa solution. Estimer la différence en \(t=50\) entre les solutions \(y_{\varepsilon}\) avec \(\varepsilon=10^{-17}\) et \(\varepsilon=0\) :
\( \begin{cases} y_\varepsilon'(t)=5y_\varepsilon(t)-5, &t\in[0;50]\\ y_\varepsilon(0)=1+\varepsilon. \end{cases} \)
Notons que \(\varepsilon=10^{-17}\) est une valeur très proche de \(0\) et correspond à l’erreur d’arrondi d’un calculateur qui travaille avec \(17\) chiffres significatifs.
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é
Exercice : problème numériquement MAL posé
L’objectif de cet exercice est de bien mettre en évidence l’influence d’une petite perturbation de la condition initiale sur la solution d’un problème numériquement mal posé et d’appliquer la méthode d’Euler explicite pour en approcher la solution exacte.
Considérons le problème de Cauchy
\(\begin{cases} y'(t)=3\dfrac{y(t)}{t}-\dfrac{5}{t^3},\\ y(1)=a \end{cases}\)
- Calculer l’unique solution exacte à ce problème de Cauchy en fonction de \(a\) pour \(t>0\). Remarquer qu’avec le changement de variable \(u(t)=\frac{y(t)}{t}\) on obtient une EDO linéaire d’ordre 1.
- Tracer la solution exacte sur l’intervalle \([1;5]\) pour \(a=1\), \(a=1\pm\frac{1}{1000}\) et \(a=1\pm\frac{1}{100}\).
- Pour \(a=1\), sur un autre graphique tracer la solution exacte et les solutions obtenues par la méthode d’Euler explicite avec \(N=20\), \(N=100\), \(N=300\), \(N=5000\). Comment expliquez-vous ce comportement du schéma?
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()
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é).