#+OPTIONS:   H:3 num:nil toc:nil \n:nil @:t ::t |:t ^:nil -:t f:t *:t TeX:t LaTeX:t skip:t d:t tags:not-in-toc creator:t timestamp:nil author:nil title:nil html5-fancy:t
#+HTML_DOCTYPE: html5
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../css/stylesheet-blog.css">
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../css/theorems-fr.css">
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../fontawesome/css/all.css" />
#+HTML_HEAD: <link rel="stylesheet" href="../../highlight/styles/base16/solarized-light.min.css">
#+HTML_HEAD: <script src="../../highlight/highlight.min.js"></script>
#+HTML_HEAD: <script>hljs.highlightAll();</script>
#+LANGUAGE: fr
#+STARTUP: latexpreview
#+TITLE: Équations différentielles du premier ordre
#+DATE: Septembre 2026
#+AUTHOR: Stéphane Adjemian
#+EMAIL: stephane.adjemian@univ-lemans.fr
#+PROPERTY: header-args:python :python /tmp/blog-edo1/bin/python

#+BEGIN_QUOTE
Une équation différentielle du premier ordre est une équation fonctionnelle où
n'apparaissent que la fonction inconnue et sa dérivée première. On commence
par les problèmes linéaires, pour lesquels un traitement analytique est toujours
possible, avant d'aborder les équations non linéaires. Dans ce dernier cas on
ne sait pas toujours exhiber une solution, et il faut recourir à l'astuce (un
changement de variable peut rendre le problème linéaire), à l'approximation
(un développement de Taylor), à l'analyse graphique (pour déduire des
propriétés de la solution sans la calculer) ou à l'ordinateur. Les
calculs numériques sont faits en Python, et les corrigés des exercices sont
donnés dans des blocs dépliables. La note se termine par une application au modèle de Solow, où l'on obtient la solution
explicite de l'approximation à l'ordre deux de la dynamique de transition.
#+END_QUOTE

\\
\\
\\

#+BEGIN_SRC bash :results silent :exports none :async t
  python3 -m venv /tmp/blog-edo1
  source /tmp/blog-edo1/bin/activate
  pip install numpy scipy matplotlib
#+END_SRC

#+begin_src python :session edo1 :exports none :results none
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.integrate import solve_ivp
#+end_src

* Équations linéaires à coefficients constants
:PROPERTIES:
:CUSTOM_ID: constants
:END:

Dans cette section, nous considérons une équation différentielle de la
forme\nbsp{}:

\[
\dot y(t) + a\,y(t) = b
\]

où $a$ et $b$ sont des paramètres réels et où $\dot y$ désigne la dérivée de
$y$ par rapport au temps. Nous distinguons deux cas selon que $b$ est nul ou
non.

** Le cas homogène
:PROPERTIES:
:CUSTOM_ID: homogene
:END:

Nous commençons par le cas le plus simple, où la constante $b$ est nulle\nbsp{}:

\[
\dot y(t) + a\,y(t) = 0
\]

On dit qu'une équation différentielle est homogène lorsque la dynamique n'est
pas affectée si l'on multiplie $y(t)$ et $\dot y(t)$ par une même constante.
On peut écrire cette équation de façon équivalente\nbsp{}:

\[
\frac{\dot y(t)}{y(t)} = -a
\]

La dynamique correspond donc à une hypothèse de croissance à taux constant.
Pour résoudre l'équation différentielle, c'est-à-dire exhiber une expression
de $y$ en fonction du paramètre $a$ et de $t$, il suffit de savoir dériver la
fonction logarithme. En effet, puisque $\frac{\mathrm d}{\mathrm dt}\log u(t) =
\dot u(t)/u(t)$, l'équation précédente s'écrit aussi\nbsp{}:

\[
\frac{\mathrm d}{\mathrm dt}\log y(t) = -a
\]

Cette équation nous dit que les variations du logarithme de $y$ sont
constantes. Pour obtenir $\log y(T)$ nous sommons ces variations entre $0$ et
$T$, c'est-à-dire que nous intégrons les deux membres\nbsp{}:

\begin{equation*}
\begin{split}
\int_0^T \frac{\mathrm d}{\mathrm dt}\log y(t)\,\mathrm dt = -\int_0^T a\,\mathrm dt
&\Leftrightarrow \bigl[\log y(t)\bigr]_0^T = -aT\\
&\Leftrightarrow \log y(T) - \log y(0) = -aT\\
&\Leftrightarrow \frac{y(T)}{y(0)} = e^{-aT}
\end{split}
\end{equation*}

où la dernière équivalence est obtenue en appliquant la fonction réciproque du
logarithme népérien, l'exponentielle. Cette transformation permet de revenir
à la variable qui nous intéresse, $y(t)$ et non $\log y(t)$. Finalement\nbsp{}:

\[
y(t) = y(0)\,e^{-at} \qquad \forall t\in\mathbb R^+
\]

Il s'agit de la solution de l'équation différentielle linéaire du premier
ordre homogène à coefficient constant. Elle dépend de la condition initiale,
de la variable $t$ et du paramètre $a$.

#+BEGIN_remarque
Si la condition initiale est nulle, alors $y(t)=0$ pour tout $t$ dans
$\mathbb R^+$. On dit que zéro est l'/état stationnaire/, ou le /point fixe/,
de cette dynamique.
#+END_remarque

#+BEGIN_remarque
Les propriétés de la dynamique sont liées au signe du paramètre $a$. Si $a>0$
alors, pour toute condition initiale $y(0)$, \(\lim_{t\to\infty}y(t)=0\)\nbsp{}: la
dynamique est /stable/, à long terme la variable $y$ rejoint l'état
stationnaire, et le niveau de long terme ne dépend pas de la condition
initiale. À l'inverse, si $a < 0$ alors, pour toute condition initiale
$y(0)\neq0$, \(\lim_{t\to\infty}|y(t)|=\infty\)\nbsp{}: la variable /diverge/, vers
$-\infty$ ou $+\infty$ selon le signe de la condition initiale. La figure
[[fig:homogene][ci-dessous]] illustre ces différents cas.
#+END_remarque

#+begin_src python :session edo1 :exports none :results none
  fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
  t = np.linspace(0, 6, 300)
  for ax, a in ((ax1, 0.5), (ax2, -0.5)):
      for y0 in (-2, -1, -0.5, 0.5, 1, 2):
          ax.plot(t, y0*np.exp(-a*t), 'b', linewidth=1)
      ax.axhline(y=0, color='r', linewidth=1, linestyle='--')
      ax.set_xlabel(r'$t$')
      ax.set_ylabel(r'$y(t)$')
      ax.set_ylim(-4, 4)
  ax1.set_title(r'$a = 0{,}5 > 0$ : dynamique stable')
  ax2.set_title(r'$a = -0{,}5 < 0$ : dynamique instable')
  fig.tight_layout()
  fig.savefig("edo-homogene.svg", transparent=True)
#+end_src

#+CAPTION: *Solutions de l'équation homogène pour plusieurs conditions initiales. À gauche les trajectoires convergent vers l'état stationnaire, à droite elles s'en éloignent.*
#+LABEL: fig:homogene
[[file:edo-homogene.svg]]

** Le cas non homogène
:PROPERTIES:
:CUSTOM_ID: non-homogene
:END:

On s'intéresse maintenant au cas où $b\neq0$, c'est-à-dire à un problème de la
forme\nbsp{}:

\[
\dot y(t) + a\,y(t) = b
\]

La résolution s'obtient en quatre étapes.

