#+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-en.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: en-GB
#+STARTUP: latexpreview
#+TITLE: First-order differential equations
#+DATE: September 2026
#+AUTHOR: Stéphane Adjemian
#+EMAIL: stephane.adjemian@univ-lemans.fr
#+PROPERTY: header-args:python :python /tmp/blog-edo1/bin/python

#+BEGIN_QUOTE
A first-order differential equation is a functional equation in which only
the unknown function and its first derivative appear. We start with linear
problems, for which an analytical treatment is always possible, before
turning to nonlinear equations. In the latter case a solution cannot always
be exhibited, and one has to resort to a trick (a change of variable may make
the problem linear), to an approximation (a Taylor expansion), to graphical
analysis (to derive properties of the solution without computing it) or to
the computer. The numerical computations are done in Python, and the
solutions of the exercises are given in collapsible blocks. The note ends
with an application to the Solow model, where the explicit solution of the
second-order approximation of the transition dynamics is obtained.
#+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-en :exports none :results none
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.integrate import solve_ivp
#+end_src

* Linear equations with constant coefficients
:PROPERTIES:
:CUSTOM_ID: constants
:END:

In this section, we consider a differential equation of the form:

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

where $a$ and $b$ are real parameters and $\dot y$ denotes the derivative of
$y$ with respect to time. We distinguish two cases according to whether $b$
is zero or not.

** The homogeneous case
:PROPERTIES:
:CUSTOM_ID: homogene
:END:

We start with the simplest case, where the constant $b$ is zero:

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

A differential equation is said to be homogeneous when the dynamics is not
affected if $y(t)$ and $\dot y(t)$ are multiplied by the same constant. This
equation can be written equivalently as:

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

The dynamics therefore corresponds to an assumption of growth at a constant
rate. To solve the differential equation, that is to exhibit an expression
of $y$ in terms of the parameter $a$ and of $t$, it suffices to know how to
differentiate the logarithm. Indeed, since $\frac{\mathrm d}{\mathrm dt}\log u(t) =
\dot u(t)/u(t)$, the previous equation also reads:

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

This equation tells us that the changes in the logarithm of $y$ are
constant. To obtain $\log y(T)$ we sum these changes between $0$ and $T$,
that is, we integrate both sides:

\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*}

where the last equivalence is obtained by applying the inverse of the natural
logarithm, the exponential. This transformation takes us back to the
variable of interest, $y(t)$ rather than $\log y(t)$. Finally:

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

This is the solution of the homogeneous first-order linear differential
equation with constant coefficient. It depends on the initial condition, on
the variable $t$ and on the parameter $a$.

#+BEGIN_remarque
If the initial condition is zero, then $y(t)=0$ for all $t$ in $\mathbb R^+$.
Zero is said to be the /steady state/, or the /fixed point/, of this
dynamics.
#+END_remarque

#+BEGIN_remarque
The properties of the dynamics are tied to the sign of the parameter $a$. If
$a>0$ then, for any initial condition $y(0)$, \(\lim_{t\to\infty}y(t)=0\): the
dynamics is /stable/, in the long run the variable $y$ reaches the steady
state, and the long-run level does not depend on the initial condition.
Conversely, if $a < 0$ then, for any initial condition $y(0)\neq0$,
\(\lim_{t\to\infty}|y(t)|=\infty\): the variable /diverges/, towards $-\infty$ or
$+\infty$ depending on the sign of the initial condition. The figure
[[fig:homogene][below]] illustrates these cases.
#+END_remarque

#+begin_src python :session edo1-en :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$: stable dynamics')
  ax2.set_title(r'$a = -0.5 < 0$: unstable dynamics')
  fig.tight_layout()
  fig.savefig("edo-homogene.svg", transparent=True)
#+end_src

#+CAPTION: *Solutions of the homogeneous equation for several initial conditions. On the left the paths converge to the steady state, on the right they move away from it.*
#+LABEL: fig:homogene
[[file:edo-homogene.svg]]

** The non-homogeneous case
:PROPERTIES:
:CUSTOM_ID: non-homogene
:END:

We now turn to the case $b\neq0$, that is to a problem of the form:

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

The solution is obtained in four steps.

