On considère le système différentiel linéaire homogène
\( \mathbf{y}'(t)=-\mathbb{A}\mathbf{y}(t), \qquad t\ge 0, \)
où \(\mathbb{A}\in\mathbb{R}^{n\times n}\) est une matrice constante. On suppose que \(\mathbb{A}\) est diagonalisable et que ses valeurs propres \(\lambda_1,\dots,\lambda_n\) sont réelles strictement positives. Dans ce cas, la solution générale s’écrit
\( \mathbf{y}(t)=\sum_{i=1}^{n} c_i\,\mathbf{v}_i\,e^{-\lambda_i t}, \)
où \(\mathbf{v}_i\) est un vecteur propre associé à \(\lambda_i\) et les constantes \(c_i\) sont fixés par la condition initiale.
Chaque mode exponentiel \(e^{-\lambda_i t}\) est associé à une échelle de temps caractéristique
\( \tau_i = \frac{1}{\lambda_i}. \)
Cette quantité mesure le temps sur lequel la composante correspondante décroît significativement :
La solution peut donc contenir simultanément plusieurs échelles de temps.
On dit que le système est raide lorsque ces échelles sont très différentes, c’est-à-dire lorsque \(\lambda_{\max} \gg \lambda_{\min}\), ou équivalemment lorsque le rapport de raideur
\( \kappa = \frac{\lambda_{\max}}{\lambda_{\min}} \)
est très grand.
Dans un problème raide, la solution présente généralement deux régimes : un régime transitoire rapide, gouverné par les grandes valeurs propres ; un régime lent, gouverné par les petites valeurs propres. Les composantes rapides disparaissent après un temps très court, puis la dynamique devient essentiellement lente. C’est cette coexistence de dynamiques rapides et lentes qui définit la raideur.
Pour un schéma explicite appliqué au système précédent, la stabilité impose typiquement :
\( h \lambda_{\max} \le C, \)
où \(h\) est le pas de temps et \(C\) une constante dépendant du schéma. Ainsi, le pas admissible est déterminé par la dynamique la plus rapide. Après disparition du transitoire rapide, cette contrainte reste pourtant active, ce qui entraîne un coût de calcul élevé.
Considérons l’équation scalaire \(y'(t) = \varphi(t, y(t))\). Même si la solution \(y(t)\) varie peu, la raideur peut apparaître si la fonction \(\varphi\) est très sensible à de petites variations de \(y\), autrement dit lorsque la dérivée partielle
\( \left| \frac{\partial \varphi}{\partial y} (t, y) \right| \)
est grande, même si la solution \(y(t)\) varie peu sur de grands intervalles.
Dans la suite, nous adoptons la convention de J.-P. Demailly.
Définition (Jean-Pierre Demailly) :
un problème est raide s’il impose à une méthode numérique d’utiliser un pas de temps excessivement petit, même lorsque la solution globale varie peu.
Idée essentielle à retenir : la raideur n’est pas une difficulté analytique, mais une inadéquation entre les échelles de temps du problème et la région de stabilité de la méthode numérique. Les problèmes raides nécessitent souvent des méthodes implicites, permettant l’usage de pas de temps plus grands tout en conservant la stabilité.
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
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\) :
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
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. Cependant, 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, c’est les méthodes inconditionnellement A-stables qu’on a vu au chapitre précédent.
La A-stabilité, bien qu’elle soit une propriété définie par le comportement asymptotique (\(t \to \infty\)), est l’outil fondamental pour approcher la solution des problèmes raides sur des intervalles de temps borné. Nous pouvons voir cette équivalence comme un changement d’échelle temporelle.
La raideur comme dilatation du temps : considérons le problème de Cauchy “stiff” sur l’intervalle borné \(I = [0, 1]\) avec \(\gamma \gg 1\) :
\( \begin{cases} y'(t) = -\gamma y(t), & t \in [0, 1], \\ y(0) = 1. \end{cases} \)
Pour analyser le comportement d’un schéma numérique à l’échelle de la variation réelle de la solution, effectuons le changement de variable (dilatation du temps) suivant : \(\tau = \gamma t\) et posons \(\tilde y(\tau) = y(t=\frac{\tau}{\gamma})\).
Par dérivation composée, on obtient \(\tilde y'(\tau) = y'(t) \cdot t'(\tau) = - \gamma y(\frac{\tau}{\gamma}) \cdot \frac{1}{\gamma} = -\tilde y(\tau)\). Le problème se reformule alors en termes de \(\tilde y\) et \(\tau\) :
\( \begin{cases} \tilde y'(\tau) = -\tilde y(\tau), & \tau \in [0, \gamma], \\ \tilde y(0) = 1. \end{cases} \)
Du problème local au comportement asymptotique : ce changement de variable révèle une équivalence qui lie la raideur à la stabilité.
Pour intégrer numériquement sur un intervalle aussi grand sans que les erreurs de discrétisation ne s’accumulent ou n’explosent, le schéma doit simuler correctement le comportement de la solution pour \(\tau \to \infty\). C’est précisément la définition de la A-stabilité.
En résumé : la raideur d’une EDO contraint le schéma numérique à “résoudre” des dynamiques internes extrêmement rapides, même si l’intervalle de temps global est court. Sans A-stabilité inconditionnelle, le pas de temps \(h\) est dicté par la stabilité (il doit être inversement proportionnel à \(\gamma\)), ce qui rend le calcul prohibitivement long. Avec la A-stabilité, on s’affranchit de cette contrainte. Le pas de temps peut être choisi librement en fonction de la précision souhaitée, permettant ainsi de traverser efficacement les zones de forte raideur.
Exercice : résolution numérique d’un problème raide
Soient \(b,g>0\) et considérons le problème de Cauchy
\(\begin{cases} y'(t)=-by(t)+g,& t\in[0;1]\\ y(0)=y_0. \end{cases}\)
Solution exacte
Méthode d’Euler explicite
Méthode d’Euler implicite
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:
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é.
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.
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='-*')
Exercice : problème MAL conditionné (=raide)
\( \begin{cases} y'(t)=-150y(t)+30\\ y(0)=\frac{1}{5}. \end{cases} \)
Vérifier qu’il est bien posé mathématiquement et numériquement.
Que se passe-t-il si on utilise la méthode d’Euler explicite ?
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.
# 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}\).
Exercice : problème raide
Il existe des problèmes qui résistent à toutes les méthodes usuelles. On trouve par exemple ces phénomènes en chimie où deux réactions ont des échelles de temps très différentes. Considérons par exemple \( \begin{cases} y'(t)=-1000\big(y(t)-e^{-t}\big)-e^{-t}, & t\in[0;1]\\ y(0)=0. \end{cases} \)
Vérifier, en calculant la solution exacte, que le problème est bien posé mathématiquement et numériquement et qu’il est raide.
# 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