*Étape 1. La solution générale de l'équation homogène associée.* L'équation
sans second membre associée est $\dot y(t) = -a\,y(t)$. Nous l'avons déjà
rencontrée dans la section [[#homogene][précédente]]\nbsp{}; nous allons chercher sa solution
générale en suivant une démarche légèrement différente. En supposant $y\neq0$,
nous avons de façon équivalente $\frac{\mathrm d}{\mathrm dt}\log y(t) = -a$. La
dérivée par rapport à $t$ du logarithme de $y(t)$ est égale à $-a$ pour tout
$t$, donc les primitives des deux membres coïncident à une constante
près\nbsp{}:

\[
\log y(t) + \gamma_1 = -at + \gamma_2
\]

ou encore $\log y(t) = -at + \gamma$, puisque nous ne pourrons pas identifier
séparément les deux constantes $\gamma_1$ et $\gamma_2$ (nous comprendrons
plus loin pourquoi). En appliquant l'exponentielle, on obtient
finalement\nbsp{}:

\[
y_1(t) = \Gamma e^{-at}
\]

où $\Gamma = e^{\gamma}$. Cette équation définit un continuum de fonctions
indexées par la constante \(\Gamma\)\nbsp{}; chacune est solution de l'équation
sans second membre. C'est en ce sens qu'il s'agit d'une solution /générale/. À
ce stade, la constante $\Gamma$ n'est pas déterminée.

*Étape 2. Une solution particulière de l'équation complète.* Nous choisissons
la plus simple possible, la solution constante\nbsp{}:

\[
y_2(t) = \frac{b}{a} \qquad \forall t
\]

Il s'agit de l'état stationnaire de l'équation complète\nbsp{}: si $y(t)=b/a$
alors $y(t+\Delta)=b/a$ pour tout $\Delta>0$. Notons en passant que nous
avons implicitement supposé $a\neq0$. Si cette hypothèse n'est pas satisfaite,
une fonction constante n'est pas une solution particulière et il faut en
chercher une autre (voir l'exercice [[ex-a-nul][ci-dessous]]).

*Étape 3. La solution générale de l'équation complète.* En sommant la solution
générale de l'équation sans second membre et une solution particulière de
l'équation complète, nous obtenons la solution générale de l'équation
complète\nbsp{}:

\[
y(t) = y_1(t) + y_2(t) = \Gamma e^{-at} + \frac{b}{a}
\]

Cette solution est générale en ce que le paramètre $\Gamma$ n'est toujours pas
déterminé. Elle fixe la forme de la fonction solution, mais pas son niveau.
Choisir une valeur pour $\Gamma$ revient à choisir une trajectoire parmi une
infinité de trajectoires envisageables.

*Étape 4. Solution de l'équation complète.* Choisir $\Gamma$ est aisé si nous
connaissons la valeur de $y$ à un instant quelconque. Par exemple, nous
pouvons connaître la condition initiale $y(0)=y_0$ (on parle aussi de
condition au bord). À l'instant $0$ nous avons $y(0) = \Gamma + b/a$, soit
$\Gamma = y_0 - b/a$. Ainsi la solution de l'équation différentielle non
homogène est\nbsp{}:

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

#+BEGIN_remarque
Cette solution n'est valable que si $a\neq0$. Dans le cas contraire, on obtient
la solution beaucoup plus simplement. Nous aurions en effet $\dot y(t) = b$
et, en intégrant par rapport à $t$, $y(t) = \gamma + bt$ où $\gamma$ est une
constante arbitraire que l'on fixe à l'aide d'une condition initiale. La
démarche en quatre étapes n'est pas conseillée dans ce cas, même si elle
reste praticable (exercice [[ex-a-nul][ci-dessous]]).
#+END_remarque

#+BEGIN_remarque
La solution reçoit la même interprétation que dans le cas homogène. La
dynamique est stable si et seulement si le paramètre $a$ est positif. Dans ce
cas, pour toute condition initiale $y(0)$, on a $\lim_{t\to\infty}y(t)=b/a$.
Si le paramètre $a$ est négatif, la variable $y(t)$ diverge vers $+\infty$ ou
$-\infty$ selon le signe de $y(0)-b/a$, dès lors que $y(0)\neq b/a$. Si la
condition initiale est égale à l'état stationnaire, alors $y(t)=b/a$ pour tout
$t$, indépendamment du signe de $a$. Il est possible d'écrire la solution sous
la forme suivante\nbsp{}:

\[
y(t) = e^{-at}\,y_0 + \left(1-e^{-at}\right)\frac{b}{a}
\]

Le niveau de $y$ à l'instant $t$ est une combinaison convexe de la condition
initiale et de l'état stationnaire. Si $a>0$, l'influence de la condition
initiale tend à disparaître alors que la pondération de l'état stationnaire
(la cible) tend vers un.
#+END_remarque

#+NAME: ex-a-nul
#+BEGIN_exercice
Montrez que, même dans le cas où $a=0$, il est possible de résoudre l'équation
différentielle non homogène $\dot y(t) + a\,y(t) = b$ en suivant l'approche en
quatre étapes décrite plus haut. Indice\nbsp{}: cherchez une solution /non
constante/ de l'équation complète.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
Avec $a=0$ l'équation s'écrit $\dot y(t) = b$. /Étape 1/\nbsp{}: l'équation sans
second membre est $\dot y(t)=0$, dont la solution générale est la fonction
constante $y_1(t)=\Gamma$. /Étape 2/\nbsp{}: une fonction constante ne peut pas
être solution de l'équation complète puisque sa dérivée est nulle et non
égale à \(b\)\nbsp{}; on cherche donc une solution particulière linéaire en $t$,
$y_2(t)=ct$, et l'équation impose $c=b$. /Étape 3/\nbsp{}: la solution générale
est $y(t) = \Gamma + bt$. /Étape 4/\nbsp{}: la condition initiale donne
$\Gamma=y_0$, d'où $y(t) = y_0 + bt$. On retrouve le résultat de l'intégration
directe. Notons que l'on ne peut pas obtenir cette solution comme limite de
$\left(y_0 - b/a\right)e^{-at} + b/a$ lorsque $a$ tend vers zéro terme à
terme, chacun des deux termes divergeant\nbsp{}; mais la forme
$e^{-at}y_0 + \frac{1-e^{-at}}{a}\,b$ a bien pour limite $y_0+bt$, puisque
$\frac{1-e^{-at}}{a}\to t$.
#+END_details

** Dynamique d'un prix de marché
:PROPERTIES:
:CUSTOM_ID: prix
:END:

Supposons que les fonctions de demande et d'offre d'un bien soient données
par\nbsp{}:

\begin{cases}
Q_d = \alpha - \beta P & (\alpha,\beta>0)\\
Q_s = \gamma + \delta P & (\gamma,\delta>0)
\end{cases}

où $P$ est le prix du bien. Les paramètres $\beta$ et $\delta$ mesurent la
sensibilité de la demande et de l'offre à une variation du prix (il ne s'agit
pas d'élasticités). En égalisant les quantités offerte et demandée, on obtient
le prix qui apure le marché\nbsp{}:

\[
P^{\star} = \frac{\alpha-\gamma}{\beta+\delta}
\]

que l'on supposera strictement positif, ce qui exige \(\alpha>\gamma\)\nbsp{}: à
prix nul, la demande excède l'offre. A priori, le prix effectif $P$ est
différent de $P^{\star}$. Nous allons montrer que si le prix augmente lorsque
la demande est supérieure à l'offre, et diminue lorsque l'offre est supérieure
à la demande, alors le prix effectif converge vers le prix qui apure le
marché. Nous formalisons cette hypothèse en supposant que\nbsp{}:

\[
\dot P(t) = j\,\bigl(Q_d(t) - Q_s(t)\bigr)
\]

avec $j>0$. Le prix est invariant si et seulement si l'offre égalise la
demande. En substituant les fonctions d'offre et de demande, on obtient
l'équation différentielle suivante pour le prix du bien\nbsp{}:

\[
\dot P(t) + j(\beta+\delta)\,P(t) = j(\alpha-\gamma)
\]

En appliquant la formule de la section précédente, on obtient
directement\nbsp{}:

\[
P(t) = \bigl[P(0) - P^{\star}\bigr]e^{-\kappa t} + P^{\star}
\]

avec $\kappa \equiv j(\beta+\delta)$. Puisque $\kappa>0$, le prix converge à
long terme vers $P^{\star}$. La convergence est monotone croissante si $P(0)$
est inférieur à $P^{\star}$, c'est-à-dire s'il y a initialement excès de
demande, et monotone décroissante si $P(0)$ est supérieur à $P^{\star}$,
c'est-à-dire s'il y a initialement excès d'offre. $P^{\star}$ est une
solution particulière de l'équation différentielle\nbsp{}; elle représente le
niveau d'équilibre intertemporel de la variable d'intérêt. Le terme
$[P(0)-P^{\star}]e^{-\kappa t}$ rend compte des déviations à ce niveau
d'équilibre.

#+begin_src python :session edo1 :exports none :results none
  al, be, ga, de, j = 10.0, 1.0, 2.0, 1.0, 0.5
  Pstar = (al-ga)/(be+de)
  kappa = j*(be+de)
  fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
  P = np.linspace(0, 8, 100)
  ax1.plot(al-be*P, P, 'b', label=r'$Q_d = \alpha-\beta P$')
  ax1.plot(ga+de*P, P, 'g', label=r'$Q_s = \gamma+\delta P$')
  ax1.axhline(y=Pstar, color='r', linewidth=1, linestyle='--')
  ax1.annotate(r'$P^\star$', xy=(0.3, Pstar+0.2))
  ax1.set_xlabel(r'$Q$')
  ax1.set_ylabel(r'$P$')
  ax1.legend()
  t = np.linspace(0, 6, 300)
  for P0, couleur, nom in ((1.0, 'b', "excès de demande initial"), (7.0, 'g', "excès d'offre initial")):
      ax2.plot(t, (P0-Pstar)*np.exp(-kappa*t)+Pstar, couleur, label=nom)
  ax2.axhline(y=Pstar, color='r', linewidth=1, linestyle='--')
  ax2.annotate(r'$P^\star$', xy=(5.5, Pstar+0.2))
  ax2.set_xlabel(r'$t$')
  ax2.set_ylabel(r'$P(t)$')
  ax2.legend()
  fig.tight_layout()
  fig.savefig("edo-prix.svg", transparent=True)
#+end_src

#+CAPTION: *À gauche, offre et demande et prix d'équilibre ($\alpha=10$, $\beta=\delta=1$, $\gamma=2$). À droite, ajustement du prix depuis une situation d'excès de demande et depuis une situation d'excès d'offre ($j=0{,}5$).*
#+LABEL: fig:prix
[[file:edo-prix.svg]]

#+BEGIN_exercice
Résoudre les équations différentielles suivantes\nbsp{}:
- (i) $\dot y(t) + 4y(t) = 12$ avec $y(0)=2$.
- (ii) $\dot y(t) - 2y(t) = 0$ avec $y(0)=9$.
- (iii) $\dot y(t) + 10y(t) = 15$ avec $y(0)=0$.
- (iv) $2\dot y(t) + 4y(t) = 6$ avec $y(0)=3/2$.
- (v) $\dot y(t) + y(t) = 4$ avec $y(0)=0$.
- (vi) $\dot y(t) = 23$ avec $y(0)=1$.
- (vii) $3\dot y(t) + 6y(t) = 5$ avec $y(0)=0$.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
On applique la formule $y(t) = (y_0 - b/a)e^{-at} + b/a$, après avoir
éventuellement divisé l'équation par le coefficient de $\dot y$.
- (i) $a=4$, $b=12$, état stationnaire \(3\)\nbsp{}: $y(t) = 3 - e^{-4t}$.
- (ii) Équation homogène avec \(a=-2\)\nbsp{}: $y(t) = 9e^{2t}$, la dynamique est
  instable.
- (iii) $a=10$, $b=15$, état stationnaire \(3/2\)\nbsp{}:
  $y(t) = \frac32\left(1-e^{-10t}\right)$.
- (iv) En divisant par $2$, $\dot y + 2y = 3$, dont l'état stationnaire est
  $3/2$. La condition initiale est égale à l'état stationnaire\nbsp{}: $y(t)=3/2$
  pour tout $t$.
- (v) $a=1$, \(b=4\)\nbsp{}: $y(t) = 4\left(1-e^{-t}\right)$.
- (vi) \(a=0\)\nbsp{}: on intègre directement, $y(t) = 1 + 23t$.
- (vii) En divisant par $3$, $\dot y + 2y = 5/3$, dont l'état stationnaire est
  \(5/6\)\nbsp{}: $y(t) = \frac56\left(1-e^{-2t}\right)$.
#+END_details

* Équations linéaires à coefficients variables
:PROPERTIES:
:CUSTOM_ID: variables
:END:

** Le cas homogène
:PROPERTIES:
:CUSTOM_ID: homogene-variable
:END:

On s'intéresse à un problème de la forme\nbsp{}:

\[
\dot y(t) + a(t)\,y(t) = 0
\]

où $a(t)$ est une fonction réelle continue. Nous pouvons résoudre cette
équation en suivant la même démarche que dans le cas des [[#homogene][coefficients
constants]]. Le taux de croissance de $y$ à l'instant $t$ est $\dot y(t)/y(t) =
-a(t)$, ou encore, en utilisant les propriétés de la fonction logarithme\nbsp{}:

\[
\frac{\mathrm d}{\mathrm dt}\log y(t) = -a(t)
\]

En intégrant les deux membres entre $0$ et $t$, la condition initiale étant
supposée connue, il vient $\log y(t) - \log y(0) = -\int_0^t a(\tau)\,\mathrm d\tau$,
soit de façon équivalente\nbsp{}:

\[
y(t) = y(0)\,e^{-\int_0^t a(\tau)\,\mathrm d\tau}
\]

On vérifie que si $a(t)=a$ pour tout $t$, alors $\int_0^t a(\tau)\,\mathrm d\tau = at$
et l'on retrouve la solution obtenue dans le cas des coefficients constants.
Les conditions de stabilité sont ici moins triviales, puisque tout dépend de
la forme de la fonction $a(t)$. En notant $A(t) = \int_0^t a(\tau)\,\mathrm d\tau$,
on a $|y(t)| = |y(0)|e^{-A(t)}$ et l'on distingue trois cas\nbsp{}:
- la variable $y$ converge vers l'état stationnaire nul si et seulement si
  $A(t)$ tend vers \(+\infty\)\nbsp{};
- elle converge vers une limite finie non nulle si et seulement si $A(t)$
  admet une limite finie\nbsp{};
- elle reste bornée, sans nécessairement converger, si et seulement si $A(t)$
  est minorée.
La dynamique diverge dans le cas restant, c'est-à-dire si $A(t)$ n'est pas
minorée. Une fonction $a$ minorée par une constante strictement positive
garantit la convergence vers zéro, comme dans le cas constant. Mais
$a(t) = (1+t)^{-2}$, pourtant strictement positive, donne $A(t) = 1-(1+t)^{-1}$
et $y(t)$ converge vers $y(0)e^{-1}$ sans jamais rejoindre l'état
stationnaire\nbsp{}; et $a(t) = \cos t$ donne $A(t) = \sin t$, de sorte que
$y(t)$ oscille indéfiniment entre $y(0)e^{-1}$ et $y(0)e$. Le signe de $a$ à
un instant donné ne renseigne que sur le sens de variation de $y$ à cet
instant.

** Équations différentielles exactes
:PROPERTIES:
:CUSTOM_ID: exactes
:END:

Nous quittons un instant le monde des équations différentielles linéaires
afin d'introduire un outil qui nous permettra de construire la solution d'une
équation différentielle linéaire non homogène à coefficients variables.

Soit une fonction de deux variables $F(y,t)$, continûment différentiable. Sa
différentielle totale est\nbsp{}:

\[
\mathrm dF(y,t) = \frac{\partial F}{\partial y}\,\mathrm dy + \frac{\partial F}{\partial t}\,\mathrm dt
\]

#+BEGIN_definition
L'équation
\[
\frac{\partial F}{\partial y}\,\mathrm dy + \frac{\partial F}{\partial t}\,\mathrm dt = 0
\]
est une /équation différentielle exacte/, car le membre de gauche est
exactement la différentielle de $F(y,t)$.
#+END_definition

#+BEGIN_exemple
Soit $F(y,t) = y^2t + k$, où $k$ est une constante. Nous avons
$\mathrm dF = 2yt\,\mathrm dy + y^2\,\mathrm dt$. Ainsi $2yt\,\mathrm dy + y^2\,\mathrm dt = 0$, ou
de façon équivalente $\dot y + \frac{y^2}{2yt} = 0$, est une équation
différentielle exacte.
#+END_exemple

De façon générale, l'équation différentielle\nbsp{}:

\[
M\,\mathrm dy + N\,\mathrm dt = 0
\]

est exacte si et seulement s'il existe une fonction $F(y,t)$ telle que
$M = \partial F/\partial y$ et $N = \partial F/\partial t$. Il nous manque un
test pour savoir si une équation différentielle est exacte, c'est-à-dire
pour savoir s'il existe une fonction $F$ à l'origine des fonctions $M$ et
$N$. Nous savons qu'une matrice hessienne est symétrique (théorème de
Young), c'est-à-dire que $\frac{\partial^2F}{\partial t\partial y} =
\frac{\partial^2F}{\partial y\partial t}$. Ainsi l'équation différentielle est
exacte si et seulement si\nbsp{}:

\[
\frac{\partial M}{\partial t} = \frac{\partial N}{\partial y}
\]

Cette condition nous donne un test pour évaluer si une équation
différentielle est exacte[fn:: La condition est nécessaire. Elle est
suffisante sur un domaine simplement connexe, ce que nous supposerons
toujours.].

#+BEGIN_exemple
Si nous appliquons ce test à l'exemple précédent, où $M=2yt$ et $N=y^2$, il
vient $\partial M/\partial t = 2y$ et $\partial N/\partial y = 2y$. Le test
conclut qu'il s'agit bien d'une équation différentielle exacte.
#+END_exemple

Une équation différentielle exacte, par définition, nous dit que
$\mathrm dF(y,t) = 0$. Ainsi sa solution générale doit être de la forme\nbsp{}:

\[
F(y,t) = c
\]

où $c$ est une constante réelle. Résoudre une équation différentielle exacte,
c'est exhiber une primitive $F(y,t)$ et l'égaliser à une constante.

*Méthode de résolution.* Puisque $M = \partial F/\partial y$, la fonction $F$
doit contenir une intégrale de $M$ par rapport à la variable $y$. Nous devrions
donc avoir\nbsp{}:

\[
F(y,t) = \int M\,\mathrm dy + \psi(t)
\]

La dérivée partielle $M$ est intégrée seulement par rapport à $y$, en
traitant $t$ comme une constante. Puisqu'en différentiant $F(y,t)$
partiellement par rapport à $y$ tout terme additif ne dépendant pas de $y$
disparaît, on doit prendre soin de réintroduire ces termes dans le processus
d'intégration. C'est exactement le rôle du terme $\psi(t)$. Il est
relativement aisé d'évaluer \(\int M\,\mathrm dy\)\nbsp{}; déterminer la fonction
$\psi(t)$ est souvent moins évident.

#+BEGIN_exemple
On veut résoudre l'équation différentielle $\dot y + \frac{y}{2t} = 0$. En
multipliant les deux membres par $2yt\,\mathrm dt$, il vient
$2yt\,\mathrm dy + y^2\,\mathrm dt = 0$. Nous avons donc $M=2yt$ et $N=y^2$, et nous
résolvons cette équation en quatre étapes.
- (i) Posons $F(y,t) = \int 2yt\,\mathrm dy + \psi(t) = y^2t + \psi(t)$, où
  $\psi(t)$ reste à déterminer et où nous avons redéfini $\psi(t)$ afin
  d'inclure la constante d'intégration.
- (ii) La dérivée partielle par rapport à $t$ est $\partial F/\partial t =
  y^2 + \psi'(t)$. En comparant avec $N=y^2$, on en déduit que la fonction
  $\psi(t)$ doit être telle que $\psi'(t)=0$ pour tout $t$.
- (iii) Nous savons donc que $\psi(t) = \kappa\in\mathbb R$ pour tout $t$.
- (iv) Finalement, $F(y,t) = y^2t + \kappa$ et la solution de l'équation
  différentielle exacte doit être de la forme $y^2t + \kappa = c$. Puisque
  les constantes $\kappa$ et $c$ ne sont pas individuellement identifiables,
  nous avons encore $y^2t = \tilde c$, avec $\tilde c\geq0$, d'où finalement,
  pour \(t>0\)\nbsp{}:
  \[
  y(t) = \bar c\,t^{-\frac12}
  \]
  où $\bar c = \pm\sqrt{\tilde c}$ est une constante réelle, du signe de $y$,
  qui pourra être déterminée à l'aide d'une condition initiale.
#+END_exemple

#+BEGIN_exemple
Soit l'équation différentielle\nbsp{}:
\[
\dot y + \frac{y+3t^2}{t+2y} = 0
\]
De façon équivalente, nous avons $(t+2y)\,\mathrm dy + (y+3t^2)\,\mathrm dt = 0$.
S'agit-il d'une équation différentielle exacte\nbsp{}? Nous avons $M = t+2y$ et
$N = y+3t^2$, et nous vérifions que $\partial M/\partial t = 1 = \partial
N/\partial y$. Nous sommes donc bien en présence d'une équation
différentielle exacte, que nous résolvons en quatre étapes.
- (i) Posons $F(y,t) = \int(t+2y)\,\mathrm dy + \psi(t) = yt + y^2 + \psi(t)$,
  où $\psi(t)$ reste à déterminer.
- (ii) La dérivée partielle par rapport à $t$ est $\partial F/\partial t =
  y + \psi'(t)$. En comparant avec $N = y+3t^2$, on en déduit que
  $\psi'(t) = 3t^2$ pour tout $t$.
- (iii) Nous savons donc que $\psi(t) = t^3 + \kappa$ avec $\kappa\in\mathbb R$.
- (iv) Finalement, $F(y,t) = yt + y^2 + t^3 + \kappa$ et la solution de
  l'équation différentielle est de la forme\nbsp{}:
  \[
  yt + y^2 + t^3 = c
  \]
  où $c$ est une constante réelle. On vérifie qu'il s'agit bien de la solution
  en différentiant cette relation par rapport à $t$. Elle définit
  implicitement $y$ comme une fonction de \(t\)\nbsp{}; ici on peut même
  l'expliciter, puisqu'il s'agit d'une équation du second degré en $y$.
#+END_exemple

La procédure décrite dans ces deux exemples peut être appliquée à toute
équation différentielle exacte. Dans certains cas, la procédure peut aussi
être appliquée à une équation différentielle non exacte, si l'on peut trouver
une équation différentielle exacte équivalente.

#+BEGIN_exemple
Soit l'équation différentielle $2t\,\mathrm dy + y\,\mathrm dt = 0$. On vérifie
facilement qu'elle n'est pas exacte\nbsp{}: $\partial M/\partial t = 2$ et
$\partial N/\partial y = 1$. Néanmoins, en multipliant chaque terme par $y$,
on revient à l'exemple précédent et l'on obtient donc une équation
différentielle exacte. On dit que $y$ est un /facteur d'intégration/ de
l'équation.
#+END_exemple

** Le cas non homogène
:PROPERTIES:
:CUSTOM_ID: non-homogene-variable
:END:

On s'intéresse à un problème de la forme\nbsp{}:

\[
\dot y(t) + a(t)\,y(t) = b(t)
\]

où $a(t)$ et $b(t)$ sont des fonctions réelles continues. Ce problème peut
s'écrire de façon équivalente\nbsp{}:

\[
\mathrm dy + \bigl(a(t)y - b(t)\bigr)\mathrm dt = 0
\]

Nous avons $M=1$ et \(N = ay-b\)\nbsp{}; on voit immédiatement qu'il ne s'agit
pas d'une équation différentielle exacte, sauf si $a$ est identiquement
nulle. Néanmoins il est possible d'exhiber un facteur d'intégration, que nous
noterons $\mathcal I$, de façon à obtenir une équation différentielle exacte.
Le facteur d'intégration est tel que\nbsp{}:

\[
\mathcal I\,\mathrm dy + \mathcal I\bigl(a(t)y - b(t)\bigr)\mathrm dt = 0
\]

est une équation différentielle exacte. Pour cela, il faut et il suffit que
la condition $\partial M/\partial t = \partial N/\partial y$, avec
$M = \mathcal I$ et $N = \mathcal I(ay-b)$, soit satisfaite. Le facteur
d'intégration est donc tel que $\dot{\mathcal I} = \mathcal I a$, soit de façon
équivalente\nbsp{}:

\[
\frac{\dot{\mathcal I}(t)}{\mathcal I(t)} = a(t)
\]

Le taux de croissance du facteur d'intégration doit être égal à $a(t)$. Le
facteur d'intégration est donc défini par une équation différentielle que
nous savons résoudre\nbsp{}:

\[
\mathcal I(t) = A\,e^{\int_0^t a(\tau)\,\mathrm d\tau}
\]

pour toute valeur réelle non nulle de $A$. Sans perte de généralité nous
poserons $A=1$.

#+BEGIN_property
La solution générale de l'équation différentielle linéaire du premier ordre à
coefficients variables $\dot y(t) + a(t)\,y(t) = b(t)$ est\nbsp{}:
\[
y(t) = e^{-\int_0^t a(\tau)\,\mathrm d\tau}\left(c + \int_0^t b(s)\,e^{\int_0^s a(\tau)\,\mathrm d\tau}\,\mathrm ds\right)
\]
où la constante $c$ est déterminée par la condition initiale\nbsp{}: si $y(0)=y_0$
est connu, alors $c = y_0$.
#+END_property

#+BEGIN_proof
L'équation différentielle transformée\nbsp{}:
\[
e^{\int_0^t a(\tau)\,\mathrm d\tau}\,\mathrm dy + e^{\int_0^t a(\tau)\,\mathrm d\tau}\bigl(a(t)y - b(t)\bigr)\mathrm dt = 0
\]
est exacte par construction. Nous la résolvons en suivant les quatre étapes
décrites plus haut.
- (i) Posons $F(y,t) = \int e^{\int_0^t a(\tau)\,\mathrm d\tau}\,\mathrm dy + \psi(t)
  = y\,e^{\int_0^t a(\tau)\,\mathrm d\tau} + \psi(t)$, où $\psi(t)$ reste à
  déterminer.
- (ii) La dérivée partielle par rapport à $t$ est\nbsp{}:
  \[
  \frac{\partial F}{\partial t} = y\,a(t)\,e^{\int_0^t a(\tau)\,\mathrm d\tau} + \psi'(t)
  \]
  En comparant avec $N$ nous obtenons une restriction sur la fonction
  \(\psi\)\nbsp{}: $\psi'(t) = -b(t)\,e^{\int_0^t a(\tau)\,\mathrm d\tau}$.
- (iii) La fonction $\psi$ est donc, l'instant initial étant en zéro\nbsp{}:
  \[
  \psi(t) = \int_0^t\psi'(s)\,\mathrm ds = -\int_0^t b(s)\,e^{\int_0^s a(\tau)\,\mathrm d\tau}\,\mathrm ds
  \]
  Nous ne pouvons aller plus loin ici, car les fonctions $a$ et $b$ ne sont
  pas spécifiées.
- (iv) Finalement, en substituant l'expression de $\psi(t)$ dans la fonction
  $F(y,t)$ postulée, la solution vérifie\nbsp{}:
  \[
  y(t)\,e^{\int_0^t a(\tau)\,\mathrm d\tau} - \int_0^t b(s)\,e^{\int_0^s a(\tau)\,\mathrm d\tau}\,\mathrm ds = c
  \]
  d'où l'expression annoncée. En $t=0$ les deux intégrales sont nulles et
  l'on obtient $y(0)=c$.
#+END_proof

#+BEGIN_exemple
Soit l'équation différentielle $\dot y(t) + 2t\,y(t) = t$, la condition
initiale $y(0)$ étant connue. Nous avons $a(t)=2t$ et $b(t)=t$, d'où
$\int_0^t a(\tau)\,\mathrm d\tau = t^2$. En appliquant le résultat précédent, il
vient\nbsp{}:
\[
y(t) = e^{-t^2}\left(y(0) + \int_0^t s\,e^{s^2}\,\mathrm ds\right)
\]
En notant que $\frac{\mathrm d}{\mathrm ds}e^{s^2} = 2s\,e^{s^2}$, l'intégrale vaut
$\frac12\left(e^{t^2}-1\right)$ et\nbsp{}:
\[
y(t) = e^{-t^2}\left(y(0) - \frac12\right) + \frac12
\]
Notons que $1/2$ est une solution particulière constante de l'équation
différentielle\nbsp{}: il s'agit d'un état stationnaire, stable dans cet exemple
puisque le terme exponentiel tend vers zéro lorsque $t$ tend vers l'infini.
On aurait pu l'obtenir directement en cherchant, comme dans le cas des
coefficients constants, une solution particulière constante, puis en lui
ajoutant la solution générale $\Gamma e^{-t^2}$ de l'équation homogène.
#+END_exemple

#+BEGIN_exercice
Résoudre les équations différentielles suivantes\nbsp{}:
- (a) $\dot y + 5y = 15$.
- (b) $\dot y + 2ty = 0$.
- (c) $\dot y + 2ty = t$ avec $y(0)=3/2$.
- (d) $\dot y + t^2y = 5t^2$ avec $y(0)=6$.
- (e) $2\dot y + 12y + 2e^t = 0$ avec $y(0)=6/7$.
- (f) $\dot y + y = t$.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
Lorsque la condition initiale n'est pas donnée, on se contente de la solution
générale, indexée par une constante $\Gamma$.
- (a) Coefficients constants, état stationnaire \(3\)\nbsp{}: $y(t) = 3 + \Gamma e^{-5t}$,
  avec $\Gamma = y(0)-3$.
- (b) Équation homogène avec $a(t)=2t$, donc \(\int_0^t a = t^2\)\nbsp{}:
  $y(t) = \Gamma e^{-t^2}$, avec $\Gamma = y(0)$.
- (c) C'est l'exemple traité plus haut\nbsp{}: $y(t) = \frac12 + \left(\frac32-\frac12\right)e^{-t^2}
  = \frac12 + e^{-t^2}$.
- (d) $a(t)=t^2$, $b(t)=5t^2$, donc $\int_0^t a = t^3/3$ et le facteur
  d'intégration est $e^{t^3/3}$. On observe que $5$ est une solution
  particulière constante, d'où $y(t) = 5 + \Gamma e^{-t^3/3}$ et, avec la
  condition initiale, $y(t) = 5 + e^{-t^3/3}$. La formule générale donne le
  même résultat, l'intégrale $\int_0^t 5s^2e^{s^3/3}\,\mathrm ds = 5\left(e^{t^3/3}-1\right)$
  se calculant sans difficulté.
- (e) En divisant par $2$, $\dot y + 6y = -e^t$. L'équation homogène a pour
  solution générale $\Gamma e^{-6t}$. Le second membre n'étant pas constant,
  on cherche une solution particulière de la même forme que lui, \(y_2(t) =
  Ae^t\)\nbsp{}: en substituant, $A + 6A = -1$, soit $A = -1/7$. La solution
  générale est $y(t) = \Gamma e^{-6t} - \frac17e^t$ et la condition initiale
  donne $\Gamma - \frac17 = \frac67$, soit \(\Gamma=1\)\nbsp{}:
  $y(t) = e^{-6t} - \frac17 e^t$. La solution diverge vers $-\infty$, bien que
  le coefficient $a=6$ soit positif\nbsp{}: c'est le second membre qui diverge.
- (f) On cherche une solution particulière affine, \(y_2(t) = ct + d\)\nbsp{}: en
  substituant, $c + ct + d = t$ impose $c=1$ et $d=-1$. La solution générale
  est $y(t) = t - 1 + \Gamma e^{-t}$, avec $\Gamma = y(0)+1$. À long terme $y$
  suit la droite $t-1$ avec un retard unitaire sur $t$.
#+END_details

* Équations non linéaires
:PROPERTIES:
:CUSTOM_ID: non-lineaires
:END:

Dans cette section nous nous intéressons aux équations différentielles qui
peuvent s'écrire sous la forme\nbsp{}:

\[
f(y,t)\,\mathrm dy + g(y,t)\,\mathrm dt = 0
\]

ou, de façon équivalente, $\dot y(t) = h\bigl(y(t),t\bigr)$ avec $h(y,t) =
-g(y,t)/f(y,t)$. On peut éventuellement reconnaître une équation
différentielle exacte, si la condition $\partial f/\partial t = \partial
g/\partial y$ est vérifiée\nbsp{}; nous savons déjà [[#exactes][résoudre]] ce type
d'équation.

** Problèmes séparables
:PROPERTIES:
:CUSTOM_ID: separables
:END:

Nous envisageons ici le cas où $f$ ne dépend pas de $t$ et $g$ ne dépend pas de
$y$. Cette classe de problèmes est très simple à résoudre\nbsp{}: chaque membre
s'intègre séparément.

#+BEGIN_exemple
Soit l'équation différentielle $3y^2\,\mathrm dy = t\,\mathrm dt$. En intégrant les
deux membres, $\int 3y^2\,\mathrm dy = \int t\,\mathrm dt$, on obtient directement
$y^3 + \gamma_1 = \frac12t^2 + \gamma_2$, ou encore $y^3 = \frac12t^2 + \gamma$,
puisque les deux constantes d'intégration ne sont pas identifiables
séparément. Nous avons donc\nbsp{}:
\[
y(t) = \left(\frac12t^2 + \gamma\right)^{\frac13}
\]
On peut alors fixer $\gamma$ à l'aide de la condition initiale, en évaluant la
dernière équation en $t=0$. Il vient $\gamma = y(0)^3$ et\nbsp{}:
\[
y(t) = \left(\frac12t^2 + y(0)^3\right)^{\frac13}
\]
#+END_exemple

#+BEGIN_exemple
Soit l'équation différentielle $2t\,\mathrm dy + y\,\mathrm dt = 0$. A priori, elle
n'appartient pas à la classe discutée ici, puisque $\mathrm dy$ est associé à une
fonction de $t$ et non de $y$. Mais en divisant les deux membres par $2yt$,
nous obtenons\nbsp{}:
\[
\frac{\mathrm dy}{y} + \frac{\mathrm dt}{2t} = 0
\]
Notons au passage que cette transformation rend l'équation différentielle
exacte, en plus de la rendre séparable. En intégrant, il vient
$\log y + \frac12\log t = \gamma$, soit $y\,t^{\frac12} = e^{\gamma}$ et
finalement\nbsp{}:
\[
y(t) = \Gamma\,t^{-\frac12}
\]
Cette solution n'est pas définie en zéro, ce que nous savions depuis la
transformation du problème initial, qui suppose implicitement $y\neq0$ et
\(t\neq0\)\nbsp{}: pour fixer la constante $\Gamma$, il faut choisir un instant
initial strictement positif. On retrouve, par une troisième voie, la solution
de l'[[#exactes][exemple]] traité avec un facteur d'intégration.
#+END_exemple

** Réduction à une dynamique linéaire : l'équation de Bernoulli
:PROPERTIES:
:CUSTOM_ID: bernoulli
:END:

#+BEGIN_definition
L'équation différentielle de la forme\nbsp{}:
\[
\dot y(t) + R(t)\,y(t) = Q(t)\,y(t)^m
\]
avec $m\notin\{0,1\}$, est une /équation de Bernoulli/.
#+END_definition

Pour $m=0$ on retrouve l'équation linéaire non homogène, et pour $m=1$
l'équation linéaire homogène de coefficient $R-Q$. Dans les autres cas
l'équation est non linéaire, mais elle peut toujours se réduire à une équation
différentielle linéaire par changement de variable. Divisons les deux membres
par \(y^m\)\nbsp{}:

\[
y(t)^{-m}\dot y(t) + R(t)\,y(t)^{1-m} = Q(t)
\]

et posons $z = y^{1-m}$, de sorte que $\dot z = (1-m)y^{-m}\dot y$. L'équation
s'écrit alors $\frac{1}{1-m}\dot z(t) + R(t)\,z(t) = Q(t)$, ou encore\nbsp{}:

\[
\dot z(t) + (1-m)\bigl[R(t)\,z(t) - Q(t)\bigr] = 0
\]

Il s'agit d'une équation différentielle linéaire à coefficients variables.
Nous pouvons donc la résoudre, puis exprimer $y$ en fonction de $z$ afin
d'obtenir la solution de l'équation de départ.

#+NAME: ex-bernoulli
#+BEGIN_exemple
Soit l'équation différentielle\nbsp{}:
\[
\dot y(t) + t\,y(t) = 3t\,y(t)^2
\]
En divisant par $y^2$, il vient $\dot y\,y^{-2} + t\,y^{-1} - 3t = 0$. Posons
$z = y^{-1}$, de sorte que \(\dot z = -y^{-2}\dot y\)\nbsp{}; la dynamique de $z$
s'écrit $-\dot z + tz - 3t = 0$, ou de façon équivalente\nbsp{}:
\[
\dot z - t\,z + 3t = 0
\]
Il s'agit d'une équation linéaire à coefficients variables, avec $a(t)=-t$
et $b(t)=-3t$. En appliquant la [[#non-homogene-variable][formule générale]], il vient
directement\nbsp{}:
\[
z(t) = e^{\frac{t^2}{2}}\left(A - 3\int_0^t s\,e^{-\frac{s^2}{2}}\,\mathrm ds\right)
= e^{\frac{t^2}{2}}\left(A-3\right) + 3
\]
Il nous reste à renverser la transformation, puisque la variable d'intérêt
est $y$ et non $z$. Nous avons $y = z^{-1}$ et, en utilisant la condition
initiale pour déterminer la constante, \(A = 1/y(0)\)\nbsp{}:
\[
y(t) = \frac{1}{e^{\frac{t^2}{2}}\left(\frac{1}{y(0)}-3\right) + 3}
\]
#+END_exemple

#+BEGIN_exercice
Cherchez la solution de l'équation différentielle suivante\nbsp{}:
\[
\dot y(t) + \frac{1}{t}\,y(t) = y(t)^3
\]
avec $t_0>0$ l'instant initial, la condition initiale $y(t_0)$ étant connue.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Corrigé
#+END_summary
C'est une équation de Bernoulli avec $m=3$, $R(t)=1/t$ et $Q(t)=1$. On divise
par $y^3$ et l'on pose $z = y^{-2}$, de sorte que $\dot z = -2y^{-3}\dot y$. Il
vient $-\frac12\dot z + \frac{z}{t} = 1$, soit\nbsp{}:
\[
\dot z - \frac{2}{t}\,z = -2
\]
Le facteur d'intégration est $e^{-\int 2/t} = t^{-2}$, et l'équation s'écrit
$\frac{\mathrm d}{\mathrm dt}\left(t^{-2}z\right) = -2t^{-2}$, d'où $t^{-2}z = \frac{2}{t} + C$
et $z(t) = 2t + Ct^2$. On revient à \(y\)\nbsp{}:
\[
y(t) = \pm\left(2t + Ct^2\right)^{-\frac12}
\]
le signe étant celui de $y(t_0)$, et la constante étant fixée par la
condition initiale, $C = \left(y(t_0)^{-2} - 2t_0\right)/t_0^2$. La solution
n'est définie que tant que $2t + Ct^2$ reste strictement positif\nbsp{}: si
$C < 0$, c'est-à-dire si $y(t_0)^2 > 1/(2t_0)$, elle explose en un temps fini
$t^{\star} = -2/C$. On vérifie que $y = 0$ est un état stationnaire, et que
pour $C\geq0$ la solution converge vers zéro comme $t^{-1}$ ou $t^{-1/2}$.
#+END_details

** Approche qualitative
:PROPERTIES:
:CUSTOM_ID: qualitative
:END:

On ne peut pas toujours résoudre analytiquement une équation différentielle
non linéaire. Dans de nombreux cas nous devons nous contenter de solutions
numériques, pour une approche quantitative (voir la section
[[#numerique][suivante]]), ou graphiques, si une approche qualitative peut suffire. À
l'aide d'une approche graphique, on peut par exemple s'interroger sur la
stabilité de la dynamique.

On s'intéresse ici à une dynamique de la forme\nbsp{}:

\[
\dot y = f(y)
\]

où $f$ est une fonction continue ne dépendant pas du temps. On parle alors
d'une équation différentielle /autonome/. On peut représenter graphiquement
la dynamique dans le plan $(y,\dot y)$. La figure [[fig:phase][ci-dessous]] donne deux
exemples\nbsp{}: dans le graphique de gauche la fonction $f$ est monotone
décroissante, dans celui de droite elle est monotone croissante. La lecture
de ces graphiques repose sur trois observations\nbsp{}:
1. lorsque $f(y)>0$, la variation $\dot y$ est positive et donc $y$
   augmente\nbsp{};
2. lorsque $f(y) < 0$, la variation $\dot y$ est négative et donc $y$
   diminue\nbsp{};
3. lorsque $f(y)=0$, la variation $\dot y$ est nulle et donc $y$ ne bouge pas.
Ces trois observations expliquent comment nous avons orienté les flèches sur
l'axe des abscisses. Ces flèches décrivent le sens de variation de la variable
$y$. Par exemple, dans le graphique de gauche, $f(y)$ est positif lorsque
$y < y_a^{\star}$ et négatif lorsque $y>y_a^{\star}$. Ainsi $y$ augmente lorsque
son niveau est faible relativement à $y_a^{\star}$, et les flèches sont
dirigées vers la droite\nbsp{}; $y$ diminue lorsque son niveau est élevé, et
les flèches sont dirigées vers la gauche.

#+begin_src python :session edo1 :exports none :results none
  def fleches(ax, segments, y=0.0):
      """Flèches sur l'axe des abscisses, de x0 vers x1 pour chaque segment."""
      for x0, x1 in segments:
          ax.annotate('', xy=(x1, y), xytext=(x0, y),
                      arrowprops=dict(arrowstyle='-|>', color='g', linewidth=1.5))

  fa = lambda y: 2*(1-y) + 0.3*(1-y)**3
  grille = np.linspace(0, 2, 200)
  fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
  for ax, signe, nom in ((ax1, 1, 'a'), (ax2, -1, 'b')):
      ax.plot(grille, signe*fa(grille), 'b')
      ax.axhline(y=0, color='k', linewidth=0.8)
      ax.plot([1], [0], 'ko', markersize=5)
      ax.annotate(r'$y_%s^\star$' % nom, xy=(1, 0), xytext=(1.03, 0.25))
      if signe > 0:
          fleches(ax, [(0.2, 0.75), (1.8, 1.25)])
      else:
          fleches(ax, [(0.8, 0.25), (1.2, 1.75)])
      ax.set_xlabel(r'$y$')
      ax.set_ylabel(r'$\dot y$')
      ax.set_ylim(-2.5, 2.5)
  ax1.set_title('état stationnaire stable')
  ax2.set_title('état stationnaire instable')
  fig.tight_layout()
  fig.savefig("edo-phase.svg", transparent=True)
#+end_src

#+CAPTION: *Diagramme de phase. À gauche $f$ est décroissante et l'état stationnaire est stable, à droite $f$ est croissante et l'état stationnaire est instable.*
#+LABEL: fig:phase
[[file:edo-phase.svg]]

$y_a^{\star}$ et $y_b^{\star}$ sont les états stationnaires, les points fixes de
la dynamique. Un seul regard sur ces graphiques nous renseigne sur leurs
propriétés\nbsp{}: $y_a^{\star}$ est stable (si pour une raison quelconque on
s'éloigne de $y_a^{\star}$, on y revient), contrairement à $y_b^{\star}$ qui est
instable (si l'on s'écarte de $y_b^{\star}$, on n'y revient jamais).

Ces graphiques suggèrent que si $f$ est monotone décroissante alors la
dynamique est stable, et que si $f$ est monotone croissante alors la dynamique
est instable. En fait nous devrions distinguer la stabilité globale et la
stabilité locale. La figure [[fig:phase-multiple][ci-dessous]] montre que la fonction $f$ peut
être beaucoup plus «\nbsp{}tordue\nbsp{}». Dans ce cas $y_a^{\star}$ et $y_c^{\star}$
sont des états stationnaires localement stables et $y_b^{\star}$ est un état
stationnaire instable. $y_a^{\star}$ est /localement/ stable, au sens où si
l'on s'écarte modérément de $y_a^{\star}$ on y revient, mais si l'on s'en
écarte trop, au-delà de $y_b^{\star}$, on ne revient jamais vers $y_a^{\star}$.

#+begin_src python :session edo1 :exports none :results none
  fm = lambda y: -0.5*(y-1)*(y-2)*(y-3)
  grille = np.linspace(0.45, 3.55, 300)
  fig, ax = plt.subplots(figsize=(7, 4))
  ax.plot(grille, fm(grille), 'b')
  ax.axhline(y=0, color='k', linewidth=0.8)
  for yst, nom in ((1, 'a'), (2, 'b'), (3, 'c')):
      ax.plot([yst], [0], 'ko', markersize=5)
      ax.annotate(r'$y_%s^\star$' % nom, xy=(yst, 0), xytext=(yst+0.04, 0.12))
  fleches(ax, [(0.5, 0.85), (1.5, 1.15), (2.5, 2.15), (3.5, 3.15)])
  ax.set_xlabel(r'$y$')
  ax.set_ylabel(r'$\dot y$')
  ax.set_ylim(-1.4, 1.4)
  fig.tight_layout()
  fig.savefig("edo-phase-multiple.svg", transparent=True)
#+end_src

#+CAPTION: *Diagramme de phase avec plusieurs états stationnaires. $y_a^{\star}$ et $y_c^{\star}$ sont localement stables, $y_b^{\star}$ est instable.*
#+LABEL: fig:phase-multiple
[[file:edo-phase-multiple.svg]]

On retient deux idées de ces exemples\nbsp{}:
1. l'état stationnaire, s'il existe, n'est pas toujours unique. On peut aussi
   imaginer des cas où il n'existe pas\nbsp{};
2. la stabilité, locale ou globale, ne peut apparaître que si la fonction
   $f$, dite /de transition/, est, au moins localement, décroissante.

Rétrospectivement, le point 2 est cohérent avec ce que nous avons vu dans le
cas linéaire. Soit l'équation différentielle linéaire $\dot y = b - ay$. Les
propriétés de stabilité dépendent du signe du paramètre \(a\)\nbsp{}: la dynamique
est stable si et seulement si $a$ est positif, c'est-à-dire s'il y a une
relation décroissante entre la variation et le niveau. Dans ce cas, pour
toute condition initiale, $y(t)$ converge vers $b/a$.

** Approche numérique
:PROPERTIES:
:CUSTOM_ID: numerique
:END:

L'approche décrite dans la section précédente peut ne pas suffire si l'on
désire obtenir des informations plus précises sur la solution. On recourt
alors au calcul numérique. Cette tâche est relativement simple avec un
logiciel scientifique\nbsp{}; nous utilisons ici la fonction =solve_ivp= de la
bibliothèque =scipy=, qui met en œuvre par défaut une méthode de Runge-Kutta
d'ordre 4(5) à pas adaptatif, la même que la fonction =ode45= de Matlab. Ces
fonctions résolvent des problèmes de la forme\nbsp{}:

\[
\dot y(t) = f\bigl(t, y(t)\bigr)
\]

où la condition initiale $y(t_0)$ est donnée, sur l'intervalle de temps
$[t_0,t_1]$.

À titre d'exemple nous allons résoudre numériquement l'équation
différentielle non linéaire de l'[[ex-bernoulli][exemple]] traité plus haut\nbsp{}; comme nous
disposons d'une solution analytique, nous pourrons évaluer la précision de la
solution numérique. Pour rappel, la solution analytique est\nbsp{}:

\[
y(t) = \frac{1}{e^{\frac{t^2}{2}}\left(\frac{1}{y(0)}-3\right) + 3}
\]

Avant de se lancer dans l'étude numérique de la dynamique, il convient de
discuter les propriétés de cette solution. On observe les points suivants.
1. $y^{\star} = 1/3$ est un état stationnaire. Pour le vérifier, on peut
   substituer $y(t) = 1/3$ dans l'équation différentielle. Si $y(0) = 1/3$
   alors $y(t) = 1/3$ pour tout $t$.
2. Si $y(0) < 1/3$ alors $\lim_{t\to\infty}y(t) = 0^+$. L'état stationnaire
   $y^{\star}$ n'est donc pas stable\nbsp{}; c'est zéro, l'autre état
   stationnaire, qui attire ces trajectoires.
3. Le cas $y(0) > 1/3$ mérite un traitement particulier. Le dénominateur
   $3 - e^{t^2/2}\left(3 - 1/y(0)\right)$ est alors une fonction décroissante
   de $t$, positive en $t=0$, qui s'annule en\nbsp{}:
   \[
   t^{\star} = \sqrt{2\log\frac{3}{3-\frac{1}{y(0)}}}
   \]
   La solution admet une asymptote verticale en \(t^{\star}\)\nbsp{}: pour tout
   $t < t^{\star}$, $y(t)$ est monotone croissante avec $\lim_{t\to t^{\star}}y(t) = +\infty$.
   La trajectoire cesse d'exister en $t^{\star}$. La formule reste définie
   au-delà, mais elle décrit alors une autre solution de l'équation, négative,
   croissante de $-\infty$ vers $0^-$, qui n'est pas reliée à la condition
   initiale.
Les propriétés de la solution dépendent radicalement de la condition
initiale. L'approche numérique du cas $y(0) < 1/3$ ne pose pas de problème,
contrairement au cas complémentaire, à cause de l'asymptote verticale.

Nous commençons par définir la fonction $f(t,y)$ associée à l'équation
différentielle, écrite sous la forme $\dot y(t) = 3t\,y(t)^2 - t\,y(t)$, puis la
solution exacte\nbsp{}:

#+begin_src python :session edo1 :exports code :results none
  def bernoulli(t, y):
      """Membre de droite de l'équation différentielle dy/dt = 3ty² - ty."""
      return 3*t*y**2 - t*y

  def exacte(t, y0):
      """Solution analytique pour la condition initiale y0."""
      return 1/(np.exp(t**2/2)*(1/y0 - 3) + 3)
#+end_src

Le script suivant résout l'équation pour deux conditions initiales, de part et
d'autre de l'état stationnaire. Dans le cas $y_0>y^{\star}$, l'instant terminal
est fixé en deçà de $t^{\star}$, car aucun solveur ne peut suivre la solution
jusqu'à une asymptote verticale. La fonction =solve_ivp= renvoie un objet
dont l'attribut =t= contient les instants où le solveur a évalué la solution,
et dont l'attribut =sol=, obtenu grâce à l'option =dense_output=, permet
d'interpoler la solution en tout point de l'intervalle.

#+begin_src python :session edo1 :exports code :results none
  ystar = 1/3
  solutions = {}
  for y0 in (0.9*ystar, 1.1*ystar):
      if y0 > ystar:
          tstar = np.sqrt(2*np.log(3/(3 - 1/y0)))
          t1 = 0.9*tstar
      else:
          t1 = 10.0
      solutions[y0] = solve_ivp(bernoulli, (0, t1), [y0], dense_output=True)
#+end_src

#+begin_src python :session edo1 :exports none :results none
  fig, axes = plt.subplots(2, 2, figsize=(11, 8))
  for j, (y0, sol) in enumerate(solutions.items()):
      t = np.linspace(0, sol.t[-1], 1000)
      y = sol.sol(t)[0]
      ax = axes[0, j]
      ax.plot(t, y, 'b', label='solution numérique')
      ax.plot(sol.t, sol.y[0], 'bo', markersize=3)
      ax.axhline(y=ystar, color='r', linewidth=1, linestyle='--')
      ax.annotate(r'$y^\star$', xy=(0.9*t[-1], ystar-0.03 if y0 < ystar else ystar+0.01))
      ax.set_xlabel(r'$t$')
      ax.set_ylabel(r'$y(t)$')
      ax.set_title(r'$y_0 = %s\,y^\star$' % ('0{,}9' if y0 < ystar else '1{,}1'))
      ax = axes[1, j]
      ax.semilogy(t, np.abs(y - exacte(t, y0)), 'b', label='tolérances par défaut')
      precis = solve_ivp(bernoulli, (0, sol.t[-1]), [y0], dense_output=True, rtol=1e-8, atol=1e-11)
      ax.semilogy(t, np.abs(precis.sol(t)[0] - exacte(t, y0)), 'g', label=r'$\mathtt{rtol}=10^{-8}$')
      ax.set_xlabel(r'$t$')
      ax.set_ylabel('erreur absolue')
      ax.legend(loc='center right')
  fig.tight_layout()
  fig.savefig("edo-bernoulli.svg", transparent=True)
#+end_src

#+CAPTION: *Solutions numériques. En haut, la solution interpolée (trait) et les points calculés par le solveur (ronds) pour $y_0 = 0{,}9\,y^{\star}$ (à gauche) et $y_0 = 1{,}1\,y^{\star}$ (à droite, l'asymptote verticale est en $t^{\star} = 2{,}1899$). En bas, l'erreur absolue par rapport à la solution analytique, avec les tolérances par défaut et avec une tolérance relative de $10^{-8}$.*
#+LABEL: fig:bernoulli
[[file:edo-bernoulli.svg]]

La figure [[fig:bernoulli][ci-dessus]] représente les résultats obtenus pour $y_0 = 0{,}9\,y^{\star}$
et $y_0 = 1{,}1\,y^{\star}$. Avec les tolérances par défaut, une tolérance
relative de $10^{-3}$ et une tolérance absolue de $10^{-6}$, le solveur se
contente de $29$ pas dans le premier cas et de $8$ pas dans le second\nbsp{}:
l'erreur est au plus de $4{,}3\times10^{-5}$ en valeur absolue dans le premier
cas, et de $5{,}7\times10^{-4}$ dans le second. On constate que les erreurs
sont quasi nulles lorsque la trajectoire de $y$ est plate et deviennent plus
importantes lorsque la trajectoire est plus pentue, en particulier à
l'approche de l'asymptote. Ces erreurs se réduisent en demandant au solveur
une approximation plus précise\nbsp{}: avec =rtol=1e-8= et =atol=1e-11=, le
nombre de pas passe à $127$ et $32$, et l'erreur maximale tombe à
$1{,}2\times10^{-8}$ et $3{,}8\times10^{-8}$. La précision se paie en
évaluations de la fonction =bernoulli=, ce qui est sans conséquence ici mais
peut ne pas l'être pour des systèmes de grande dimension.

Cette application nous apprend qu'il est toujours utile de s'interroger sur
les propriétés de la solution avant de se lancer dans un calcul numérique,
avec la solution analytique lorsqu'on en dispose (mais alors la solution
numérique devient sans intérêt, sauf à titre de vérification) ou avec
l'approche qualitative. Une réflexion préalable est toujours une bonne idée.
Dans le cas $y_0>y^{\star}$, le solveur échoue si l'on essaie d'obtenir la
solution pour $t\in[0,T]$ avec $T$ trop proche de $t^{\star}$, là où la
trajectoire est très pentue, ou a fortiori au-delà de \(t^{\star}\)\nbsp{}; et il
n'est pas garanti qu'il signale son échec autrement que par des valeurs
aberrantes.

#+BEGIN_remarque
Les solveurs à pas adaptatif comme =solve_ivp= choisissent eux-mêmes la
longueur des pas de temps de façon à respecter les tolérances demandées.
Les schémas à pas fixe, en particulier le schéma d'Euler, sont présentés dans
la note sur [[https://stephane-adjemian.fr/posts/modele-de-solow-en-temps-discret/][le modèle de Solow en temps discret]]\nbsp{}; la note sur la
[[https://stephane-adjemian.fr/posts/simulation-du-modele-de-solow/][simulation du modèle de Solow]] utilise la fonction =odeint= de =scipy=, qui
repose sur une autre famille de méthodes.
#+END_remarque

** Approximation locale
:PROPERTIES:
:CUSTOM_ID: locale
:END:

Comme dans la section sur l'[[#qualitative][approche qualitative]], on s'intéresse ici à une
dynamique autonome $\dot y = f(y)$, où $f$ est une fonction continue et
dérivable ne dépendant pas du temps. Jusqu'à présent nous n'avons considéré
que des approches globales. L'approche qualitative et l'approche analytique
sont globales, car nous caractérisons ou obtenons la solution pour toutes les
valeurs possibles de $t$ et indépendamment du niveau de $y$. L'approche
numérique est elle aussi globale, car elle est valable sur un intervalle de
valeurs de $t$ indépendamment du niveau de $y$. Ici nous abordons une approche
locale, au sens où ce que nous pourrons dire de $y$ ne sera pertinent que dans
un voisinage d'un niveau spécifique de la variable.

Si la fonction $f$ est dérivable en $\bar y$, un développement de Taylor à
l'ordre un donne\nbsp{}:

\[
\dot y = f(\bar y) + f'(\bar y)(y-\bar y) + O\bigl(|y-\bar y|^2\bigr)
\]

Si nous omettons le terme résiduel, qui tend vers zéro plus vite que
$|y-\bar y|$ lorsque $y$ se rapproche de $\bar y$, nous avons\nbsp{}:

\[
\dot y \approx f(\bar y) + f'(\bar y)(y-\bar y)
\]

Dans un voisinage de $\bar y$, la dynamique est approximativement linéaire.
La qualité de l'approximation dépend de la distance de $y$ à $\bar y$. Nous
savons résoudre cette équation différentielle, il s'agit d'une équation à
coefficients constants. Mais il faut garder à l'esprit que la solution que
nous obtiendrons en appliquant les résultats de la [[#constants][première section]] ne sera
acceptable que pour des valeurs de $y$ dans un voisinage de $\bar y$. C'est en
ce sens que cette approche est locale.

A priori nous pouvons choisir n'importe quel point $\bar y$, pourvu que la
fonction $f$ y soit dérivable. Habituellement on approche le modèle autour
d'un état stationnaire de la dynamique, $y^{\star}$. Cela permet d'éliminer la
constante, puisque par définition $\dot y$ est nul à l'état stationnaire.
Ainsi nous avons\nbsp{}:

\[
\dot y \approx f'(y^{\star})(y-y^{\star})
\]

Une autre motivation est que si cet état stationnaire est stable, alors nous
sommes sûrs que si $y(0)$ est dans un voisinage de $y^{\star}$, $y(t)$ restera
dans ce même voisinage, ce qui nous assure de la qualité de l'approximation.
Notons que dans le cas d'états stationnaires multiples, il y a un degré de
liberté sur le choix du point au voisinage duquel $f$ est approchée. Les
propriétés de la solution approchée peuvent être fort différentes, par
exemple en termes de stabilité.

L'interprétation géométrique de cette approximation est directe. Faire une
approximation à l'ordre un, on parle aussi de /linéarisation/, c'est
remplacer la fonction $f$ par sa tangente en $y^{\star}$. La figure
[[fig:linearisation][ci-dessous]] illustre ce remplacement, pour la fonction de transition du
modèle de Solow que nous retrouverons dans la [[#solow][section suivante]].

#+begin_src python :session edo1 :exports none :results none
  s, al, mu = 0.20, 0.36, 0.14
  fs = lambda k: s*k**al - mu*k
  kst = (s/mu)**(1/(1-al))
  pente = mu*(al-1)
  grille = np.linspace(0, 5, 300)
  fig, ax = plt.subplots(figsize=(7, 4))
  ax.plot(grille, fs(grille), 'b', label=r'$f(y)$')
  ax.plot(grille, pente*(grille-kst), 'r', linewidth=1, linestyle='--', label=r"tangente en $y^\star$")
  ax.axhline(y=0, color='k', linewidth=0.8)
  ax.plot([kst], [0], 'ko', markersize=5)
  ax.annotate(r'$y^\star$', xy=(kst, 0), xytext=(kst+0.08, 0.02))
  ax.plot([0], [-pente*kst], 'ro', markersize=4)
  ax.annotate(r"$-f'(y^\star)\,y^\star$", xy=(0, -pente*kst), xytext=(0.15, -pente*kst+0.01), color='r')
  ax.set_xlabel(r'$y$')
  ax.set_ylabel(r'$\dot y$')
  ax.legend(loc='lower left')
  fig.tight_layout()
  fig.savefig("edo-linearisation.svg", transparent=True)
#+end_src

#+CAPTION: *Interprétation graphique de la linéarisation. La fonction de transition $f(y) = sy^{\alpha} - (n+g+\delta)y$ du modèle de Solow et sa tangente à l'état stationnaire.*
#+LABEL: fig:linearisation
[[file:edo-linearisation.svg]]

Pour finir, notons qu'a priori rien ne nous empêche de considérer des ordres
d'approximation supérieurs à un. Dans certains cas, il peut être nécessaire,
afin de ne pas omettre des propriétés intéressantes de la variable étudiée,
d'aller chercher une approximation à l'ordre deux ou trois. Le développement
à l'ordre deux autour de l'état stationnaire\nbsp{}:

\[
\dot y \approx f'(y^{\star})(y-y^{\star}) + \frac12f''(y^{\star})(y-y^{\star})^2
\]

est une équation de Bernoulli à coefficients constants, avec $m=2$, dont nous
savons donc exhiber la solution explicite. C'est ce que nous ferons pour le
modèle de Solow.

** Application : le modèle de Solow
:PROPERTIES:
:CUSTOM_ID: solow
:END:

Le modèle de Solow est un modèle dynamique en temps continu décrivant
l'évolution du stock de capital physique, ou de la production, dans une
économie fermée. Nous ne l'aborderons pas ici en détail, il est présenté dans
la note sur [[https://stephane-adjemian.fr/posts/modele-de-solow/][le modèle de Solow]]\nbsp{}; le but est simplement d'appliquer les
méthodes présentées plus haut.

On se donne une fonction de production néoclassique[fn:: Une fonction
homogène de degré un, croissante et concave en chacun de ses arguments, et
vérifiant les conditions d'Inada. Voir la note sur le modèle de Solow.],
$Y(t) = F\bigl(K(t), A(t)L(t)\bigr)$, avec $L(t) = e^{nt}$ la population et
$A(t) = e^{gt}$ l'indice d'efficience du travail[fn:: Les autres notes de ce
site notent $x$ le taux de croissance du progrès technique.]. La loi d'évolution du stock de
capital physique est\nbsp{}:

\[
\dot K(t) = sY(t) - \delta K(t)
\]

avec $s\in\,]0,1[$ le taux d'épargne exogène et $\delta\in[0,1]$ le taux de
dépréciation du capital. Cette équation nous dit simplement que le stock de
capital augmente si et seulement si l'investissement domine la
dépréciation. En substituant la fonction de production, on voit bien qu'il
s'agit d'une équation différentielle non linéaire du premier ordre\nbsp{}:

\[
\dot K(t) = sF\bigl(K(t), e^{(n+g)t}\bigr) - \delta K(t)
\]

Cette équation différentielle est non autonome\nbsp{}: la relation entre la
variation du stock de capital et son niveau dépend du temps, via le terme
$e^{(n+g)t}$ qui représente la croissance du travail efficace. La première
chose à faire est de nous ramener à une équation différentielle autonome.
Pour cela il suffit d'éliminer les tendances démographique et technologique.
On pose $\hat k(t) = K(t)/\bigl(A(t)L(t)\bigr)$, le stock de capital par tête
efficace, «\nbsp{}purgé\nbsp{}» de la croissance de la population et du progrès
technique. On vérifie facilement que la dynamique du capital par tête
efficace est donnée par\nbsp{}:

\[
\dot{\hat k}(t) = sf\bigl(\hat k(t)\bigr) - (n+g+\delta)\,\hat k(t)
\]

où $f(\hat k) = F(\hat k, 1) = \hat y$ est la production par tête efficace. Il
s'agit bien d'une équation différentielle autonome. Si la fonction de
production est néoclassique, il existe un unique état stationnaire
$\hat k^{\star}$ strictement positif. En passant, notons que $\hat k=0$ est
aussi un état stationnaire, dit trivial, que nous laisserons de côté. Pour se
convaincre de ce résultat d'existence et d'unicité, il suffit de reprendre
l'équation sous la forme\nbsp{}:

\[
\frac{\dot{\hat k}(t)}{\hat k(t)} = s\,\frac{f\bigl(\hat k(t)\bigr)}{\hat k(t)} - (n+g+\delta)
\]

Le taux de croissance du capital par tête efficace est strictement positif si
et seulement si l'investissement brut par unité de capital est strictement
supérieur au taux de dépréciation effectif $n+g+\delta$. Le premier terme du
membre de droite est égal à la productivité moyenne du capital multipliée par
le taux d'épargne. Il est strictement positif, monotone décroissant, parce
que les rendements du capital sont décroissants, et tend vers l'infini en zéro
et vers zéro à l'infini par les conditions d'Inada. Ainsi la courbe
représentative de l'investissement brut par unité de capital croise
nécessairement une seule fois la droite horizontale d'ordonnée
\(n+g+\delta\)\nbsp{}: l'état stationnaire est unique, et il vérifie\nbsp{}:

\[
\frac{\hat y^{\star}}{\hat k^{\star}} = \frac{n+g+\delta}{s}
\]

Dans la suite nous chercherons à décrire la dynamique du produit par tête
efficace, et plus spécialement la dynamique d'ajustement vers l'état
stationnaire, en calculant la /vitesse de convergence/. En dérivant
$\hat y(t) = f(\hat k(t))$ par rapport à $t$, on obtient\nbsp{}:

\[
\frac{\dot{\hat y}(t)}{\hat y(t)} = s\,r\bigl(\hat k(t)\bigr) - (n+g+\delta)\,\alpha\bigl(\hat k(t)\bigr)
\]

avec $r(\hat k) = f'(\hat k)$ la productivité marginale du capital, qui dans un
environnement parfaitement concurrentiel correspond au taux d'intérêt réel,
et $\alpha(\hat k) = \hat kf'(\hat k)/f(\hat k)$ l'élasticité du produit par
rapport au capital, ou encore, dans un environnement parfaitement
concurrentiel, la part des revenus du capital dans le revenu total. Par
définition de $\alpha$, l'état stationnaire vérifie aussi
$\hat y^{\star}/\hat k^{\star} = r^{\star}/\alpha^{\star}$. Enfin, nous
utiliserons l'élasticité de substitution entre le capital et le travail,
définie [[https://stephane-adjemian.fr/posts/fonction-de-production-ces/][ici]], qui s'écrit en fonction de la technologie intensive\nbsp{}:

\[
\sigma(\hat k) = -\frac{f'(\hat k)\bigl[f(\hat k) - \hat kf'(\hat k)\bigr]}{\hat k\,f(\hat k)\,f''(\hat k)}
\qquad\Leftrightarrow\qquad
\frac{\hat kf''(\hat k)}{f'(\hat k)} = -\frac{1-\alpha(\hat k)}{\sigma(\hat k)}
\]

Nous supposerons parfois que la technologie est de type CES (pour /constant
elasticity of substitution/). Cette fonction de production est plus générale
que la fonction Cobb-Douglas, mais elle n'est pas néoclassique, les
conditions d'Inada n'étant pas satisfaites. Dans ce cas l'existence de l'état
stationnaire $\hat k^{\star}>0$ n'est plus assurée et dépend des valeurs des
paramètres. On posera\nbsp{}:

\[
\hat y(t) = \left(\gamma_1\,\hat k(t)^{\rho} + \gamma_2\right)^{\frac1\rho}
\]

avec $\rho\in\,]-\infty,1]$, $\sigma = (1-\rho)^{-1}$, soit $\rho = 1-1/\sigma$,
l'élasticité de substitution constante entre les facteurs, et
$\gamma_1+\gamma_2=1$. On retrouve la fonction Cobb-Douglas, qui est
néoclassique, lorsque $\rho\to0$, c'est-à-dire $\sigma\to1$. La productivité
moyenne $f(\hat k)/\hat k = \left(\gamma_1 + \gamma_2\hat k^{-\rho}\right)^{1/\rho}$
est toujours décroissante, mais elle ne parcourt plus $]0,\infty[$ tout
entier, d'où les conditions d'existence de l'état stationnaire\nbsp{}:
1. si $\rho\in\,]0,1]$, les facteurs étant plus substituables que dans le cas
   Cobb-Douglas, la productivité moyenne décroît de l'infini vers
   $\gamma_1^{1/\rho}$, et l'unique état stationnaire $\hat k^{\star}>0$
   existe si et seulement si $s\gamma_1^{1/\rho} < n+g+\delta$, c'est-à-dire
   si le taux d'épargne n'est pas trop élevé\nbsp{};
2. si $\rho < 0$, les facteurs étant moins substituables que dans le cas
   Cobb-Douglas, la productivité moyenne décroît de $\gamma_1^{1/\rho}$ vers
   zéro, et l'unique état stationnaire $\hat k^{\star}>0$ existe si et
   seulement si $s\gamma_1^{1/\rho} > n+g+\delta$, c'est-à-dire si le taux
   d'épargne n'est pas trop faible\nbsp{};
3. si $\rho=0$ la fonction de production est Cobb-Douglas, et l'existence et
   l'unicité de l'état stationnaire sont assurées indépendamment des valeurs
   des paramètres.
Avec cette technologie, le taux d'intérêt et la part des revenus du capital
s'écrivent en fonction de \(\hat y\)\nbsp{}:

\[
r(\hat y) = \gamma_1^{\frac1\rho}\left(1-\gamma_2\,\hat y^{-\rho}\right)^{-\frac{1-\rho}{\rho}}
\qquad\text{et}\qquad
\alpha(\hat y) = 1-\gamma_2\,\hat y^{-\rho}
\]

La part du capital est une fonction croissante du produit par tête efficace
si et seulement si $\rho>0$, c'est-à-dire si et seulement si l'élasticité de
substitution est supérieure à un. À l'état stationnaire, on a
$\alpha^{\star} = \gamma_1\left(s/(n+g+\delta)\right)^{\rho}$.

Dans toute la suite nous utilisons l'étalonnage des notes précédentes\nbsp{}:
$s=0{,}20$, $n=0{,}02$, $g=0{,}02$, $\delta=0{,}10$, et une part du capital à
l'état stationnaire $\alpha^{\star}=0{,}36$. Pour comparer les technologies à
état stationnaire donné, nous fixons $\gamma_1$ de façon que
$\alpha^{\star}=0{,}36$ quelle que soit l'élasticité de substitution\nbsp{}:
$\gamma_1 = \alpha^{\star}\left((n+g+\delta)/s\right)^{\rho}$ et
$\gamma_2 = 1-\gamma_1$.

#+begin_src python :session edo1 :exports code :results none
  s, n, g, delta, alphastar = 0.20, 0.02, 0.02, 0.10, 0.36
  mu = n + g + delta

  def technologie(sigma):
      """Fonction de production CES intensive, sa dérivée, l'élasticité du produit
      par rapport au capital et l'état stationnaire, pour une élasticité de
      substitution sigma, avec alpha* = 0,36 quelle que soit sigma."""
      rho = 1 - 1/sigma
      if abs(rho) < 1e-12:                     # Cobb-Douglas
          f = lambda k: k**alphastar
          fprime = lambda k: alphastar*k**(alphastar-1)
          kstar = (s/mu)**(1/(1-alphastar))
      else:
          g1 = alphastar*(mu/s)**rho
          g2 = 1 - g1
          f = lambda k: (g1*k**rho + g2)**(1/rho)
          fprime = lambda k: g1*k**(rho-1)*(g1*k**rho + g2)**(1/rho-1)
          kstar = ((mu/s)**rho*(1-alphastar)/g2)**(-1/rho)
      alpha = lambda k: k*fprime(k)/f(k)
      return f, fprime, alpha, kstar

  def transition(sigma, k0, T=80.0):
      """Trajectoire exacte du capital par tête efficace, par intégration numérique."""
      f, fprime, alpha, kstar = technologie(sigma)
      dk = lambda t, k: s*f(k) - mu*k
      return solve_ivp(dk, (0, T), [k0], dense_output=True, rtol=1e-10, atol=1e-12)
#+end_src

*** Caractérisation de l'ajustement sans approximation
:PROPERTIES:
:CUSTOM_ID: sans-approximation
:END:

La vitesse de convergence est définie comme l'opposé du taux de croissance du
taux de croissance du produit par tête efficace\nbsp{}:

\[
\beta(t) = -\frac{\dot g_{\hat y}(t)}{g_{\hat y}(t)}
\qquad\text{avec}\qquad
g_{\hat y}(t) = \frac{\dot{\hat y}(t)}{\hat y(t)}
\]

Elle mesure la vitesse à laquelle le taux de croissance s'éteint, à mesure
que l'économie s'approche de son état stationnaire. Il sera commode, ici et
dans les sections suivantes, de raisonner sur l'écart logarithmique du
capital à son niveau stationnaire, $u(t) = \log\bigl(\hat k(t)/\hat k^{\star}\bigr)$,
dont la dynamique s'écrit\nbsp{}:

\[
\dot u = G(u) \equiv s\,\frac{f(\hat k^{\star}e^{u})}{\hat k^{\star}e^{u}} - (n+g+\delta)
\]

avec $G(0)=0$. Deux dérivées reviendront sans cesse. D'abord, en notant que
\(\mathrm d\hat k/\mathrm du = \hat k\)\nbsp{}:

\[
G'(u) = \hat k\,\frac{\mathrm d}{\mathrm d\hat k}\left[s\frac{f(\hat k)}{\hat k}\right]
= s\left[f'(\hat k) - \frac{f(\hat k)}{\hat k}\right]
= -s\frac{f(\hat k)}{\hat k}\bigl[1-\alpha(\hat k)\bigr]
\]

Ensuite, en dérivant $\log\alpha = \log\hat k + \log f' - \log f$ par rapport
à $u$ et en utilisant la définition de l'élasticité de substitution\nbsp{}:

\[
\frac{\mathrm d\alpha}{\mathrm du} = \alpha\left[1 + \frac{\hat kf''}{f'} - \alpha\right]
= \alpha(\hat k)\bigl[1-\alpha(\hat k)\bigr]\left[1-\frac{1}{\sigma(\hat k)}\right]
\]

#+BEGIN_property
La vitesse de convergence du produit par tête efficace est, à tout instant de
la transition\nbsp{}:
\[
\beta(t) = (n+g+\delta)\bigl[1-\alpha(\hat k)\bigr]\left[1 + \frac{1}{\sigma(\hat k)}\left(\frac{\hat y/\hat y^{\star}}{\hat k/\hat k^{\star}} - 1\right)\right]
\]
#+END_property

#+BEGIN_proof
Le taux de croissance du produit est $g_{\hat y} = \alpha(\hat k)\,g_{\hat k}$ et
$g_{\hat k} = \dot u = G(u)$. En dérivant par rapport au temps,
$\dot g_{\hat y} = \left[\frac{\mathrm d\alpha}{\mathrm du}G(u) + \alpha G'(u)\right]G(u)$,
d'où\nbsp{}:
\[
\beta = -G'(u) - \frac{1}{\alpha}\frac{\mathrm d\alpha}{\mathrm du}G(u)
= s\frac{f}{\hat k}(1-\alpha) - (1-\alpha)\left(1-\frac1\sigma\right)\left(s\frac{f}{\hat k} - (n+g+\delta)\right)
\]
soit, en regroupant les termes en \(sf/\hat k\)\nbsp{}:
\[
\beta = (1-\alpha)\left[\frac{1}{\sigma}\,s\frac{f}{\hat k} + (n+g+\delta)\left(1-\frac1\sigma\right)\right]
= (n+g+\delta)(1-\alpha)\left[1 + \frac1\sigma\left(\frac{sf/\hat k}{n+g+\delta} - 1\right)\right]
\]
Il reste à noter que, d'après la condition d'état stationnaire,
$\frac{sf(\hat k)/\hat k}{n+g+\delta} = \frac{\hat y/\hat k}{\hat y^{\star}/\hat k^{\star}}$.
#+END_proof

Dans le cas d'une technologie Cobb-Douglas, l'élasticité de substitution est
unitaire et $\hat y/\hat y^{\star} = (\hat k/\hat k^{\star})^{\alpha}$, de sorte
que l'on retrouve le résultat de Barro et Sala-i-Martin (1995, annexe du
chapitre 1)\nbsp{}:

\[
\beta_1(t) = (n+g+\delta)(1-\alpha)\left(\frac{\hat y(t)}{\hat y^{\star}}\right)^{-\frac{1-\alpha}{\alpha}} > 0
\]

La vitesse de convergence diminue le long de la transition si et seulement
si l'économie rejoint son état stationnaire par dessous, c'est-à-dire si
$g_{\hat y}>0$. À l'état stationnaire on retrouve le résultat standard, établi
dans la [[#ordre-1][section suivante]] en linéarisant le modèle\nbsp{}:
$\beta_1(t)\to\beta^{\star} = (1-\alpha)(n+g+\delta)$. La note sur la
[[https://stephane-adjemian.fr/posts/simulation-du-modele-de-solow/][simulation du modèle de Solow]] représente cette vitesse le long de la
transition.

Plus généralement, lorsque la fonction de production est de type CES, on a
$\hat y/\hat k = \gamma_1^{1/\rho}\alpha(\hat y)^{-1/\rho}$, puisque
$\alpha = \gamma_1(\hat k/\hat y)^{\rho}$, et la vitesse de convergence
s'exprime en fonction de la seule part du capital\nbsp{}:

\[
\beta_{\sigma}(t) = (n+g+\delta)\bigl[1-\alpha(\hat y)\bigr]\left[1 + \frac1\sigma\left(\left(\frac{\alpha(\hat y)}{\alpha^{\star}}\right)^{\frac{\sigma}{1-\sigma}} - 1\right)\right]
\]

avec $\alpha(\hat y) = 1-\gamma_2\hat y^{-\rho}$. Le sens de variation de la
part du capital au cours de la transition dépend du signe de \(\rho\)\nbsp{}: elle
augmente avec le produit si $\sigma>1$ et diminue si $\sigma < 1$. Il serait
tentant d'en déduire le sens de variation de la vitesse de convergence à
partir du seul facteur \(1-\alpha(\hat y)\)\nbsp{}: une part du capital croissante
ralentirait la convergence, une part décroissante l'accélérerait, de sorte que le sens de
variation de $\beta_{\sigma}$ dépendrait du signe de $\rho$. Mais le terme
entre crochets varie lui aussi, en sens inverse, et c'est lui qui l'emporte.

#+BEGIN_property
Au voisinage de l'état stationnaire, la vitesse de convergence varie avec
l'écart $u = \log(\hat k/\hat k^{\star})$ selon\nbsp{}:
\[
\left.\frac{\mathrm d\beta_{\sigma}}{\mathrm du}\right|_{u=0} = -(n+g+\delta)(1-\alpha^{\star})\left[\alpha^{\star} + \frac{1-2\alpha^{\star}}{\sigma}\right]
\]
Cette dérivée est négative dès que $\alpha^{\star} < 1/2$, quelle que soit
l'élasticité de substitution. Une économie située au-dessous de son état
stationnaire converge donc plus vite que $\beta^{\star}$, et sa vitesse de
convergence diminue au cours de la transition\nbsp{}; une économie située
au-dessus converge moins vite que $\beta^{\star}$, et sa vitesse augmente.
#+END_property

#+BEGIN_proof
Écrivons $\beta = (1-\alpha)\left[\frac{\phi}{\sigma} + (n+g+\delta)\left(1-\frac1\sigma\right)\right]$
avec $\phi = sf(\hat k)/\hat k$. On a $\mathrm d\phi/\mathrm du = G'(u) = -\phi(1-\alpha)$ et
$\mathrm d(1-\alpha)/\mathrm du = -\alpha(1-\alpha)(1-1/\sigma)$. À l'état stationnaire
$\phi = n+g+\delta$, et le crochet vaut $n+g+\delta$, d'où\nbsp{}:
\[
\frac{\mathrm d\beta}{\mathrm du} = -\alpha(1-\alpha)\left(1-\frac1\sigma\right)(n+g+\delta) - (1-\alpha)\frac{(n+g+\delta)(1-\alpha)}{\sigma}
= -(n+g+\delta)(1-\alpha)\left[\alpha + \frac{1-2\alpha}{\sigma}\right]
\]
Le premier terme est l'effet de la part du capital, positif lorsque
\(\sigma < 1\)\nbsp{}; le second est l'effet du crochet, toujours négatif, et il
domine le premier si et seulement si $1-2\alpha+\alpha\sigma>0$, ce qui est
acquis pour $\alpha < 1/2$.
#+END_proof

La figure [[fig:vitesse][ci-dessous]] confirme ce résultat, et montre ce qui se passe loin de
l'état stationnaire. Les trois technologies partagent le même état
stationnaire et la même vitesse asymptotique $\beta^{\star} = 0{,}0896$, par
construction de l'étalonnage. Par dessous, la vitesse de convergence part
d'un niveau d'autant plus élevé que les facteurs sont substituables, et
décroît vers \(\beta^{\star}\)\nbsp{}; pour $\sigma=0{,}5$ elle est presque plate
pendant les premières années, l'effet de la part du capital compensant
alors presque exactement celui du crochet. Par dessus, la vitesse croît vers
$\beta^{\star}$ dans les trois cas. Pour $\sigma=0{,}5$ elle est même négative
pendant les premières années\nbsp{}: la production décroît de plus en plus
vite avant de ralentir. Rien de paradoxal, puisque $g_{\hat y} = \alpha\,g_{\hat k}$
et que, loin au-dessus de l'état stationnaire, la part du capital est très
faible ($0{,}12$ contre $0{,}36$ à l'état stationnaire) et remonte vite,
alors que le taux de décroissance du capital ne s'atténue que lentement.

#+begin_src python :session edo1 :exports none :results none
  def vitesse(sigma, k):
      """Vitesse de convergence du produit par tête efficace le long d'une trajectoire."""
      f, fprime, alpha, kstar = technologie(sigma)
      return mu*(1-alpha(k))*(1 + (1/sigma)*((f(k)/f(kstar))/(k/kstar) - 1))

  betastar = mu*(1-alphastar)
  tg = np.linspace(0, 60, 400)
  fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
  for ax, ratio in ((ax1, 0.25), (ax2, 4.0)):
      for sigma, couleur in ((0.5, 'r'), (1.0, 'b'), (2.0, 'g')):
          kstar = technologie(sigma)[3]
          k = transition(sigma, ratio*kstar).sol(tg)[0]
          ax.plot(tg, vitesse(sigma, k), couleur, label=r'$\sigma = %s$' % str(sigma).replace('.', '{,}'))
      ax.axhline(y=betastar, color='k', linewidth=0.8, linestyle=':')
      ax.annotate(r'$\beta^\star$', xy=(57, betastar+0.004))
      ax.set_xlabel(r'$t$ (années)')
      ax.set_ylabel(r'$\beta(t)$')
      ax.legend()
  ax1.set_title(r'départ à $0{,}25\,\hat k^\star$')
  ax2.set_title(r'départ à $4\,\hat k^\star$')
  fig.tight_layout()
  fig.savefig("solow-vitesse.svg", transparent=True)
#+end_src

#+CAPTION: *Vitesse de convergence du produit par tête efficace le long de la transition, pour trois élasticités de substitution, en partant au-dessous (à gauche) et au-dessus (à droite) de l'état stationnaire.*
#+LABEL: fig:vitesse
[[file:solow-vitesse.svg]]

On retient que la vitesse d'ajustement vers l'état stationnaire n'est pas
constante, mais change le long de la transition, et que le sens de variation
dépend d'abord de la position initiale vis-à-vis de l'état stationnaire\nbsp{};
l'élasticité de substitution entre les facteurs module l'ampleur de ces
variations, et peut en changer la forme loin de l'état stationnaire. Dans
cette section nous n'avons pas résolu d'équation différentielle, pour
calculer par exemple le niveau du produit par tête efficace à l'instant $t$.
Cela n'est généralement pas possible[fn:: Sauf pour certaines fonctions de
production. Dans le cas d'une technologie Cobb-Douglas ou Leontief il est
possible de mener les calculs jusqu'au bout\nbsp{}; nous y reviendrons.], et il
faut alors recourir à une approche numérique, comme nous l'avons fait pour
tracer la figure, ou considérer une approximation de l'équation
différentielle. Notons néanmoins que, sans même chercher à résoudre
l'équation différentielle, nous avons pu apprendre des choses relativement
précises sur la dynamique de transition.

*** Caractérisation de l'ajustement avec une linéarisation de la dynamique
:PROPERTIES:
:CUSTOM_ID: ordre-1
:END:

Appliquons maintenant l'[[#locale][approximation locale]] à la dynamique de l'écart
$u = \log(\hat k/\hat k^{\star})$, dont l'état stationnaire est $u=0$.

#+BEGIN_property
Au voisinage de l'état stationnaire, la dynamique du capital par tête
efficace est approximativement linéaire\nbsp{}:
\[
\dot u \approx -\beta\,u
\qquad\text{avec}\qquad
\beta = (n+g+\delta)(1-\alpha^{\star})
\]
et il en va de même du produit par tête efficace\nbsp{}:
\[
\frac{\mathrm d\log\hat y(t)}{\mathrm dt} \approx -\beta\log\frac{\hat y(t)}{\hat y^{\star}}
\]
Dans le cas d'une fonction de production CES,
$\beta = (n+g+\delta)\left[1-\gamma_1\left(\frac{s}{n+g+\delta}\right)^{\rho}\right]$.
#+END_property

#+BEGIN_proof
Le développement de Taylor à l'ordre un de $G$ autour de zéro donne
$\dot u \approx G'(0)\,u$, et $G'(0) = -s\frac{f(\hat k^{\star})}{\hat k^{\star}}(1-\alpha^{\star}) = -(n+g+\delta)(1-\alpha^{\star})$
d'après la condition d'état stationnaire. Pour le produit, le développement
de $\log f(\hat k^{\star}e^{u})$ à l'ordre un donne
$\log(\hat y/\hat y^{\star}) \approx \alpha^{\star}u$, puisque la dérivée de
$\log f(\hat k^{\star}e^{u})$ par rapport à $u$ est \(\alpha(\hat k)\)\nbsp{}; ainsi
$\frac{\mathrm d}{\mathrm dt}\log\hat y = \alpha(\hat k)\dot u \approx -\beta\alpha^{\star}u \approx -\beta\log(\hat y/\hat y^{\star})$.
Enfin, pour la CES, $\alpha^{\star} = \gamma_1(\hat k^{\star}/\hat y^{\star})^{\rho} = \gamma_1\left(s/(n+g+\delta)\right)^{\rho}$.
#+END_proof

En considérant un développement de Taylor d'ordre un autour de l'état
stationnaire, nous perdons de l'information sur la dynamique de
transition\nbsp{}: ici la vitesse d'ajustement vers l'état stationnaire est
constante. On note néanmoins que la vitesse de convergence obtenue ici est la
limite, lorsque $t$ tend vers l'infini, de la vitesse de convergence obtenue
dans la section précédente. Cette équivalence asymptotique n'est pas
étonnante, puisque la stabilité de l'état stationnaire nous assure que,
lorsque $t$ devient assez grand, le produit par tête efficace se trouve
arbitrairement proche de l'état stationnaire. Remarquons enfin que nous avons
pu dérouler les calculs sans spécifier la fonction de production\nbsp{}: c'est
généralement le cas lorsque l'on considère une approximation du modèle.

Nous avons perdu une partie des propriétés de la transition, mais nous
pouvons maintenant résoudre la dynamique, c'est-à-dire calculer le niveau du
produit par tête efficace à un instant quelconque. Évidemment, cette
résolution ne sera valable que dans un voisinage de l'état stationnaire.
L'équation linéaire à coefficient constant nous dit que le taux de croissance
de la distance à l'état stationnaire est négatif et constant\nbsp{}: pour toute
condition initiale, la distance se résorbe en un temps infini, $u(t) = u(0)e^{-\beta t}$,
et de même pour le produit. En substituant la définition de l'écart, on
obtient\nbsp{}:

\[
\hat y(t) = \hat y^{\star\,\omega(t)}\,\hat y(0)^{1-\omega(t)}
\qquad\text{avec}\qquad
0\leq\omega(t) = 1-e^{-\beta t}\xrightarrow[t\to\infty]{}1
\]

Le produit par tête efficace à l'instant $t$ est une moyenne géométrique
pondérée de sa valeur initiale et de sa valeur stationnaire, le poids de
cette dernière tendant vers un. Pour notre étalonnage, $\beta^{\star} = 0{,}0896$
et la demi-vie de l'écart à l'état stationnaire vaut $\log 2/\beta^{\star} = 7{,}7$
années.

*** Caractérisation de l'ajustement avec une approximation à l'ordre deux
:PROPERTIES:
:CUSTOM_ID: ordre-2
:END:

On pourrait développer à l'ordre deux la dynamique du produit, comme on
vient de le faire à l'ordre un, mais la dérivée seconde est alors très
lourde et ne se prête pas à une résolution. En raisonnant sur l'écart du
capital $u$ plutôt que sur celui du produit, le calcul tient en quelques
lignes, et l'équation approchée se résout explicitement.

#+BEGIN_property
Au voisinage de l'état stationnaire, la dynamique de l'écart
$u = \log(\hat k/\hat k^{\star})$ est, à l'ordre deux\nbsp{}:
\[
\dot u \approx -\beta\,u + \gamma\,u^2
\qquad\text{avec}\qquad
\gamma = \frac{\beta}{2}\left(1-\frac{\alpha^{\star}}{\sigma^{\star}}\right)
\]
où $\sigma^{\star} = \sigma(\hat k^{\star})$ est l'élasticité de substitution à
l'état stationnaire.
#+END_property

#+BEGIN_proof
Le développement de Taylor à l'ordre deux de $G$ autour de zéro s'écrit
$\dot u \approx G'(0)u + \frac12G''(0)u^2$. Nous connaissons $G'(0) = -\beta$. En
dérivant $G'(u) = -\phi(\hat k)\bigl[1-\alpha(\hat k)\bigr]$, avec
$\phi = sf/\hat k$, il vient\nbsp{}:
\[
G''(u) = -\frac{\mathrm d\phi}{\mathrm du}(1-\alpha) + \phi\,\frac{\mathrm d\alpha}{\mathrm du}
= \phi(1-\alpha)^2 + \phi\,\alpha(1-\alpha)\left(1-\frac1\sigma\right)
= \phi(1-\alpha)\left[1-\frac{\alpha}{\sigma}\right]
\]
À l'état stationnaire, $\phi = n+g+\delta$ et
$G''(0) = (n+g+\delta)(1-\alpha^{\star})\left(1-\alpha^{\star}/\sigma^{\star}\right) = 2\gamma$.
#+END_proof

Le terme quadratique est absent si et seulement si $\sigma^{\star}=\alpha^{\star}$.
Pour une technologie Cobb-Douglas, \(\gamma = \frac{\beta}{2}(1-\alpha)>0\)\nbsp{}; plus
généralement, $\gamma$ est positif dès que $\sigma^{\star}>\alpha^{\star}$, ce
qui est le cas de tous les étalonnages usuels. L'équation approchée est une
[[#bernoulli][équation de Bernoulli]] à coefficients constants, avec $m=2$, $R=\beta$ et
$Q=\gamma$, et nous savons la résoudre.

#+BEGIN_property
La solution de l'approximation à l'ordre deux est\nbsp{}:
\[
u(t) = \frac{u(0)\,e^{-\beta t}}{1 - \frac{\gamma}{\beta}\,u(0)\left(1-e^{-\beta t}\right)}
\]
et le produit par tête efficace s'en déduit, au même ordre, par\nbsp{}:
\[
\log\frac{\hat y(t)}{\hat y^{\star}} \approx \alpha^{\star}u(t) + \frac12\alpha^{\star}(1-\alpha^{\star})\left(1-\frac{1}{\sigma^{\star}}\right)u(t)^2
\]
#+END_property

#+BEGIN_proof
On divise l'équation par $u^2$ et l'on pose $z = 1/u$, de sorte que
$\dot z = -\dot u/u^2$. Il vient $\dot z = \beta z - \gamma$, équation linéaire
à coefficients constants dont l'état stationnaire est $\gamma/\beta$ et la
solution $z(t) = \left(z(0) - \frac{\gamma}{\beta}\right)e^{\beta t} + \frac{\gamma}{\beta}$.
En revenant à $u = 1/z$ avec $z(0) = 1/u(0)$, et en multipliant numérateur
et dénominateur par $u(0)e^{-\beta t}$, on obtient l'expression annoncée.
Pour le produit, la dérivée seconde de $\log f(\hat k^{\star}e^{u})$ par
rapport à $u$ est $\mathrm d\alpha/\mathrm du = \alpha(1-\alpha)(1-1/\sigma)$, ce qui
donne le développement à l'ordre deux.
#+END_proof

On peut aller un cran plus loin et écrire directement, à l'ordre deux, la
dynamique de l'écart du produit $v = \log(\hat y/\hat y^{\star})$. Elle est de
la même forme que celle du capital, et sa vitesse de convergence se compare
directement à celle de la section [[#sans-approximation][sans approximation]].

#+BEGIN_property
Au voisinage de l'état stationnaire, l'écart du produit $v = \log(\hat y/\hat y^{\star})$
suit, à l'ordre deux, l'équation de Bernoulli\nbsp{}:
\[
\dot v \approx -\beta\,v + \gamma_y\,v^2
\qquad\text{avec}\qquad
\gamma_y = \frac{\beta}{2\alpha^{\star}}\left[\alpha^{\star} + \frac{1-2\alpha^{\star}}{\sigma^{\star}}\right]
\]
dont la solution est $v(t) = v(0)e^{-\beta t}\big/\bigl[1-\frac{\gamma_y}{\beta}v(0)(1-e^{-\beta t})\bigr]$.
La vitesse de convergence du produit impliquée par l'approximation est\nbsp{}:
\[
-\frac{\ddot v}{\dot v} = \beta - 2\gamma_y\,v
\]
#+END_property

#+BEGIN_proof
Notons $c = \frac12\alpha^{\star}(1-\alpha^{\star})(1-1/\sigma^{\star})$, de sorte que
$v \approx \alpha^{\star}u + cu^2$ et $\dot v \approx (\alpha^{\star} + 2cu)\dot u$. En
substituant $\dot u \approx -\beta u + \gamma u^2$ et en ne conservant que les
termes d'ordre au plus deux, $\dot v \approx -\alpha^{\star}\beta\,u + (\alpha^{\star}\gamma - 2c\beta)u^2$.
Il reste à exprimer $u$ en fonction de $v$ au même ordre, en inversant le
développement\nbsp{}: $u \approx v/\alpha^{\star} - (c/\alpha^{\star 3})v^2$. Il vient
$\dot v \approx -\beta v + \bigl[\gamma/\alpha^{\star} - \beta c/\alpha^{\star 2}\bigr]v^2$
et, en remplaçant $\gamma$ et $c$ par leurs expressions\nbsp{}:
\[
\frac{\gamma}{\alpha^{\star}} - \frac{\beta c}{\alpha^{\star 2}}
= \frac{\beta}{2\alpha^{\star}}\left[\left(1-\frac{\alpha^{\star}}{\sigma^{\star}}\right) - (1-\alpha^{\star})\left(1-\frac{1}{\sigma^{\star}}\right)\right]
= \frac{\beta}{2\alpha^{\star}}\left[\alpha^{\star} + \frac{1-2\alpha^{\star}}{\sigma^{\star}}\right]
\]
L'équation obtenue est celle du capital avec $\gamma_y$ à la place de
$\gamma$, d'où la solution. Enfin $\ddot v = (-\beta + 2\gamma_y v)\dot v$.
#+END_proof

Deux vérifications de cohérence. D'abord, comme $v \approx \alpha^{\star}u$, la
pente de cette vitesse par rapport à $u$ à l'état stationnaire est
\(-2\gamma_y\alpha^{\star} = -(n+g+\delta)(1-\alpha^{\star})\bigl[\alpha^{\star} + (1-2\alpha^{\star})/\sigma^{\star}\bigr]\)\nbsp{}:
c'est exactement la dérivée établie dans la section sans approximation. La
vitesse de l'ordre deux est la tangente de la vitesse exacte $\beta_{\sigma}$ à
l'état stationnaire, là où l'ordre un n'en retient que la valeur $\beta$.
Ensuite, dans le cas Cobb-Douglas, $\gamma_y = \beta(1-\alpha)/(2\alpha)$, et la
vitesse exacte $\beta_1 = \beta\,e^{-\frac{1-\alpha}{\alpha}v}$ a bien pour
développement $\beta - 2\gamma_y v$ à l'ordre un en $v$.

#+begin_src python :session edo1 :exports code :results none
  def gamma_y(sigma):
      return betastar/(2*alphastar)*(alphastar + (1 - 2*alphastar)/sigma)

  def ordre2_produit(t, v0, sigma):
      gy = gamma_y(sigma)
      return v0*np.exp(-betastar*t)/(1 - (gy/betastar)*v0*(1 - np.exp(-betastar*t)))
#+end_src

#+begin_src python :session edo1 :exports none :results none
  tg = np.linspace(0, 60, 400)
  fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
  for ax, ratio in ((ax1, 0.25), (ax2, 4.0)):
      for sigma, couleur in ((0.5, 'r'), (1.0, 'b'), (2.0, 'g')):
          f, fprime, alpha, kstar = technologie(sigma)
          k = transition(sigma, ratio*kstar).sol(tg)[0]
          v0 = np.log(f(ratio*kstar)/f(kstar))
          ax.plot(tg, vitesse(sigma, k), couleur, label=r'$\sigma = %s$' % str(sigma).replace('.', '{,}'))
          ax.plot(tg, betastar - 2*gamma_y(sigma)*ordre2_produit(tg, v0, sigma), couleur,
                  linewidth=1, linestyle='--')
      ax.axhline(y=betastar, color='k', linewidth=0.8, linestyle=':')
      ax.annotate(r'$\beta^\star$', xy=(57, betastar+0.004))
      ax.set_xlabel(r'$t$ (années)')
      ax.set_ylabel(r'$\beta(t)$')
      ax.legend()
  ax1.set_title(r'départ à $0{,}25\,\hat k^\star$')
  ax2.set_title(r'départ à $4\,\hat k^\star$')
  fig.tight_layout()
  fig.savefig("solow-vitesse-ordre2.svg", transparent=True)
#+end_src

#+CAPTION: *Vitesse de convergence du produit par tête efficace\nbsp{}: dynamique exacte (trait plein, comme dans la figure de la section sans approximation) et approximation à l'ordre deux, $\beta - 2\gamma_y v(t)$ évaluée le long de la solution approchée (tirets). L'ordre un donne la constante $\beta^{\star}$ (pointillés).*
#+LABEL: fig:vitesse-ordre2
[[file:solow-vitesse-ordre2.svg]]

La figure [[fig:vitesse-ordre2][ci-dessus]] superpose les deux vitesses le long des transitions de la
section sans approximation. Près de l'état stationnaire, l'ordre deux
reproduit le sens et l'ordre de grandeur des écarts à $\beta^{\star}$, dans les
deux directions et pour les trois technologies, ce que l'ordre un ne peut pas
faire. Loin de l'état stationnaire, l'approximation, qui n'est que la
tangente de la vitesse exacte, s'en écarte. Pour $\sigma=1$ et $\sigma=2$
elle sous-estime la vitesse dans les deux directions, la vitesse exacte
étant convexe en $v$ et la tangente au-dessous d'elle. Pour $\sigma=0{,}5$
elle exagère fortement la vitesse par dessous, là où la vitesse exacte
s'aplatit parce que la fonction de transition change de courbure le long de
la trajectoire, comme on le verra plus bas\nbsp{}; et elle ne peut pas
reproduire la vitesse négative observée par dessus, qui est un effet de la
variation de la part du capital loin de son niveau stationnaire.

Trois remarques sur cette solution.

#+BEGIN_remarque
Pour $\gamma=0$ on retrouve la solution de l'ordre un, et de façon générale
le dénominateur corrige la solution de l'ordre un d'un facteur qui dépend du
signe de $u(0)$. Si $\gamma>0$ et que l'économie part au-dessous de son état
stationnaire ($u(0) < 0$), le dénominateur est supérieur à un et $|u(t)|$ est
plus petit qu'à l'ordre un\nbsp{}: la convergence est plus rapide. Si elle part
au-dessus, le dénominateur est inférieur à un et la convergence est plus
lente. C'est exactement ce que dit la figure [[fig:vitesse][précédente]] pour la dynamique
exacte. On le lit aussi sur la vitesse d'ajustement du capital. Puisque
$\dot u$ est le taux de croissance du capital par tête efficace, $g_{\hat k}$,
le rapport $-\ddot u/\dot u = -\dot g_{\hat k}/g_{\hat k}$ est l'analogue, pour
le capital, de la vitesse de convergence définie pour le produit dans la
section [[#sans-approximation][sans approximation]]\nbsp{}: il mesure la vitesse à laquelle le taux
de croissance du capital s'éteint, et vaut exactement $\beta$ pour la
dynamique linéaire $\dot u = -\beta u$. Comme $\dot u = G(u)$, on a
$\ddot u = G'(u)\dot u$ et donc \(-\ddot u/\dot u = -G'(u)\)\nbsp{}: c'est la pente
de la fonction de transition au point courant, la vitesse que donnerait une
linéarisation faite en $u$ plutôt qu'en zéro. Pour l'approximation d'ordre
deux, cette vitesse est $\beta - 2\gamma u$, qui n'est plus constante\nbsp{}:
supérieure à $\beta$ au-dessous de l'état stationnaire, inférieure au-dessus.
Pour la dynamique exacte, la dérivée de $G$ calculée plus haut donne\nbsp{}:
\[
-\frac{\ddot u}{\dot u} = (n+g+\delta)\bigl[1-\alpha(\hat k)\bigr]\frac{\hat y/\hat y^{\star}}{\hat k/\hat k^{\star}}
\]
dont $\beta-2\gamma u$ est le développement à l'ordre un en $u$. Cette
vitesse du capital n'est pas la vitesse du produit $\beta(t)$ de la propriété
établie dans la section sans approximation\nbsp{}: comme
$g_{\hat y} = \alpha(\hat k)\,g_{\hat k}$, celle-ci contient en plus le terme dû
à la variation de la part du capital, qui est à l'origine du crochet. Les deux
coïncident lorsque $\alpha$ est constante, c'est-à-dire dans le cas
Cobb-Douglas. La propriété précédente donne la contrepartie de $\beta - 2\gamma u$
pour le produit, $\beta - 2\gamma_y v$, qui se compare directement à $\beta_{\sigma}$.
#+END_remarque

#+BEGIN_remarque
L'équation approchée possède un second état stationnaire, $u = \beta/\gamma$,
instable, qui n'a pas de contrepartie dans le modèle\nbsp{}: c'est un artefact de
l'approximation, comme la figure sur l'[[#locale][approximation locale]] le laissait
prévoir pour toute approximation polynomiale de la fonction de transition.
Si $u(0)>\beta/\gamma$, le dénominateur s'annule en un temps fini, et la
solution approchée explose exactement comme celle de l'[[ex-bernoulli][exemple]] traité plus
haut. Pour une technologie Cobb-Douglas, $\beta/\gamma = 2/(1-\alpha)$, soit
$3{,}125$ pour \(\alpha=0{,}36\)\nbsp{}: il faudrait partir de plus de vingt-deux
fois le capital stationnaire, ce qui est très loin du voisinage où
l'approximation a un sens. L'approximation à l'ordre deux s'utilise donc sans
crainte, mais elle ne s'utilise pas n'importe où.
#+END_remarque

#+BEGIN_remarque
Dans le cas Cobb-Douglas, le modèle exact est lui-même une équation de
Bernoulli\nbsp{}: $\dot{\hat k} = s\hat k^{\alpha} - (n+g+\delta)\hat k$, avec
$m=\alpha$. Le changement de variable $z = \hat k^{1-\alpha}$ la rend linéaire,
$\dot z = (1-\alpha)\left[s - (n+g+\delta)z\right]$, et l'on obtient la solution
explicite de la dynamique de transition, utilisée dans la note sur la
[[https://stephane-adjemian.fr/posts/simulation-du-modele-de-solow/][simulation du modèle de Solow]]\nbsp{}:
\[
u(t) = \frac{1}{1-\alpha}\log\left[1 + \left(e^{(1-\alpha)u(0)} - 1\right)e^{-\beta t}\right]
\]
De façon équivalente, la variable $w = e^{-(1-\alpha)u} = (\hat k/\hat k^{\star})^{-(1-\alpha)}$
suit la dynamique logistique $\dot w = \beta w(1-w)$. L'approximation à
l'ordre deux remplace donc une équation de Bernoulli en $\hat k$ par une
équation de Bernoulli en $u$, et l'approximation à l'ordre un par une
équation linéaire en $u$. C'est la raison pour laquelle cette technologie
permet de mesurer exactement les erreurs d'approximation.
#+END_remarque

La figure [[fig:ordre2][ci-dessous]] compare la dynamique exacte de l'écart $u$, obtenue par
la formule précédente dans le cas Cobb-Douglas et par intégration numérique
dans le cas $\sigma=0{,}5$, aux approximations d'ordre un et d'ordre deux, en
partant de $0{,}25\,\hat k^{\star}$ et de $4\,\hat k^{\star}$.

#+begin_src python :session edo1 :exports code :results none
  def ordre1(t, u0):
      return u0*np.exp(-betastar*t)

  def ordre2(t, u0, sigma):
      gamma = betastar/2*(1 - alphastar/sigma)
      return u0*np.exp(-betastar*t)/(1 - (gamma/betastar)*u0*(1 - np.exp(-betastar*t)))
#+end_src

#+begin_src python :session edo1 :exports none :results none
  tg = np.linspace(0, 60, 400)
  fig, axes = plt.subplots(2, 2, figsize=(11, 8))
  erreurs = {}
  for j, sigma in enumerate((1.0, 0.5)):
      kstar = technologie(sigma)[3]
      for ratio, style in ((0.25, '-'), (4.0, '--')):
          u0 = np.log(ratio)
          if sigma == 1.0:
              u = np.log(1 + (np.exp((1-alphastar)*u0) - 1)*np.exp(-betastar*tg))/(1-alphastar)
          else:
              u = np.log(transition(sigma, ratio*kstar).sol(tg)[0]/kstar)
          u1, u2 = ordre1(tg, u0), ordre2(tg, u0, sigma)
          erreurs[(sigma, ratio)] = (np.abs(u1-u).max(), np.abs(u2-u).max())
          axes[0, j].plot(tg, u, 'k', linestyle=style, linewidth=1.8, label='exacte' if ratio < 1 else None)
          axes[0, j].plot(tg, u1, 'b', linestyle=style, linewidth=1, label='ordre un' if ratio < 1 else None)
          axes[0, j].plot(tg, u2, 'r', linestyle=style, linewidth=1, label='ordre deux' if ratio < 1 else None)
          axes[1, j].plot(tg, u1-u, 'b', linestyle=style, linewidth=1,
                          label=r'ordre un, départ à $%s\,\hat k^\star$' % ('0{,}25' if ratio < 1 else '4'))
          axes[1, j].plot(tg, u2-u, 'r', linestyle=style, linewidth=1,
                          label=r'ordre deux, départ à $%s\,\hat k^\star$' % ('0{,}25' if ratio < 1 else '4'))
      axes[0, j].axhline(y=0, color='k', linewidth=0.8, linestyle=':')
      axes[0, j].set_title(r'$\sigma = %s$' % str(sigma).replace('.', '{,}'))
      axes[0, j].set_ylabel(r'$u(t) = \log(\hat k/\hat k^\star)$')
      axes[0, j].legend()
      axes[1, j].axhline(y=0, color='k', linewidth=0.8, linestyle=':')
      axes[1, j].set_ylabel('erreur')
      axes[1, j].legend(fontsize=8)
      for ax in axes[:, j]:
          ax.set_xlabel(r'$t$ (années)')
  fig.tight_layout()
  fig.savefig("solow-ordre2.svg", transparent=True)
#+end_src

#+CAPTION: *Dynamique exacte de l'écart à l'état stationnaire (noir) et approximations d'ordre un (bleu) et d'ordre deux (rouge), en partant de $0{,}25\,\hat k^{\star}$ (trait plein) et de $4\,\hat k^{\star}$ (tirets), pour une technologie Cobb-Douglas (à gauche) et une CES avec $\sigma=0{,}5$ (à droite). En bas, les erreurs d'approximation.*
#+LABEL: fig:ordre2
[[file:solow-ordre2.svg]]

L'erreur maximale sur l'écart logarithmique, en valeur absolue, est la
suivante\nbsp{}:

| départ                   | technologie          | ordre un  | ordre deux |
|--------------------------+----------------------+-----------+------------|
| $0{,}25\,\hat k^{\star}$ | Cobb-Douglas         | $0{,}152$ | $0{,}027$  |
| $4\,\hat k^{\star}$      | Cobb-Douglas         | $0{,}152$ | $0{,}050$  |
| $0{,}25\,\hat k^{\star}$ | CES, $\sigma=0{,}5$  | $0{,}029$ | $0{,}036$  |
| $4\,\hat k^{\star}$      | CES, $\sigma=0{,}5$  | $0{,}101$ | $0{,}026$  |

Dans le cas Cobb-Douglas, l'ordre deux divise l'erreur de l'ordre un par
plus de cinq par dessous et par trois par dessus, alors que les départs
retenus sont loin d'être petits\nbsp{}: un capital initial quatre fois
inférieur ou supérieur au capital stationnaire. L'approximation à l'ordre un
se trompe toujours dans le même sens, elle sous-estime la vitesse par
dessous et la surestime par dessus, ce qui est exactement le contenu de la
propriété établie dans la section [[#sans-approximation][sans approximation]]\nbsp{}; l'ordre deux,
en rétablissant la dépendance de la vitesse à la position, corrige cette
asymétrie. Il reste une erreur, plus importante par dessus que par dessous,
qui ne se réduirait qu'en poussant le développement à l'ordre suivant.

Le cas $\sigma=0{,}5$ réserve une surprise\nbsp{}: par dessous, l'ordre deux
fait /moins bien/ que l'ordre un. Ce n'est pas une erreur de calcul, mais
un rappel de ce qu'est une approximation locale. La courbure de la fonction
de transition est, d'après la preuve de la propriété, du signe de
$1-\alpha(\hat k)/\sigma$. Avec $\sigma=0{,}5$, la part du capital augmente
lorsque le capital diminue, et dépasse $\sigma$ dès que $\hat k$ passe
au-dessous de la moitié de \(\hat k^{\star}\)\nbsp{}: la fonction $G$ est convexe
au voisinage de l'état stationnaire, où l'approximation est construite, mais
concave sur la première partie de la trajectoire, où elle est utilisée. Le
terme quadratique, qui reproduit la courbure à l'état stationnaire, corrige
alors dans le mauvais sens. Que l'ordre un s'en tire bien tient à une
compensation, la vitesse exacte du capital étant tantôt inférieure, tantôt
supérieure à $\beta^{\star}$ le long de cette trajectoire, ce que la figure sur
la [[fig:vitesse][vitesse de convergence]] laissait entrevoir. Par dessus, où la fonction
reste convexe, l'ordre deux divise l'erreur par quatre.

Il ne faut pas en conclure que l'approximation à l'ordre deux est toujours
préférable. Elle demande de connaître la courbure de la fonction de
transition à l'état stationnaire, ici l'élasticité de substitution, là où
l'ordre un se contente de la part du capital\nbsp{}; et sa solution explicite
est un privilège de la dimension un, qui disparaît dès que l'on passe à un
système. C'est ce qui explique que la log-linéarisation reste la méthode de
référence pour les modèles d'équilibre général dynamiques, et que l'on
n'aille chercher les ordres supérieurs que lorsque la question posée
l'exige.