*Step 1. The general solution of the associated homogeneous equation.* The
associated equation without right-hand side is $\dot y(t) = -a\,y(t)$. We
already met it in the [[#homogene][previous]] section; we shall look for its general
solution following a slightly different route. Assuming $y\neq0$, we have
equivalently $\frac{\mathrm d}{\mathrm dt}\log y(t) = -a$. The derivative with
respect to $t$ of the logarithm of $y(t)$ equals $-a$ for all $t$, so the
antiderivatives of both sides coincide up to a constant:

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

or $\log y(t) = -at + \gamma$, since we will not be able to identify the two
constants $\gamma_1$ and $\gamma_2$ separately (we shall understand why
later). Applying the exponential, we finally get:

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

where $\Gamma = e^{\gamma}$. This equation defines a continuum of functions
indexed by the constant \(\Gamma\); each is a solution of the equation without
right-hand side. It is in this sense that it is a /general/ solution. At this
stage, the constant $\Gamma$ is not determined.

*Step 2. A particular solution of the complete equation.* We choose the
simplest possible one, the constant solution:

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

It is the steady state of the complete equation: if $y(t)=b/a$ then
$y(t+\Delta)=b/a$ for all $\Delta>0$. Note in passing that we have implicitly
assumed $a\neq0$. If this assumption fails, a constant function is not a
particular solution and another one must be sought (see the exercise
[[ex-a-nul][below]]).

*Step 3. The general solution of the complete equation.* Adding the general
solution of the equation without right-hand side and a particular solution
of the complete equation, we obtain the general solution of the complete
equation:

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

This solution is general in that the parameter $\Gamma$ is still not
determined. It fixes the shape of the solution function, but not its level.
Choosing a value for $\Gamma$ amounts to choosing one path among infinitely
many possible paths.

*Step 4. Solution of the complete equation.* Choosing $\Gamma$ is easy if we
know the value of $y$ at some instant. For instance, we may know the initial
condition $y(0)=y_0$ (one also speaks of a boundary condition). At instant
$0$ we have $y(0) = \Gamma + b/a$, that is $\Gamma = y_0 - b/a$. Hence the
solution of the non-homogeneous differential equation is:

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

#+BEGIN_remarque
This solution is only valid if $a\neq0$. Otherwise, the solution is obtained
much more simply. We would indeed have $\dot y(t) = b$ and, integrating with
respect to $t$, $y(t) = \gamma + bt$ where $\gamma$ is an arbitrary constant
fixed by an initial condition. The four-step approach is not recommended in
that case, even though it remains feasible (exercise [[ex-a-nul][below]]).
#+END_remarque

#+BEGIN_remarque
The solution receives the same interpretation as in the homogeneous case.
The dynamics is stable if and only if the parameter $a$ is positive. In that
case, for any initial condition $y(0)$, we have $\lim_{t\to\infty}y(t)=b/a$.
If the parameter $a$ is negative, the variable $y(t)$ diverges towards
$+\infty$ or $-\infty$ depending on the sign of $y(0)-b/a$, as long as
$y(0)\neq b/a$. If the initial condition equals the steady state, then
$y(t)=b/a$ for all $t$, whatever the sign of $a$. The solution can be
written in the following form:

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

The level of $y$ at instant $t$ is a convex combination of the initial
condition and of the steady state. If $a>0$, the influence of the initial
condition tends to vanish while the weight of the steady state (the target)
tends to one.
#+END_remarque

#+NAME: ex-a-nul
#+BEGIN_exercice
Show that, even in the case $a=0$, it is possible to solve the
non-homogeneous differential equation $\dot y(t) + a\,y(t) = b$ following the
four-step approach described above. Hint: look for a /non-constant/ solution
of the complete equation.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
With $a=0$ the equation reads $\dot y(t) = b$. /Step 1/: the equation without
right-hand side is $\dot y(t)=0$, whose general solution is the constant
function $y_1(t)=\Gamma$. /Step 2/: a constant function cannot be a solution
of the complete equation since its derivative is zero and not equal to
\(b\); we therefore look for a particular solution linear in $t$,
$y_2(t)=ct$, and the equation imposes $c=b$. /Step 3/: the general solution
is $y(t) = \Gamma + bt$. /Step 4/: the initial condition gives $\Gamma=y_0$,
hence $y(t) = y_0 + bt$. We recover the result of direct integration. Note
that this solution cannot be obtained as the limit of
$\left(y_0 - b/a\right)e^{-at} + b/a$ as $a$ tends to zero term by term,
each of the two terms diverging; but the form
$e^{-at}y_0 + \frac{1-e^{-at}}{a}\,b$ does have $y_0+bt$ as its limit, since
$\frac{1-e^{-at}}{a}\to t$.
#+END_details

** Dynamics of a market price
:PROPERTIES:
:CUSTOM_ID: prix
:END:

Suppose the demand and supply functions of a good are given by:

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

where $P$ is the price of the good. The parameters $\beta$ and $\delta$
measure the sensitivity of demand and supply to a change in the price (they
are not elasticities). Equating the quantities supplied and demanded, we
obtain the market-clearing price:

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

which we assume strictly positive, which requires \(\alpha>\gamma\): at a zero
price, demand exceeds supply. A priori, the actual price $P$ differs from
$P^{\star}$. We shall show that if the price rises when demand exceeds
supply, and falls when supply exceeds demand, then the actual price converges
to the market-clearing price. We formalise this assumption by supposing
that:

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

with $j>0$. The price is invariant if and only if supply equals demand.
Substituting the supply and demand functions, we obtain the following
differential equation for the price of the good:

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

Applying the formula of the previous section, we directly get:

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

with $\kappa \equiv j(\beta+\delta)$. Since $\kappa>0$, the price converges in
the long run to $P^{\star}$. Convergence is monotonically increasing if
$P(0)$ is below $P^{\star}$, that is if there is initially excess demand, and
monotonically decreasing if $P(0)$ is above $P^{\star}$, that is if there is
initially excess supply. $P^{\star}$ is a particular solution of the
differential equation; it represents the intertemporal equilibrium level of
the variable of interest. The term $[P(0)-P^{\star}]e^{-\kappa t}$ accounts
for the deviations from this equilibrium level.

#+begin_src python :session edo1-en :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, colour, name in ((1.0, 'b', "initial excess demand"), (7.0, 'g', "initial excess supply")):
      ax2.plot(t, (P0-Pstar)*np.exp(-kappa*t)+Pstar, colour, label=name)
  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: *On the left, supply, demand and the equilibrium price ($\alpha=10$, $\beta=\delta=1$, $\gamma=2$). On the right, price adjustment from a situation of excess demand and from a situation of excess supply ($j=0.5$).*
#+LABEL: fig:prix
[[file:edo-prix.svg]]

#+BEGIN_exercice
Solve the following differential equations:
- (i) $\dot y(t) + 4y(t) = 12$ with $y(0)=2$.
- (ii) $\dot y(t) - 2y(t) = 0$ with $y(0)=9$.
- (iii) $\dot y(t) + 10y(t) = 15$ with $y(0)=0$.
- (iv) $2\dot y(t) + 4y(t) = 6$ with $y(0)=3/2$.
- (v) $\dot y(t) + y(t) = 4$ with $y(0)=0$.
- (vi) $\dot y(t) = 23$ with $y(0)=1$.
- (vii) $3\dot y(t) + 6y(t) = 5$ with $y(0)=0$.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
We apply the formula $y(t) = (y_0 - b/a)e^{-at} + b/a$, after dividing the
equation, where necessary, by the coefficient of $\dot y$.
- (i) $a=4$, $b=12$, steady state \(3\): $y(t) = 3 - e^{-4t}$.
- (ii) Homogeneous equation with \(a=-2\): $y(t) = 9e^{2t}$, the dynamics is
  unstable.
- (iii) $a=10$, $b=15$, steady state \(3/2\):
  $y(t) = \frac32\left(1-e^{-10t}\right)$.
- (iv) Dividing by $2$, $\dot y + 2y = 3$, whose steady state is $3/2$. The
  initial condition equals the steady state: $y(t)=3/2$ for all $t$.
- (v) $a=1$, \(b=4\): $y(t) = 4\left(1-e^{-t}\right)$.
- (vi) \(a=0\): we integrate directly, $y(t) = 1 + 23t$.
- (vii) Dividing by $3$, $\dot y + 2y = 5/3$, whose steady state is \(5/6\):
  $y(t) = \frac56\left(1-e^{-2t}\right)$.
#+END_details

* Linear equations with variable coefficients
:PROPERTIES:
:CUSTOM_ID: variables
:END:

** The homogeneous case
:PROPERTIES:
:CUSTOM_ID: homogene-variable
:END:

We consider a problem of the form:

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

where $a(t)$ is a continuous real function. We can solve this equation
following the same route as in the case of [[#homogene][constant coefficients]]. The growth
rate of $y$ at instant $t$ is $\dot y(t)/y(t) = -a(t)$, or, using the
properties of the logarithm:

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

Integrating both sides between $0$ and $t$, the initial condition being
assumed known, we get $\log y(t) - \log y(0) = -\int_0^t a(\tau)\,\mathrm d\tau$,
or equivalently:

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

One checks that if $a(t)=a$ for all $t$, then $\int_0^t a(\tau)\,\mathrm d\tau = at$
and we recover the solution obtained with constant coefficients. The
stability conditions are less trivial here, since everything depends on the
shape of the function $a(t)$. Writing $A(t) = \int_0^t a(\tau)\,\mathrm d\tau$,
we have $|y(t)| = |y(0)|e^{-A(t)}$ and three cases arise:
- the variable $y$ converges to the zero steady state if and only if $A(t)$
  tends to \(+\infty\);
- it converges to a finite non-zero limit if and only if $A(t)$ has a finite
  limit;
- it remains bounded, without necessarily converging, if and only if $A(t)$
  is bounded below.
The dynamics diverges in the remaining case, that is if $A(t)$ is not bounded
below. A function $a$ bounded below by a strictly positive constant
guarantees convergence to zero, as in the constant case. But
$a(t) = (1+t)^{-2}$, though strictly positive, gives $A(t) = 1-(1+t)^{-1}$ and
$y(t)$ converges to $y(0)e^{-1}$ without ever reaching the steady state; and
$a(t) = \cos t$ gives $A(t) = \sin t$, so that $y(t)$ oscillates indefinitely
between $y(0)e^{-1}$ and $y(0)e$. The sign of $a$ at a given instant only
tells us the direction of change of $y$ at that instant.

** Exact differential equations
:PROPERTIES:
:CUSTOM_ID: exactes
:END:

We leave the world of linear differential equations for a moment to
introduce a tool which will allow us to build the solution of a
non-homogeneous linear differential equation with variable coefficients.

Let $F(y,t)$ be a continuously differentiable function of two variables. Its
total differential is:

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

#+BEGIN_definition
The equation
\[
\frac{\partial F}{\partial y}\,\mathrm dy + \frac{\partial F}{\partial t}\,\mathrm dt = 0
\]
is an /exact differential equation/, because the left-hand side is exactly
the differential of $F(y,t)$.
#+END_definition

#+BEGIN_exemple
Let $F(y,t) = y^2t + k$, where $k$ is a constant. We have
$\mathrm dF = 2yt\,\mathrm dy + y^2\,\mathrm dt$. Thus $2yt\,\mathrm dy + y^2\,\mathrm dt = 0$, or
equivalently $\dot y + \frac{y^2}{2yt} = 0$, is an exact differential
equation.
#+END_exemple

In general, the differential equation:

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

is exact if and only if there exists a function $F(y,t)$ such that
$M = \partial F/\partial y$ and $N = \partial F/\partial t$. We lack a test to
know whether a differential equation is exact, that is, whether there exists
a function $F$ from which the functions $M$ and $N$ derive. We know that a
Hessian matrix is symmetric (Young's theorem), that is
$\frac{\partial^2F}{\partial t\partial y} = \frac{\partial^2F}{\partial y\partial t}$.
Hence the differential equation is exact if and only if:

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

This condition gives us a test to assess whether a differential equation is
exact[fn:: The condition is necessary. It is sufficient on a simply
connected domain, which we shall always assume.].

#+BEGIN_exemple
Applying this test to the previous example, where $M=2yt$ and $N=y^2$, we get
$\partial M/\partial t = 2y$ and $\partial N/\partial y = 2y$. The test
concludes that it is indeed an exact differential equation.
#+END_exemple

An exact differential equation, by definition, tells us that
$\mathrm dF(y,t) = 0$. Its general solution must therefore be of the form:

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

where $c$ is a real constant. Solving an exact differential equation means
exhibiting an antiderivative $F(y,t)$ and setting it equal to a constant.

*Solution method.* Since $M = \partial F/\partial y$, the function $F$ must
contain an integral of $M$ with respect to the variable $y$. We should
therefore have:

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

The partial derivative $M$ is integrated with respect to $y$ only, treating
$t$ as a constant. Since, when differentiating $F(y,t)$ partially with
respect to $y$, any additive term not depending on $y$ disappears, care must
be taken to reintroduce these terms in the integration process. This is
exactly the role of the term $\psi(t)$. It is relatively easy to evaluate
\(\int M\,\mathrm dy\); determining the function $\psi(t)$ is often less
obvious.

#+BEGIN_exemple
We want to solve the differential equation $\dot y + \frac{y}{2t} = 0$.
Multiplying both sides by $2yt\,\mathrm dt$, we get
$2yt\,\mathrm dy + y^2\,\mathrm dt = 0$. We therefore have $M=2yt$ and $N=y^2$, and
we solve this equation in four steps.
- (i) Set $F(y,t) = \int 2yt\,\mathrm dy + \psi(t) = y^2t + \psi(t)$, where
  $\psi(t)$ remains to be determined and where we have redefined $\psi(t)$
  to include the constant of integration.
- (ii) The partial derivative with respect to $t$ is $\partial F/\partial t =
  y^2 + \psi'(t)$. Comparing with $N=y^2$, we deduce that the function
  $\psi(t)$ must satisfy $\psi'(t)=0$ for all $t$.
- (iii) We therefore know that $\psi(t) = \kappa\in\mathbb R$ for all $t$.
- (iv) Finally, $F(y,t) = y^2t + \kappa$ and the solution of the exact
  differential equation must be of the form $y^2t + \kappa = c$. Since the
  constants $\kappa$ and $c$ are not individually identifiable, we still
  have $y^2t = \tilde c$, with $\tilde c\geq0$, hence finally, for \(t>0\):
  \[
  y(t) = \bar c\,t^{-\frac12}
  \]
  where $\bar c = \pm\sqrt{\tilde c}$ is a real constant, with the sign of $y$,
  which can be determined from an initial condition.
#+END_exemple

#+BEGIN_exemple
Consider the differential equation:
\[
\dot y + \frac{y+3t^2}{t+2y} = 0
\]
Equivalently, we have $(t+2y)\,\mathrm dy + (y+3t^2)\,\mathrm dt = 0$. Is it an
exact differential equation? We have $M = t+2y$ and $N = y+3t^2$, and we
check that $\partial M/\partial t = 1 = \partial N/\partial y$. We are
therefore indeed dealing with an exact differential equation, which we solve
in four steps.
- (i) Set $F(y,t) = \int(t+2y)\,\mathrm dy + \psi(t) = yt + y^2 + \psi(t)$,
  where $\psi(t)$ remains to be determined.
- (ii) The partial derivative with respect to $t$ is $\partial F/\partial t =
  y + \psi'(t)$. Comparing with $N = y+3t^2$, we deduce that $\psi'(t) = 3t^2$
  for all $t$.
- (iii) We therefore know that $\psi(t) = t^3 + \kappa$ with $\kappa\in\mathbb R$.
- (iv) Finally, $F(y,t) = yt + y^2 + t^3 + \kappa$ and the solution of the
  differential equation is of the form:
  \[
  yt + y^2 + t^3 = c
  \]
  where $c$ is a real constant. One checks that it is indeed the solution by
  differentiating this relation with respect to $t$. It defines $y$
  implicitly as a function of \(t\); here it can even be made explicit, since
  it is a quadratic equation in $y$.
#+END_exemple

The procedure described in these two examples can be applied to any exact
differential equation. In some cases, the procedure can also be applied to a
non-exact differential equation, if an equivalent exact differential
equation can be found.

#+BEGIN_exemple
Consider the differential equation $2t\,\mathrm dy + y\,\mathrm dt = 0$. One easily
checks that it is not exact: $\partial M/\partial t = 2$ and
$\partial N/\partial y = 1$. However, multiplying each term by $y$, we are
back to the previous example and therefore obtain an exact differential
equation. $y$ is said to be an /integrating factor/ of the equation.
#+END_exemple

** The non-homogeneous case
:PROPERTIES:
:CUSTOM_ID: non-homogene-variable
:END:

We consider a problem of the form:

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

where $a(t)$ and $b(t)$ are continuous real functions. This problem can be
written equivalently as:

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

We have $M=1$ and \(N = ay-b\); one immediately sees that this is not an exact
differential equation, unless $a$ is identically zero. Nevertheless it is
possible to exhibit an integrating factor, which we denote $\mathcal I$, so
as to obtain an exact differential equation. The integrating factor is such
that:

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

is an exact differential equation. For this, it is necessary and sufficient
that the condition $\partial M/\partial t = \partial N/\partial y$, with
$M = \mathcal I$ and $N = \mathcal I(ay-b)$, be satisfied. The integrating
factor is therefore such that $\dot{\mathcal I} = \mathcal I a$, or
equivalently:

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

The growth rate of the integrating factor must equal $a(t)$. The integrating
factor is therefore defined by a differential equation we know how to solve:

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

for any non-zero real value of $A$. Without loss of generality we set $A=1$.

#+BEGIN_property
The general solution of the first-order linear differential equation with
variable coefficients $\dot y(t) + a(t)\,y(t) = b(t)$ is:
\[
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)
\]
where the constant $c$ is determined by the initial condition: if $y(0)=y_0$
is known, then $c = y_0$.
#+END_property

#+BEGIN_proof
The transformed differential equation:
\[
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
\]
is exact by construction. We solve it following the four steps described
above.
- (i) Set $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)$, where $\psi(t)$ remains to
  be determined.
- (ii) The partial derivative with respect to $t$ is:
  \[
  \frac{\partial F}{\partial t} = y\,a(t)\,e^{\int_0^t a(\tau)\,\mathrm d\tau} + \psi'(t)
  \]
  Comparing with $N$ we obtain a restriction on the function \(\psi\):
  $\psi'(t) = -b(t)\,e^{\int_0^t a(\tau)\,\mathrm d\tau}$.
- (iii) The function $\psi$ is therefore, the initial instant being zero:
  \[
  \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
  \]
  We cannot go further here, because the functions $a$ and $b$ are not
  specified.
- (iv) Finally, substituting the expression of $\psi(t)$ into the postulated
  function $F(y,t)$, the solution satisfies:
  \[
  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
  \]
  which gives the stated expression. At $t=0$ both integrals vanish and we
  get $y(0)=c$.
#+END_proof

#+BEGIN_exemple
Consider the differential equation $\dot y(t) + 2t\,y(t) = t$, the initial
condition $y(0)$ being known. We have $a(t)=2t$ and $b(t)=t$, hence
$\int_0^t a(\tau)\,\mathrm d\tau = t^2$. Applying the previous result, we get:
\[
y(t) = e^{-t^2}\left(y(0) + \int_0^t s\,e^{s^2}\,\mathrm ds\right)
\]
Noting that $\frac{\mathrm d}{\mathrm ds}e^{s^2} = 2s\,e^{s^2}$, the integral equals
$\frac12\left(e^{t^2}-1\right)$ and:
\[
y(t) = e^{-t^2}\left(y(0) - \frac12\right) + \frac12
\]
Note that $1/2$ is a constant particular solution of the differential
equation: it is a steady state, stable in this example since the
exponential term tends to zero as $t$ tends to infinity. It could have been
obtained directly by looking, as in the constant-coefficient case, for a
constant particular solution, then adding the general solution
$\Gamma e^{-t^2}$ of the homogeneous equation.
#+END_exemple

#+BEGIN_exercice
Solve the following differential equations:
- (a) $\dot y + 5y = 15$.
- (b) $\dot y + 2ty = 0$.
- (c) $\dot y + 2ty = t$ with $y(0)=3/2$.
- (d) $\dot y + t^2y = 5t^2$ with $y(0)=6$.
- (e) $2\dot y + 12y + 2e^t = 0$ with $y(0)=6/7$.
- (f) $\dot y + y = t$.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
When the initial condition is not given, we settle for the general solution,
indexed by a constant $\Gamma$.
- (a) Constant coefficients, steady state \(3\): $y(t) = 3 + \Gamma e^{-5t}$,
  with $\Gamma = y(0)-3$.
- (b) Homogeneous equation with $a(t)=2t$, hence \(\int_0^t a = t^2\):
  $y(t) = \Gamma e^{-t^2}$, with $\Gamma = y(0)$.
- (c) This is the example treated above: $y(t) = \frac12 + \left(\frac32-\frac12\right)e^{-t^2}
  = \frac12 + e^{-t^2}$.
- (d) $a(t)=t^2$, $b(t)=5t^2$, hence $\int_0^t a = t^3/3$ and the integrating
  factor is $e^{t^3/3}$. One observes that $5$ is a constant particular
  solution, hence $y(t) = 5 + \Gamma e^{-t^3/3}$ and, with the initial
  condition, $y(t) = 5 + e^{-t^3/3}$. The general formula gives the same
  result, the integral $\int_0^t 5s^2e^{s^3/3}\,\mathrm ds = 5\left(e^{t^3/3}-1\right)$
  being computed without difficulty.
- (e) Dividing by $2$, $\dot y + 6y = -e^t$. The homogeneous equation has
  general solution $\Gamma e^{-6t}$. The right-hand side not being constant,
  we look for a particular solution of the same form, \(y_2(t) = Ae^t\):
  substituting, $A + 6A = -1$, that is $A = -1/7$. The general solution is
  $y(t) = \Gamma e^{-6t} - \frac17e^t$ and the initial condition gives
  $\Gamma - \frac17 = \frac67$, that is \(\Gamma=1\):
  $y(t) = e^{-6t} - \frac17 e^t$. The solution diverges towards $-\infty$,
  even though the coefficient $a=6$ is positive: it is the right-hand side
  which diverges.
- (f) We look for an affine particular solution, \(y_2(t) = ct + d\):
  substituting, $c + ct + d = t$ imposes $c=1$ and $d=-1$. The general
  solution is $y(t) = t - 1 + \Gamma e^{-t}$, with $\Gamma = y(0)+1$. In the
  long run $y$ follows the line $t-1$, lagging $t$ by one unit.
#+END_details

* Nonlinear equations
:PROPERTIES:
:CUSTOM_ID: non-lineaires
:END:

In this section we are interested in differential equations which can be
written in the form:

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

or, equivalently, $\dot y(t) = h\bigl(y(t),t\bigr)$ with $h(y,t) =
-g(y,t)/f(y,t)$. One may recognise an exact differential equation, if the
condition $\partial f/\partial t = \partial g/\partial y$ holds; we already
know how to [[#exactes][solve]] this type of equation.

** Separable problems
:PROPERTIES:
:CUSTOM_ID: separables
:END:

We consider here the case where $f$ does not depend on $t$ and $g$ does not
depend on $y$. This class of problems is very easy to solve: each side is
integrated separately.

#+BEGIN_exemple
Consider the differential equation $3y^2\,\mathrm dy = t\,\mathrm dt$. Integrating
both sides, $\int 3y^2\,\mathrm dy = \int t\,\mathrm dt$, we directly obtain
$y^3 + \gamma_1 = \frac12t^2 + \gamma_2$, or $y^3 = \frac12t^2 + \gamma$, since
the two constants of integration are not separately identifiable. We
therefore have:
\[
y(t) = \left(\frac12t^2 + \gamma\right)^{\frac13}
\]
We can then fix $\gamma$ from the initial condition, evaluating the last
equation at $t=0$. We get $\gamma = y(0)^3$ and:
\[
y(t) = \left(\frac12t^2 + y(0)^3\right)^{\frac13}
\]
#+END_exemple

#+BEGIN_exemple
Consider the differential equation $2t\,\mathrm dy + y\,\mathrm dt = 0$. A priori,
it does not belong to the class discussed here, since $\mathrm dy$ is
associated with a function of $t$ and not of $y$. But dividing both sides by
$2yt$, we obtain:
\[
\frac{\mathrm dy}{y} + \frac{\mathrm dt}{2t} = 0
\]
Note in passing that this transformation makes the differential equation
exact, besides making it separable. Integrating, we get
$\log y + \frac12\log t = \gamma$, that is $y\,t^{\frac12} = e^{\gamma}$ and
finally:
\[
y(t) = \Gamma\,t^{-\frac12}
\]
This solution is not defined at zero, which we knew from the transformation
of the initial problem, which implicitly assumes $y\neq0$ and \(t\neq0\): to
fix the constant $\Gamma$, one must choose a strictly positive initial
instant. We recover, by a third route, the solution of the [[#exactes][example]]
treated with an integrating factor.
#+END_exemple

** Reduction to a linear dynamics: Bernoulli's equation
:PROPERTIES:
:CUSTOM_ID: bernoulli
:END:

#+BEGIN_definition
The differential equation of the form:
\[
\dot y(t) + R(t)\,y(t) = Q(t)\,y(t)^m
\]
with $m\notin\{0,1\}$, is a /Bernoulli equation/.
#+END_definition

For $m=0$ we recover the non-homogeneous linear equation, and for $m=1$ the
homogeneous linear equation with coefficient $R-Q$. In the other cases the
equation is nonlinear, but it can always be reduced to a linear differential
equation by a change of variable. Divide both sides by \(y^m\):

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

and set $z = y^{1-m}$, so that $\dot z = (1-m)y^{-m}\dot y$. The equation then
reads $\frac{1}{1-m}\dot z(t) + R(t)\,z(t) = Q(t)$, or:

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

This is a linear differential equation with variable coefficients. We can
therefore solve it, then express $y$ in terms of $z$ to obtain the solution
of the original equation.

#+NAME: ex-bernoulli
#+BEGIN_exemple
Consider the differential equation:
\[
\dot y(t) + t\,y(t) = 3t\,y(t)^2
\]
Dividing by $y^2$, we get $\dot y\,y^{-2} + t\,y^{-1} - 3t = 0$. Set
$z = y^{-1}$, so that \(\dot z = -y^{-2}\dot y\); the dynamics of $z$ reads
$-\dot z + tz - 3t = 0$, or equivalently:
\[
\dot z - t\,z + 3t = 0
\]
This is a linear equation with variable coefficients, with $a(t)=-t$ and
$b(t)=-3t$. Applying the [[#non-homogene-variable][general formula]], we directly get:
\[
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
\]
It remains to reverse the transformation, since the variable of interest is
$y$ and not $z$. We have $y = z^{-1}$ and, using the initial condition to
determine the constant, \(A = 1/y(0)\):
\[
y(t) = \frac{1}{e^{\frac{t^2}{2}}\left(\frac{1}{y(0)}-3\right) + 3}
\]
#+END_exemple

#+BEGIN_exercice
Find the solution of the following differential equation:
\[
\dot y(t) + \frac{1}{t}\,y(t) = y(t)^3
\]
with $t_0>0$ the initial instant, the initial condition $y(t_0)$ being known.
#+END_exercice

#+ATTR_HTML: :class corrige
#+BEGIN_details
#+BEGIN_summary
Solution
#+END_summary
This is a Bernoulli equation with $m=3$, $R(t)=1/t$ and $Q(t)=1$. We divide by
$y^3$ and set $z = y^{-2}$, so that $\dot z = -2y^{-3}\dot y$. We get
$-\frac12\dot z + \frac{z}{t} = 1$, that is:
\[
\dot z - \frac{2}{t}\,z = -2
\]
The integrating factor is $e^{-\int 2/t} = t^{-2}$, and the equation reads
$\frac{\mathrm d}{\mathrm dt}\left(t^{-2}z\right) = -2t^{-2}$, hence $t^{-2}z = \frac{2}{t} + C$
and $z(t) = 2t + Ct^2$. Going back to \(y\):
\[
y(t) = \pm\left(2t + Ct^2\right)^{-\frac12}
\]
the sign being that of $y(t_0)$, and the constant being fixed by the initial
condition, $C = \left(y(t_0)^{-2} - 2t_0\right)/t_0^2$. The solution is only
defined as long as $2t + Ct^2$ remains strictly positive: if $C < 0$, that
is if $y(t_0)^2 > 1/(2t_0)$, it blows up in finite time $t^{\star} = -2/C$.
One checks that $y = 0$ is a steady state, and that for $C\geq0$ the solution
converges to zero like $t^{-1}$ or $t^{-1/2}$.
#+END_details

** Qualitative approach
:PROPERTIES:
:CUSTOM_ID: qualitative
:END:

A nonlinear differential equation cannot always be solved analytically. In
many cases we have to settle for numerical solutions, for a quantitative
approach (see the [[#numerique][next]] section), or graphical ones, if a qualitative
approach is enough. With a graphical approach, one can for instance
investigate the stability of the dynamics.

We are interested here in a dynamics of the form:

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

where $f$ is a continuous function not depending on time. One then speaks of
an /autonomous/ differential equation. The dynamics can be represented
graphically in the $(y,\dot y)$ plane. The figure [[fig:phase][below]] gives two
examples: in the left chart the function $f$ is monotonically decreasing, in
the right one it is monotonically increasing. Reading these charts rests on
three observations:
1. when $f(y)>0$, the change $\dot y$ is positive and so $y$ increases;
2. when $f(y) < 0$, the change $\dot y$ is negative and so $y$ decreases;
3. when $f(y)=0$, the change $\dot y$ is zero and so $y$ does not move.
These three observations explain how we oriented the arrows on the
horizontal axis. These arrows describe the direction of change of the
variable $y$. For instance, in the left chart, $f(y)$ is positive when
$y < y_a^{\star}$ and negative when $y>y_a^{\star}$. Thus $y$ increases when
its level is low relative to $y_a^{\star}$, and the arrows point to the
right; $y$ decreases when its level is high, and the arrows point to the
left.

#+begin_src python :session edo1-en :exports none :results none
  def arrows(ax, segments, y=0.0):
      """Arrows on the horizontal axis, from x0 to x1 for each 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
  grid = np.linspace(0, 2, 200)
  fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
  for ax, sign, name in ((ax1, 1, 'a'), (ax2, -1, 'b')):
      ax.plot(grid, sign*fa(grid), 'b')
      ax.axhline(y=0, color='k', linewidth=0.8)
      ax.plot([1], [0], 'ko', markersize=5)
      ax.annotate(r'$y_%s^\star$' % name, xy=(1, 0), xytext=(1.03, 0.25))
      if sign > 0:
          arrows(ax, [(0.2, 0.75), (1.8, 1.25)])
      else:
          arrows(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('stable steady state')
  ax2.set_title('unstable steady state')
  fig.tight_layout()
  fig.savefig("edo-phase.svg", transparent=True)
#+end_src

#+CAPTION: *Phase diagram. On the left $f$ is decreasing and the steady state is stable, on the right $f$ is increasing and the steady state is unstable.*
#+LABEL: fig:phase
[[file:edo-phase.svg]]

$y_a^{\star}$ and $y_b^{\star}$ are the steady states, the fixed points of the
dynamics. A single glance at these charts tells us about their properties:
$y_a^{\star}$ is stable (if for some reason one moves away from $y_a^{\star}$,
one comes back to it), unlike $y_b^{\star}$ which is unstable (if one moves
away from $y_b^{\star}$, one never comes back).

These charts suggest that if $f$ is monotonically decreasing then the
dynamics is stable, and that if $f$ is monotonically increasing then the
dynamics is unstable. In fact we should distinguish global stability from
local stability. The figure [[fig:phase-multiple][below]] shows that the function $f$ can be
much more "twisted". In that case $y_a^{\star}$ and $y_c^{\star}$ are locally
stable steady states and $y_b^{\star}$ is an unstable steady state.
$y_a^{\star}$ is /locally/ stable, in the sense that if one moves moderately
away from $y_a^{\star}$ one comes back to it, but if one moves too far,
beyond $y_b^{\star}$, one never returns to $y_a^{\star}$.

#+begin_src python :session edo1-en :exports none :results none
  fm = lambda y: -0.5*(y-1)*(y-2)*(y-3)
  grid = np.linspace(0.45, 3.55, 300)
  fig, ax = plt.subplots(figsize=(7, 4))
  ax.plot(grid, fm(grid), 'b')
  ax.axhline(y=0, color='k', linewidth=0.8)
  for yst, name in ((1, 'a'), (2, 'b'), (3, 'c')):
      ax.plot([yst], [0], 'ko', markersize=5)
      ax.annotate(r'$y_%s^\star$' % name, xy=(yst, 0), xytext=(yst+0.04, 0.12))
  arrows(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: *Phase diagram with several steady states. $y_a^{\star}$ and $y_c^{\star}$ are locally stable, $y_b^{\star}$ is unstable.*
#+LABEL: fig:phase-multiple
[[file:edo-phase-multiple.svg]]

Two ideas are worth retaining from these examples:
1. the steady state, if it exists, is not always unique. One can also
   imagine cases where it does not exist;
2. stability, local or global, can only arise if the function $f$, called
   the /transition function/, is, at least locally, decreasing.

In retrospect, point 2 is consistent with what we saw in the linear case.
Consider the linear differential equation $\dot y = b - ay$. The stability
properties depend on the sign of the parameter \(a\): the dynamics is stable
if and only if $a$ is positive, that is if there is a decreasing relation
between the change and the level. In that case, for any initial condition,
$y(t)$ converges to $b/a$.

** Numerical approach
:PROPERTIES:
:CUSTOM_ID: numerique
:END:

The approach described in the previous section may not suffice if one wants
more precise information about the solution. One then resorts to numerical
computation. This task is relatively simple with scientific software; we use
here the =solve_ivp= function of the =scipy= library, which implements by
default an adaptive-step Runge-Kutta method of order 4(5), the same as
Matlab's =ode45= function. These functions solve problems of the form:

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

where the initial condition $y(t_0)$ is given, over the time interval
$[t_0,t_1]$.

As an example we shall numerically solve the nonlinear differential equation
of the [[ex-bernoulli][example]] treated above; since we have an analytical solution, we
will be able to assess the accuracy of the numerical solution. As a
reminder, the analytical solution is:

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

Before embarking on the numerical study of the dynamics, the properties of
this solution deserve discussion. We observe the following points.
1. $y^{\star} = 1/3$ is a steady state. To check it, one can substitute
   $y(t) = 1/3$ into the differential equation. If $y(0) = 1/3$ then
   $y(t) = 1/3$ for all $t$.
2. If $y(0) < 1/3$ then $\lim_{t\to\infty}y(t) = 0^+$. The steady state
   $y^{\star}$ is therefore not stable; it is zero, the other steady state,
   which attracts these paths.
3. The case $y(0) > 1/3$ deserves special treatment. The denominator
   $3 - e^{t^2/2}\left(3 - 1/y(0)\right)$ is then a decreasing function of
   $t$, positive at $t=0$, which vanishes at:
   \[
   t^{\star} = \sqrt{2\log\frac{3}{3-\frac{1}{y(0)}}}
   \]
   The solution has a vertical asymptote at \(t^{\star}\): for all
   $t < t^{\star}$, $y(t)$ is monotonically increasing with $\lim_{t\to t^{\star}}y(t) = +\infty$.
   The path ceases to exist at $t^{\star}$. The formula remains defined
   beyond, but it then describes another solution of the equation,
   negative, increasing from $-\infty$ to $0^-$, which is not connected to
   the initial condition.
The properties of the solution depend radically on the initial condition.
The numerical approach of the case $y(0) < 1/3$ poses no problem, unlike the
complementary case, because of the vertical asymptote.

We start by defining the function $f(t,y)$ associated with the differential
equation, written in the form $\dot y(t) = 3t\,y(t)^2 - t\,y(t)$, then the
exact solution:

#+begin_src python :session edo1-en :exports code :results none
  def bernoulli(t, y):
      """Right-hand side of the differential equation dy/dt = 3ty² - ty."""
      return 3*t*y**2 - t*y

  def exact(t, y0):
      """Analytical solution for the initial condition y0."""
      return 1/(np.exp(t**2/2)*(1/y0 - 3) + 3)
#+end_src

The following script solves the equation for two initial conditions, on
either side of the steady state. In the case $y_0>y^{\star}$, the terminal
instant is set short of $t^{\star}$, since no solver can follow the solution
up to a vertical asymptote. The =solve_ivp= function returns an object whose
=t= attribute contains the instants at which the solver evaluated the
solution, and whose =sol= attribute, obtained thanks to the =dense_output=
option, allows the solution to be interpolated at any point of the
interval.

#+begin_src python :session edo1-en :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-en :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='numerical solution')
      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 - exact(t, y0)), 'b', label='default tolerances')
      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] - exact(t, y0)), 'g', label=r'$\mathtt{rtol}=10^{-8}$')
      ax.set_xlabel(r'$t$')
      ax.set_ylabel('absolute error')
      ax.legend(loc='center right')
  fig.tight_layout()
  fig.savefig("edo-bernoulli.svg", transparent=True)
#+end_src

#+CAPTION: *Numerical solutions. Top, the interpolated solution (line) and the points computed by the solver (dots) for $y_0 = 0.9\,y^{\star}$ (left) and $y_0 = 1.1\,y^{\star}$ (right, the vertical asymptote is at $t^{\star} = 2.1899$). Bottom, the absolute error with respect to the analytical solution, with the default tolerances and with a relative tolerance of $10^{-8}$.*
#+LABEL: fig:bernoulli
[[file:edo-bernoulli.svg]]

The figure [[fig:bernoulli][above]] shows the results obtained for $y_0 = 0.9\,y^{\star}$
and $y_0 = 1.1\,y^{\star}$. With the default tolerances, a relative tolerance
of $10^{-3}$ and an absolute tolerance of $10^{-6}$, the solver makes do with
$29$ steps in the first case and $8$ steps in the second: the error is at
most $4.3\times10^{-5}$ in absolute value in the first case, and
$5.7\times10^{-4}$ in the second. One observes that the errors are almost
zero where the path of $y$ is flat and become larger where the path is
steeper, in particular as the asymptote is approached. These errors shrink
if the solver is asked for a more accurate approximation: with =rtol=1e-8=
and =atol=1e-11=, the number of steps rises to $127$ and $32$, and the
maximum error drops to $1.2\times10^{-8}$ and $3.8\times10^{-8}$. Accuracy
is paid for in evaluations of the =bernoulli= function, which is of no
consequence here but may be for large systems.

This application teaches us that it is always useful to investigate the
properties of the solution before embarking on a numerical computation, with
the analytical solution when it is available (but then the numerical
solution becomes pointless, except as a check) or with the qualitative
approach. Some preliminary thought is always a good idea. In the case
$y_0>y^{\star}$, the solver fails if one tries to obtain the solution for
$t\in[0,T]$ with $T$ too close to $t^{\star}$, where the path is very steep,
or a fortiori beyond \(t^{\star}\); and there is no guarantee that it will
signal its failure other than by aberrant values.

#+BEGIN_remarque
Adaptive-step solvers such as =solve_ivp= choose the length of the time
steps themselves so as to meet the requested tolerances. Fixed-step schemes,
in particular the Euler scheme, are presented in the note on [[https://stephane-adjemian.fr/en/posts/the-solow-model-in-discrete-time/][the Solow model
in discrete time]]; the note on the [[https://stephane-adjemian.fr/en/posts/simulating-the-solow-model/][simulation of the Solow model]] uses the
=odeint= function of =scipy=, which relies on another family of methods.
#+END_remarque

** Local approximation
:PROPERTIES:
:CUSTOM_ID: locale
:END:

As in the section on the [[#qualitative][qualitative approach]], we are interested here in an
autonomous dynamics $\dot y = f(y)$, where $f$ is a continuous and
differentiable function not depending on time. So far we have only
considered global approaches. The qualitative approach and the analytical
approach are global, because we characterise or obtain the solution for all
possible values of $t$ and independently of the level of $y$. The numerical
approach is also global, because it is valid on an interval of values of
$t$ independently of the level of $y$. Here we take a local approach, in the
sense that what we will be able to say about $y$ will only be relevant in a
neighbourhood of a specific level of the variable.

If the function $f$ is differentiable at $\bar y$, a first-order Taylor
expansion gives:

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

If we drop the residual term, which tends to zero faster than $|y-\bar y|$
as $y$ approaches $\bar y$, we have:

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

In a neighbourhood of $\bar y$, the dynamics is approximately linear. The
quality of the approximation depends on the distance from $y$ to $\bar y$.
We know how to solve this differential equation, it is an equation with
constant coefficients. But one must bear in mind that the solution we will
obtain by applying the results of the [[#constants][first section]] will only be
acceptable for values of $y$ in a neighbourhood of $\bar y$. It is in this
sense that this approach is local.

A priori we can choose any point $\bar y$, provided the function $f$ is
differentiable there. Usually the model is approximated around a steady
state of the dynamics, $y^{\star}$. This eliminates the constant, since by
definition $\dot y$ is zero at the steady state. We thus have:

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

Another motivation is that if this steady state is stable, then we are sure
that if $y(0)$ is in a neighbourhood of $y^{\star}$, $y(t)$ will remain in
that same neighbourhood, which guarantees the quality of the approximation.
Note that in the case of multiple steady states, there is a degree of
freedom in the choice of the point around which $f$ is approximated. The
properties of the approximate solution can be very different, for instance
in terms of stability.

The geometric interpretation of this approximation is direct. Making a
first-order approximation, also called a /linearisation/, means replacing
the function $f$ by its tangent at $y^{\star}$. The figure [[fig:linearisation][below]]
illustrates this replacement, for the transition function of the Solow model
which we shall meet again in the [[#solow][next section]].

#+begin_src python :session edo1-en :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))
  slope = mu*(al-1)
  grid = np.linspace(0, 5, 300)
  fig, ax = plt.subplots(figsize=(7, 4))
  ax.plot(grid, fs(grid), 'b', label=r'$f(y)$')
  ax.plot(grid, slope*(grid-kst), 'r', linewidth=1, linestyle='--', label=r"tangent at $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], [-slope*kst], 'ro', markersize=4)
  ax.annotate(r"$-f'(y^\star)\,y^\star$", xy=(0, -slope*kst), xytext=(0.15, -slope*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: *Graphical interpretation of the linearisation. The transition function $f(y) = sy^{\alpha} - (n+g+\delta)y$ of the Solow model and its tangent at the steady state.*
#+LABEL: fig:linearisation
[[file:edo-linearisation.svg]]

Finally, note that a priori nothing prevents us from considering orders of
approximation higher than one. In some cases it may be necessary, so as not
to miss interesting properties of the variable under study, to go for a
second- or third-order approximation. The second-order expansion around the
steady state:

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

is a Bernoulli equation with constant coefficients, with $m=2$, whose
explicit solution we therefore know how to exhibit. This is what we shall do
for the Solow model.

** Application: the Solow model
:PROPERTIES:
:CUSTOM_ID: solow
:END:

The Solow model is a continuous-time dynamic model describing the evolution
of the stock of physical capital, or of output, in a closed economy. We
shall not go into it in detail here, it is presented in the note on [[https://stephane-adjemian.fr/en/posts/the-solow-model/][the Solow
model]]; the aim is simply to apply the methods presented above.

We take a neoclassical production function[fn:: A function homogeneous of
degree one, increasing and concave in each of its arguments, and satisfying
the Inada conditions. See the note on the Solow model.],
$Y(t) = F\bigl(K(t), A(t)L(t)\bigr)$, with $L(t) = e^{nt}$ the population and
$A(t) = e^{gt}$ the index of labour efficiency[fn:: The other notes on this
site denote by $x$ the growth rate of technical progress.]. The law of
motion of the stock of physical capital is:

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

with $s\in\,]0,1[$ the exogenous saving rate and $\delta\in[0,1]$ the
depreciation rate of capital. This equation simply tells us that the capital
stock increases if and only if investment dominates depreciation.
Substituting the production function, one sees clearly that this is a
first-order nonlinear differential equation:

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

This differential equation is non-autonomous: the relation between the
change in the capital stock and its level depends on time, through the term
$e^{(n+g)t}$ which represents the growth of efficient labour. The first thing
to do is to reduce it to an autonomous differential equation. To this end,
it suffices to remove the demographic and technological trends. Let
$\hat k(t) = K(t)/\bigl(A(t)L(t)\bigr)$ be the capital stock per efficient
worker, "purged" of population growth and technical progress. One easily
checks that the dynamics of capital per efficient worker is given by:

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

where $f(\hat k) = F(\hat k, 1) = \hat y$ is output per efficient worker. This
is indeed an autonomous differential equation. If the production function is
neoclassical, there exists a unique strictly positive steady state
$\hat k^{\star}$. In passing, note that $\hat k=0$ is also a steady state,
called trivial, which we shall set aside. To be convinced of this existence
and uniqueness result, it suffices to rewrite the equation in the form:

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

The growth rate of capital per efficient worker is strictly positive if and
only if gross investment per unit of capital strictly exceeds the effective
depreciation rate $n+g+\delta$. The first term on the right-hand side equals
the average product of capital multiplied by the saving rate. It is strictly
positive, monotonically decreasing, because the returns to capital are
decreasing, and tends to infinity at zero and to zero at infinity by the
Inada conditions. Hence the curve of gross investment per unit of capital
necessarily crosses the horizontal line at height \(n+g+\delta\) exactly
once: the steady state is unique, and it satisfies:

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

In what follows we shall seek to describe the dynamics of output per
efficient worker, and more specifically the dynamics of adjustment towards
the steady state, by computing the /speed of convergence/. Differentiating
$\hat y(t) = f(\hat k(t))$ with respect to $t$, we obtain:

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

with $r(\hat k) = f'(\hat k)$ the marginal product of capital, which in a
perfectly competitive environment corresponds to the real interest rate, and
$\alpha(\hat k) = \hat kf'(\hat k)/f(\hat k)$ the elasticity of output with
respect to capital, or, in a perfectly competitive environment, the share of
capital income in total income. By definition of $\alpha$, the steady state
also satisfies $\hat y^{\star}/\hat k^{\star} = r^{\star}/\alpha^{\star}$.
Finally, we shall use the elasticity of substitution between capital and
labour, defined [[https://stephane-adjemian.fr/en/posts/properties-of-the-ces-production-function/][here]], which in terms of the intensive technology reads:

\[
\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)}
\]

We shall sometimes assume that the technology is of the CES type (for
/constant elasticity of substitution/). This production function is more
general than the Cobb-Douglas function, but it is not neoclassical, the
Inada conditions not being satisfied. In that case the existence of the
steady state $\hat k^{\star}>0$ is no longer guaranteed and depends on the
values of the parameters. We set:

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

with $\rho\in\,]-\infty,1]$, $\sigma = (1-\rho)^{-1}$, that is
$\rho = 1-1/\sigma$, the constant elasticity of substitution between the
factors, and $\gamma_1+\gamma_2=1$. We recover the Cobb-Douglas function,
which is neoclassical, when $\rho\to0$, that is $\sigma\to1$. The average
product $f(\hat k)/\hat k = \left(\gamma_1 + \gamma_2\hat k^{-\rho}\right)^{1/\rho}$
is still decreasing, but it no longer spans the whole of $]0,\infty[$, hence
the existence conditions of the steady state:
1. if $\rho\in\,]0,1]$, the factors being more substitutable than in the
   Cobb-Douglas case, the average product decreases from infinity to
   $\gamma_1^{1/\rho}$, and the unique steady state $\hat k^{\star}>0$ exists
   if and only if $s\gamma_1^{1/\rho} < n+g+\delta$, that is if the saving
   rate is not too high;
2. if $\rho < 0$, the factors being less substitutable than in the
   Cobb-Douglas case, the average product decreases from $\gamma_1^{1/\rho}$
   to zero, and the unique steady state $\hat k^{\star}>0$ exists if and only
   if $s\gamma_1^{1/\rho} > n+g+\delta$, that is if the saving rate is not too
   low;
3. if $\rho=0$ the production function is Cobb-Douglas, and the existence and
   uniqueness of the steady state are guaranteed whatever the values of the
   parameters.
With this technology, the interest rate and the share of capital income are
written as functions of \(\hat y\):

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

The capital share is an increasing function of output per efficient worker
if and only if $\rho>0$, that is if and only if the elasticity of
substitution exceeds one. At the steady state,
$\alpha^{\star} = \gamma_1\left(s/(n+g+\delta)\right)^{\rho}$.

Throughout what follows we use the calibration of the previous notes:
$s=0.20$, $n=0.02$, $g=0.02$, $\delta=0.10$, and a capital share at the
steady state $\alpha^{\star}=0.36$. To compare technologies for a given
steady state, we set $\gamma_1$ so that $\alpha^{\star}=0.36$ whatever the
elasticity of substitution:
$\gamma_1 = \alpha^{\star}\left((n+g+\delta)/s\right)^{\rho}$ and
$\gamma_2 = 1-\gamma_1$.

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

  def technology(sigma):
      """Intensive CES production function, its derivative, the elasticity of
      output with respect to capital and the steady state, for an elasticity of
      substitution sigma, with alpha* = 0.36 whatever 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):
      """Exact path of capital per efficient worker, by numerical integration."""
      f, fprime, alpha, kstar = technology(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

*** Characterising the adjustment without approximation
:PROPERTIES:
:CUSTOM_ID: sans-approximation
:END:

The speed of convergence is defined as the opposite of the growth rate of
the growth rate of output per efficient worker:

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

It measures the speed at which the growth rate dies out, as the economy
approaches its steady state. It will be convenient, here and in the
following sections, to reason on the logarithmic deviation of capital from
its steady-state level, $u(t) = \log\bigl(\hat k(t)/\hat k^{\star}\bigr)$, whose
dynamics reads:

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

with $G(0)=0$. Two derivatives will keep coming back. First, noting that
\(\mathrm d\hat k/\mathrm du = \hat k\):

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

Second, differentiating $\log\alpha = \log\hat k + \log f' - \log f$ with
respect to $u$ and using the definition of the elasticity of substitution:

\[
\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
The speed of convergence of output per efficient worker is, at any instant
of the transition:
\[
\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
The growth rate of output is $g_{\hat y} = \alpha(\hat k)\,g_{\hat k}$ and
$g_{\hat k} = \dot u = G(u)$. Differentiating with respect to time,
$\dot g_{\hat y} = \left[\frac{\mathrm d\alpha}{\mathrm du}G(u) + \alpha G'(u)\right]G(u)$,
hence:
\[
\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)
\]
that is, collecting the terms in \(sf/\hat k\):
\[
\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]
\]
It remains to note that, by the steady-state condition,
$\frac{sf(\hat k)/\hat k}{n+g+\delta} = \frac{\hat y/\hat k}{\hat y^{\star}/\hat k^{\star}}$.
#+END_proof

For a Cobb-Douglas technology, the elasticity of substitution is unitary and
$\hat y/\hat y^{\star} = (\hat k/\hat k^{\star})^{\alpha}$, so that we recover
the result of Barro and Sala-i-Martin (1995, appendix to chapter 1):

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

The speed of convergence decreases along the transition if and only if the
economy reaches its steady state from below, that is if $g_{\hat y}>0$. At
the steady state we recover the standard result, established in the [[#ordre-1][next
section]] by linearising the model:
$\beta_1(t)\to\beta^{\star} = (1-\alpha)(n+g+\delta)$. The note on the
[[https://stephane-adjemian.fr/en/posts/simulating-the-solow-model/][simulation of the Solow model]] plots this speed along the transition.

More generally, when the production function is of the CES type, we have
$\hat y/\hat k = \gamma_1^{1/\rho}\alpha(\hat y)^{-1/\rho}$, since
$\alpha = \gamma_1(\hat k/\hat y)^{\rho}$, and the speed of convergence can be
expressed in terms of the capital share alone:

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

with $\alpha(\hat y) = 1-\gamma_2\hat y^{-\rho}$. The direction of change of
the capital share during the transition depends on the sign of \(\rho\): it
increases with output if $\sigma>1$ and decreases if $\sigma < 1$. It would
be tempting to deduce the direction of change of the speed of convergence
from the factor \(1-\alpha(\hat y)\) alone: an increasing capital share would
slow convergence down, a decreasing share would speed it up, so that the
direction of change of $\beta_{\sigma}$ would depend on the sign of $\rho$.
But the bracketed term also varies, in the opposite direction, and it is the
one that prevails.

#+BEGIN_property
In a neighbourhood of the steady state, the speed of convergence varies with
the deviation $u = \log(\hat k/\hat k^{\star})$ according to:
\[
\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]
\]
This derivative is negative as soon as $\alpha^{\star} < 1/2$, whatever the
elasticity of substitution. An economy below its steady state therefore
converges faster than $\beta^{\star}$, and its speed of convergence decreases
during the transition; an economy above converges more slowly than
$\beta^{\star}$, and its speed increases.
#+END_property

#+BEGIN_proof
Write $\beta = (1-\alpha)\left[\frac{\phi}{\sigma} + (n+g+\delta)\left(1-\frac1\sigma\right)\right]$
with $\phi = sf(\hat k)/\hat k$. We have $\mathrm d\phi/\mathrm du = G'(u) = -\phi(1-\alpha)$ and
$\mathrm d(1-\alpha)/\mathrm du = -\alpha(1-\alpha)(1-1/\sigma)$. At the steady state
$\phi = n+g+\delta$, and the bracket equals $n+g+\delta$, hence:
\[
\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]
\]
The first term is the effect of the capital share, positive when
\(\sigma < 1\); the second is the effect of the bracket, always negative, and
it dominates the first if and only if $1-2\alpha+\alpha\sigma>0$, which is
guaranteed for $\alpha < 1/2$.
#+END_proof

The figure [[fig:vitesse][below]] confirms this result, and shows what happens far from the
steady state. The three technologies share the same steady state and the
same asymptotic speed $\beta^{\star} = 0.0896$, by construction of the
calibration. From below, the speed of convergence starts from a level that
is higher the more substitutable the factors, and decreases towards
\(\beta^{\star}\); for $\sigma=0.5$ it is almost flat during the first years,
the effect of the capital share then almost exactly offsetting that of the
bracket. From above, the speed increases towards $\beta^{\star}$ in all three
cases. For $\sigma=0.5$ it is even negative during the first years: output
decreases faster and faster before slowing down. Nothing paradoxical, since
$g_{\hat y} = \alpha\,g_{\hat k}$ and, far above the steady state, the capital
share is very low ($0.12$ against $0.36$ at the steady state) and rises
quickly, while the rate of decline of capital fades only slowly.

#+begin_src python :session edo1-en :exports none :results none
  def speed(sigma, k):
      """Speed of convergence of output per efficient worker along a path."""
      f, fprime, alpha, kstar = technology(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, colour in ((0.5, 'r'), (1.0, 'b'), (2.0, 'g')):
          kstar = technology(sigma)[3]
          k = transition(sigma, ratio*kstar).sol(tg)[0]
          ax.plot(tg, speed(sigma, k), colour, label=r'$\sigma = %s$' % str(sigma))
      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$ (years)')
      ax.set_ylabel(r'$\beta(t)$')
      ax.legend()
  ax1.set_title(r'start at $0.25\,\hat k^\star$')
  ax2.set_title(r'start at $4\,\hat k^\star$')
  fig.tight_layout()
  fig.savefig("solow-vitesse.svg", transparent=True)
#+end_src

#+CAPTION: *Speed of convergence of output per efficient worker along the transition, for three elasticities of substitution, starting below (left) and above (right) the steady state.*
#+LABEL: fig:vitesse
[[file:solow-vitesse.svg]]

We retain that the speed of adjustment towards the steady state is not
constant, but changes along the transition, and that its direction of change
depends first of all on the initial position relative to the steady state;
the elasticity of substitution between the factors modulates the size of
these variations, and can change their shape far from the steady state. In
this section we have not solved any differential equation, for instance to
compute the level of output per efficient worker at instant $t$. This is
generally not possible[fn:: Except for some production functions. For a
Cobb-Douglas or Leontief technology the computations can be carried through;
we shall come back to it.], and one must then resort to a numerical
approach, as we did to draw the figure, or consider an approximation of the
differential equation. Note nevertheless that, without even trying to solve
the differential equation, we were able to learn fairly precise things about
the transition dynamics.

*** Characterising the adjustment with a linearisation of the dynamics
:PROPERTIES:
:CUSTOM_ID: ordre-1
:END:

Let us now apply the [[#locale][local approximation]] to the dynamics of the deviation
$u = \log(\hat k/\hat k^{\star})$, whose steady state is $u=0$.

#+BEGIN_property
In a neighbourhood of the steady state, the dynamics of capital per
efficient worker is approximately linear:
\[
\dot u \approx -\beta\,u
\qquad\text{with}\qquad
\beta = (n+g+\delta)(1-\alpha^{\star})
\]
and the same holds for output per efficient worker:
\[
\frac{\mathrm d\log\hat y(t)}{\mathrm dt} \approx -\beta\log\frac{\hat y(t)}{\hat y^{\star}}
\]
For a CES production function,
$\beta = (n+g+\delta)\left[1-\gamma_1\left(\frac{s}{n+g+\delta}\right)^{\rho}\right]$.
#+END_property

#+BEGIN_proof
The first-order Taylor expansion of $G$ around zero gives
$\dot u \approx G'(0)\,u$, and $G'(0) = -s\frac{f(\hat k^{\star})}{\hat k^{\star}}(1-\alpha^{\star}) = -(n+g+\delta)(1-\alpha^{\star})$
by the steady-state condition. For output, the first-order expansion of
$\log f(\hat k^{\star}e^{u})$ gives $\log(\hat y/\hat y^{\star}) \approx \alpha^{\star}u$,
since the derivative of $\log f(\hat k^{\star}e^{u})$ with respect to $u$ is
\(\alpha(\hat k)\); thus
$\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})$.
Finally, for the CES, $\alpha^{\star} = \gamma_1(\hat k^{\star}/\hat y^{\star})^{\rho} = \gamma_1\left(s/(n+g+\delta)\right)^{\rho}$.
#+END_proof

By considering a first-order Taylor expansion around the steady state, we
lose information about the transition dynamics: here the speed of adjustment
towards the steady state is constant. Note nevertheless that the speed of
convergence obtained here is the limit, as $t$ tends to infinity, of the
speed of convergence obtained in the previous section. This asymptotic
equivalence is not surprising, since the stability of the steady state
guarantees that, once $t$ is large enough, output per efficient worker is
arbitrarily close to the steady state. Finally, note that we were able to
carry out the computations without specifying the production function: this
is generally the case when an approximation of the model is considered.

We have lost part of the properties of the transition, but we can now solve
the dynamics, that is, compute the level of output per efficient worker at
any instant. Obviously, this solution will only be valid in a neighbourhood
of the steady state. The linear equation with constant coefficient tells us
that the growth rate of the distance to the steady state is negative and
constant: for any initial condition, the distance is absorbed in infinite
time, $u(t) = u(0)e^{-\beta t}$, and likewise for output. Substituting the
definition of the deviation, we obtain:

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

Output per efficient worker at instant $t$ is a weighted geometric mean of
its initial value and of its steady-state value, the weight of the latter
tending to one. For our calibration, $\beta^{\star} = 0.0896$ and the
half-life of the deviation from the steady state is
$\log 2/\beta^{\star} = 7.7$ years.

*** Characterising the adjustment with a second-order approximation
:PROPERTIES:
:CUSTOM_ID: ordre-2
:END:

One could expand the dynamics of output to second order, as we just did to
first order, but the second derivative is then very cumbersome and does not
lend itself to a solution. By reasoning on the capital deviation $u$ rather
than on the output deviation, the computation takes a few lines, and the
approximate equation can be solved explicitly.

#+BEGIN_property
In a neighbourhood of the steady state, the dynamics of the deviation
$u = \log(\hat k/\hat k^{\star})$ is, to second order:
\[
\dot u \approx -\beta\,u + \gamma\,u^2
\qquad\text{with}\qquad
\gamma = \frac{\beta}{2}\left(1-\frac{\alpha^{\star}}{\sigma^{\star}}\right)
\]
where $\sigma^{\star} = \sigma(\hat k^{\star})$ is the elasticity of
substitution at the steady state.
#+END_property

#+BEGIN_proof
The second-order Taylor expansion of $G$ around zero reads
$\dot u \approx G'(0)u + \frac12G''(0)u^2$. We know $G'(0) = -\beta$.
Differentiating $G'(u) = -\phi(\hat k)\bigl[1-\alpha(\hat k)\bigr]$, with
$\phi = sf/\hat k$, we get:
\[
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]
\]
At the steady state, $\phi = n+g+\delta$ and
$G''(0) = (n+g+\delta)(1-\alpha^{\star})\left(1-\alpha^{\star}/\sigma^{\star}\right) = 2\gamma$.
#+END_proof

The quadratic term vanishes if and only if $\sigma^{\star}=\alpha^{\star}$.
For a Cobb-Douglas technology, \(\gamma = \frac{\beta}{2}(1-\alpha)>0\); more
generally, $\gamma$ is positive as soon as $\sigma^{\star}>\alpha^{\star}$,
which is the case for all usual calibrations. The approximate equation is a
[[#bernoulli][Bernoulli equation]] with constant coefficients, with $m=2$, $R=\beta$ and
$Q=\gamma$, and we know how to solve it.

#+BEGIN_property
The solution of the second-order approximation is:
\[
u(t) = \frac{u(0)\,e^{-\beta t}}{1 - \frac{\gamma}{\beta}\,u(0)\left(1-e^{-\beta t}\right)}
\]
and output per efficient worker follows, to the same order, from:
\[
\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
We divide the equation by $u^2$ and set $z = 1/u$, so that
$\dot z = -\dot u/u^2$. We get $\dot z = \beta z - \gamma$, a linear equation
with constant coefficients whose steady state is $\gamma/\beta$ and whose
solution is $z(t) = \left(z(0) - \frac{\gamma}{\beta}\right)e^{\beta t} + \frac{\gamma}{\beta}$.
Going back to $u = 1/z$ with $z(0) = 1/u(0)$, and multiplying numerator and
denominator by $u(0)e^{-\beta t}$, we obtain the stated expression. For
output, the second derivative of $\log f(\hat k^{\star}e^{u})$ with respect
to $u$ is $\mathrm d\alpha/\mathrm du = \alpha(1-\alpha)(1-1/\sigma)$, which gives
the second-order expansion.
#+END_proof

One can go one step further and write directly, to second order, the
dynamics of the output deviation $v = \log(\hat y/\hat y^{\star})$. It has the
same form as that of capital, and its speed of convergence compares directly
with that of the section [[#sans-approximation][without approximation]].

#+BEGIN_property
In a neighbourhood of the steady state, the output deviation
$v = \log(\hat y/\hat y^{\star})$ follows, to second order, the Bernoulli
equation:
\[
\dot v \approx -\beta\,v + \gamma_y\,v^2
\qquad\text{with}\qquad
\gamma_y = \frac{\beta}{2\alpha^{\star}}\left[\alpha^{\star} + \frac{1-2\alpha^{\star}}{\sigma^{\star}}\right]
\]
whose solution is $v(t) = v(0)e^{-\beta t}\big/\bigl[1-\frac{\gamma_y}{\beta}v(0)(1-e^{-\beta t})\bigr]$.
The speed of convergence of output implied by the approximation is:
\[
-\frac{\ddot v}{\dot v} = \beta - 2\gamma_y\,v
\]
#+END_property

#+BEGIN_proof
Let $c = \frac12\alpha^{\star}(1-\alpha^{\star})(1-1/\sigma^{\star})$, so that
$v \approx \alpha^{\star}u + cu^2$ and $\dot v \approx (\alpha^{\star} + 2cu)\dot u$.
Substituting $\dot u \approx -\beta u + \gamma u^2$ and keeping only the terms
of order at most two, $\dot v \approx -\alpha^{\star}\beta\,u + (\alpha^{\star}\gamma - 2c\beta)u^2$.
It remains to express $u$ in terms of $v$ to the same order, by inverting the
expansion: $u \approx v/\alpha^{\star} - (c/\alpha^{\star 3})v^2$. We get
$\dot v \approx -\beta v + \bigl[\gamma/\alpha^{\star} - \beta c/\alpha^{\star 2}\bigr]v^2$
and, replacing $\gamma$ and $c$ by their expressions:
\[
\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]
\]
The equation obtained is that of capital with $\gamma_y$ in place of
$\gamma$, hence the solution. Finally $\ddot v = (-\beta + 2\gamma_y v)\dot v$.
#+END_proof

Two consistency checks. First, since $v \approx \alpha^{\star}u$, the slope of
this speed with respect to $u$ at the steady state is
\(-2\gamma_y\alpha^{\star} = -(n+g+\delta)(1-\alpha^{\star})\bigl[\alpha^{\star} + (1-2\alpha^{\star})/\sigma^{\star}\bigr]\):
this is exactly the derivative established in the section without
approximation. The second-order speed is the tangent to the exact speed
$\beta_{\sigma}$ at the steady state, where the first order only retains
the value $\beta$. Second, in the Cobb-Douglas case,
$\gamma_y = \beta(1-\alpha)/(2\alpha)$, and the exact speed
$\beta_1 = \beta\,e^{-\frac{1-\alpha}{\alpha}v}$ does have $\beta - 2\gamma_y v$
as its first-order expansion in $v$.

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

  def order2_output(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-en :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, colour in ((0.5, 'r'), (1.0, 'b'), (2.0, 'g')):
          f, fprime, alpha, kstar = technology(sigma)
          k = transition(sigma, ratio*kstar).sol(tg)[0]
          v0 = np.log(f(ratio*kstar)/f(kstar))
          ax.plot(tg, speed(sigma, k), colour, label=r'$\sigma = %s$' % str(sigma))
          ax.plot(tg, betastar - 2*gamma_y(sigma)*order2_output(tg, v0, sigma), colour,
                  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$ (years)')
      ax.set_ylabel(r'$\beta(t)$')
      ax.legend()
  ax1.set_title(r'start at $0.25\,\hat k^\star$')
  ax2.set_title(r'start at $4\,\hat k^\star$')
  fig.tight_layout()
  fig.savefig("solow-vitesse-ordre2.svg", transparent=True)
#+end_src

#+CAPTION: *Speed of convergence of output per efficient worker: exact dynamics (solid, as in the figure of the section without approximation) and second-order approximation, $\beta - 2\gamma_y v(t)$ evaluated along the approximate solution (dashed). The first order gives the constant $\beta^{\star}$ (dotted).*
#+LABEL: fig:vitesse-ordre2
[[file:solow-vitesse-ordre2.svg]]

The figure [[fig:vitesse-ordre2][above]] superimposes the two speeds along the transitions of the
section without approximation. Near the steady state, the second order
reproduces the direction and the order of magnitude of the deviations from
$\beta^{\star}$, in both directions and for all three technologies, which the
first order cannot do. Far from the steady state, the approximation, which is
only the tangent to the exact speed, departs from it. For $\sigma=1$ and
$\sigma=2$ it underestimates the speed in both directions, the exact speed
being convex in $v$ and the tangent lying below it. For $\sigma=0.5$ it
greatly exaggerates the speed from below, where the exact speed flattens out
because the transition function changes curvature along the path, as we
shall see below; and it cannot reproduce the negative speed observed from
above, which is an effect of the variation of the capital share far from its
steady-state level.

Three remarks on this solution.

#+BEGIN_remarque
For $\gamma=0$ we recover the first-order solution, and in general the
denominator corrects the first-order solution by a factor which depends on
the sign of $u(0)$. If $\gamma>0$ and the economy starts below its steady
state ($u(0) < 0$), the denominator exceeds one and $|u(t)|$ is smaller than
at first order: convergence is faster. If it starts above, the denominator
is below one and convergence is slower. This is exactly what the
[[fig:vitesse][previous]] figure says for the exact dynamics. It can also be read on the
speed of adjustment of capital. Since $\dot u$ is the growth rate of capital
per efficient worker, $g_{\hat k}$, the ratio
$-\ddot u/\dot u = -\dot g_{\hat k}/g_{\hat k}$ is the counterpart, for
capital, of the speed of convergence defined for output in the section
[[#sans-approximation][without approximation]]: it measures the speed at which the growth rate of
capital dies out, and equals exactly $\beta$ for the linear dynamics
$\dot u = -\beta u$. Since $\dot u = G(u)$, we have $\ddot u = G'(u)\dot u$ and
hence \(-\ddot u/\dot u = -G'(u)\): it is the slope of the transition function
at the current point, the speed that a linearisation performed at $u$ rather
than at zero would give. For the second-order approximation, this speed is
$\beta - 2\gamma u$, which is no longer constant: above $\beta$ below the
steady state, below it above. For the exact dynamics, the derivative of $G$
computed above gives:
\[
-\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}}
\]
of which $\beta-2\gamma u$ is the first-order expansion in $u$. This speed of
capital is not the speed of output $\beta(t)$ of the property established in
the section without approximation: since $g_{\hat y} = \alpha(\hat k)\,g_{\hat k}$,
the latter contains in addition the term due to the variation of the
capital share, which is the origin of the bracket. The two coincide when
$\alpha$ is constant, that is in the Cobb-Douglas case. The previous property
gives the counterpart of $\beta - 2\gamma u$ for output, $\beta - 2\gamma_y v$,
which compares directly with $\beta_{\sigma}$.
#+END_remarque

#+BEGIN_remarque
The approximate equation has a second steady state, $u = \beta/\gamma$,
unstable, which has no counterpart in the model: it is an artefact of the
approximation, as the figure on the [[#locale][local approximation]] led one to expect
for any polynomial approximation of the transition function. If
$u(0)>\beta/\gamma$, the denominator vanishes in finite time, and the
approximate solution blows up exactly like that of the [[ex-bernoulli][example]] treated
above. For a Cobb-Douglas technology, $\beta/\gamma = 2/(1-\alpha)$, that is
$3.125$ for \(\alpha=0.36\): one would have to start from more than
twenty-two times the steady-state capital, which is very far from the
neighbourhood where the approximation makes sense. The second-order
approximation can therefore be used without fear, but it cannot be used
anywhere.
#+END_remarque

#+BEGIN_remarque
In the Cobb-Douglas case, the exact model is itself a Bernoulli equation:
$\dot{\hat k} = s\hat k^{\alpha} - (n+g+\delta)\hat k$, with $m=\alpha$. The
change of variable $z = \hat k^{1-\alpha}$ makes it linear,
$\dot z = (1-\alpha)\left[s - (n+g+\delta)z\right]$, and one obtains the
explicit solution of the transition dynamics, used in the note on the
[[https://stephane-adjemian.fr/en/posts/simulating-the-solow-model/][simulation of the Solow model]]:
\[
u(t) = \frac{1}{1-\alpha}\log\left[1 + \left(e^{(1-\alpha)u(0)} - 1\right)e^{-\beta t}\right]
\]
Equivalently, the variable $w = e^{-(1-\alpha)u} = (\hat k/\hat k^{\star})^{-(1-\alpha)}$
follows the logistic dynamics $\dot w = \beta w(1-w)$. The second-order
approximation therefore replaces a Bernoulli equation in $\hat k$ by a
Bernoulli equation in $u$, and the first-order approximation by a linear
equation in $u$. This is why this technology allows the approximation errors
to be measured exactly.
#+END_remarque

The figure [[fig:ordre2][below]] compares the exact dynamics of the deviation $u$, obtained
from the previous formula in the Cobb-Douglas case and by numerical
integration in the case $\sigma=0.5$, with the first- and second-order
approximations, starting from $0.25\,\hat k^{\star}$ and from
$4\,\hat k^{\star}$.

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

  def order2(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-en :exports none :results none
  tg = np.linspace(0, 60, 400)
  fig, axes = plt.subplots(2, 2, figsize=(11, 8))
  errors = {}
  for j, sigma in enumerate((1.0, 0.5)):
      kstar = technology(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 = order1(tg, u0), order2(tg, u0, sigma)
          errors[(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='exact' if ratio < 1 else None)
          axes[0, j].plot(tg, u1, 'b', linestyle=style, linewidth=1, label='order one' if ratio < 1 else None)
          axes[0, j].plot(tg, u2, 'r', linestyle=style, linewidth=1, label='order two' if ratio < 1 else None)
          axes[1, j].plot(tg, u1-u, 'b', linestyle=style, linewidth=1,
                          label=r'order one, start at $%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'order two, start at $%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))
      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('error')
      axes[1, j].legend(fontsize=8)
      for ax in axes[:, j]:
          ax.set_xlabel(r'$t$ (years)')
  fig.tight_layout()
  fig.savefig("solow-ordre2.svg", transparent=True)
#+end_src

#+CAPTION: *Exact dynamics of the deviation from the steady state (black) and first-order (blue) and second-order (red) approximations, starting from $0.25\,\hat k^{\star}$ (solid) and from $4\,\hat k^{\star}$ (dashed), for a Cobb-Douglas technology (left) and a CES with $\sigma=0.5$ (right). Bottom, the approximation errors.*
#+LABEL: fig:ordre2
[[file:solow-ordre2.svg]]

The maximum error on the logarithmic deviation, in absolute value, is as
follows:

| start                    | technology           | order one | order two  |
|--------------------------+----------------------+-----------+------------|
| $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$    |

In the Cobb-Douglas case, the second order divides the first-order error by
more than five from below and by three from above, even though the starting
points retained are far from small: an initial capital four times below or
above the steady-state capital. The first-order approximation always errs in
the same direction, it underestimates the speed from below and overestimates
it from above, which is exactly the content of the property established in
the section [[#sans-approximation][without approximation]]; the second order, by restoring the
dependence of the speed on the position, corrects this asymmetry. An error
remains, larger from above than from below, which would only shrink by
pushing the expansion to the next order.

The case $\sigma=0.5$ holds a surprise: from below, the second order does
/worse/ than the first order. This is not a computational error, but a
reminder of what a local approximation is. The curvature of the transition
function has, by the proof of the property, the sign of
$1-\alpha(\hat k)/\sigma$. With $\sigma=0.5$, the capital share increases as
capital decreases, and exceeds $\sigma$ as soon as $\hat k$ falls below half
of \(\hat k^{\star}\): the function $G$ is convex in a neighbourhood of the
steady state, where the approximation is built, but concave on the first
part of the path, where it is used. The quadratic term, which reproduces the
curvature at the steady state, then corrects in the wrong direction. That
the first order does well is due to a compensation, the exact speed of
capital being sometimes below, sometimes above $\beta^{\star}$ along this
path, as the figure on the [[fig:vitesse][speed of convergence]] suggested. From above,
where the function remains convex, the second order divides the error by
four.

One should not conclude that the second-order approximation is always
preferable. It requires knowing the curvature of the transition function at
the steady state, here the elasticity of substitution, where the first order
makes do with the capital share; and its explicit solution is a privilege of
dimension one, which disappears as soon as one moves to a system. This is
why log-linearisation remains the reference method for dynamic general
equilibrium models, and why higher orders are only sought when the question
at hand requires them.