First-order differential equations

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:

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

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

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

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.

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

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

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

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.

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.

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

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.

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.

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

\(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 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}\).

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

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

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)

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.

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)

The figure 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.

As in the section on the 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 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 below illustrates this replacement, for the transition function of the Solow model which we shall meet again in the next section.

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.

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 the Solow model; the aim is simply to apply the methods presented above.

We take a neoclassical production function², \(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³. 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 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\).

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)