Équations différentielles du premier ordre

Nous commençons par le cas le plus simple, où la constante \(b\) est nulle :

\[ \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 :

\[ \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 :

\[ \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 :

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 :

\[ 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\).

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

\[ \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 précédente ; 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 :

\[ \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 :

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

où \(\Gamma = e^{\gamma}\). Cette équation définit un continuum de fonctions indexées par la constante \(\Gamma\) ; 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 :

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

Il s'agit de l'état stationnaire de l'équation complète : 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 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 :

\[ 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 :

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

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

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

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

que l'on supposera strictement positif, ce qui exige \(\alpha>\gamma\) : à 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 :

\[ \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 :

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

En appliquant la formule de la section précédente, on obtient directement :

\[ 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 ; 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.

On s'intéresse à un problème de la forme :

\[ \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 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 :

\[ \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 :

\[ 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 :

• la variable \(y\) converge vers l'état stationnaire nul si et seulement si \(A(t)\) tend vers \(+\infty\) ;
• elle converge vers une limite finie non nulle si et seulement si \(A(t)\) admet une limite finie ;
• 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 ; 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.

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 :

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

De façon générale, l'équation différentielle :

\[ 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 :

\[ \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¹.

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 :

\[ 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 :

\[ 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\) ; déterminer la fonction \(\psi(t)\) est souvent moins évident.

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.

On s'intéresse à un problème de la forme :

\[ \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 :

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

Nous avons \(M=1\) et \(N = ay-b\) ; 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 :

\[ \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 :

\[ \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 :

\[ \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\).

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 : chaque membre s'intègre séparément.

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\) :

\[ 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 :

\[ \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.

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 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 :

\[ \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 ci-dessous donne deux exemples : 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 :

1. lorsque \(f(y)>0\), la variation \(\dot y\) est positive et donc \(y\) augmente ;
2. lorsque \(f(y) < 0\), la variation \(\dot y\) est négative et donc \(y\) diminue ;
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 ; \(y\) diminue lorsque son niveau est élevé, et les flèches sont dirigées vers la gauche.

\(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 : \(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 ci-dessous montre que la fonction \(f\) peut être beaucoup plus « tordue ». 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}\).

On retient deux idées de ces exemples :

1. l'état stationnaire, s'il existe, n'est pas toujours unique. On peut aussi imaginer des cas où il n'existe pas ;
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\) : 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\).

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 ; 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 :

\[ \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'exemple traité plus haut ; comme nous disposons d'une solution analytique, nous pourrons évaluer la précision de la solution numérique. Pour rappel, la solution analytique est :

\[ 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 ; 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 : \[ t^{\star} = \sqrt{2\log\frac{3}{3-\frac{1}{y(0)}}} \] La solution admet une asymptote verticale en \(t^{\star}\) : 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 :

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)

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.

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)

La figure 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 : 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 : 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}\) ; et il n'est pas garanti qu'il signale son échec autrement que par des valeurs aberrantes.

Comme dans la section sur l'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 :

\[ \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 :

\[ \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 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 :

\[ \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 ci-dessous illustre ce remplacement, pour la fonction de transition du modèle de Solow que nous retrouverons dans la section suivante.

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 :

\[ \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.

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 le modèle de Solow ; le but est simplement d'appliquer les méthodes présentées plus haut.

On se donne une fonction de production néoclassique², \(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³. La loi d'évolution du stock de capital physique est :

\[ \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 :

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

Cette équation différentielle est non autonome : 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, « purgé » 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 :

\[ \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 :

\[ \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\) : l'état stationnaire est unique, et il vérifie :

\[ \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 :

\[ \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 ici, qui s'écrit en fonction de la technologie intensive :

\[ \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 :

\[ \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 :

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é ;
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 ;
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\) :

\[ 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 : \(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 : \(\gamma_1 = \alpha^{\star}\left((n+g+\delta)/s\right)^{\rho}\) et \(\gamma_2 = 1-\gamma_1\).

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)