Solving initial value problems (IVPs) in quantecon

David R. Pugh

1. Introduction

This notebook demonstrates how to solve initial value problems (IVPs) using the quantecon Python library. Before demonstrating how one might solve an IVP using quantecon, I provide formal definitions for ordinary differential equations (ODEs) and initial value problems (IVPs), as well as a short discussion of finite-difference methods that will be used to solve IVPs.

Ordinary differential equations (ODE)

An ordinary differential equation (ODE) is in equation of the form

\begin{equation} \textbf{y}'= \textbf{f}(t ,\textbf{y}) \tag{1.1} \end{equation}

where $\textbf{f}:[t_0,\infty) \times \mathbb{R}^n\rightarrow\mathbb{R}^n$. In the case where $n=1$, then equation 1.1 reduces to a single ODE; when $n>1$, equation 1.1 defines a system of ODEs. ODEs are one of the most basic examples of functional equations: a solution to equation 1.1 is a function $\textbf{y}(t): D \subset \mathbb{R}\rightarrow\mathbb{R}^n$. There are potentially an infinite number of solutions to the ODE defined in equation 1.1. In order to reduce the number of potentially solutions, we need to impose a bit more structure on the problem.

Initial value problems (IVPs)

An initial value problem (IVP) has the form

\begin{equation} \textbf{y}'= \textbf{f}(t ,\textbf{y}),\ t \ge t_0,\ \textbf{y}(t_0) = \textbf{y}_0 \tag{1.2} \end{equation}

where $\textbf{f}:[t_0,\infty) \times \mathbb{R}^n\rightarrow\mathbb{R}^n$ and the initial condition $\textbf{y}_0 \in \mathbb{R}^n$ is a given vector. Alternatively, I could also specify an initial value problem by imposing a terminal condition of the form $\textbf{y}(T) = \textbf{y}_T$. The key point is that the solution $\textbf{y}(t)$ is pinned down at one $t\in[t_0, T]$.

The solution to the IVP defined by equation 1.2 is the function $\textbf{y}(t): [t_0,T] \subset \mathbb{R}\rightarrow\mathbb{R}^n$ that satisfies the initial condition $\textbf{y}(t_0) = \textbf{y}_0$. So long as the function $\textbf{f}$ is reasonably well-behaved, the solution $\textbf{y}(t)$ exists and is unique.

Finite-difference methods

Finite-difference methods are perhaps the most commonly used class of numerical methods for approximating solutions to IVPs. The basic idea behind all finite-difference methods is to construct a difference equation

\begin{equation} \frac{\textbf{y}(t_i + h) - \textbf{y}_i}{h} \approx \textbf{y}'(t_i) = \textbf{f}(t_i ,\textbf{y}(t_i)) \tag{1.3} \end{equation}

which is "similar" to the differential equation at some grid of values $t_0 < \dots < t_N$. Finite-difference methods then "solve" the original differential equation by finding for each $n=0,\dots,N$ a value $\textbf{y}_n$ that approximates the value of the solution $\textbf{y}(t_n)$.

It is important to note that finite-difference methods only approximate the solution $\textbf{y}$ at the $N$ grid points. In order to approximate $\textbf{y}$ between grid points one must resort to some form of interpolation.

The literature on finite-difference methods for solving IVPs is vast and there are many excellent reference texts. Those interested in a more in-depth treatment of these topics, including formal proofs of convergence, order, and stability of the numerical methods used in this notebook, should consult (Hairer, 1993), (Butcher, 2008), (Iserles, 2009). Chapter 10 of (Judd, 1998) covers a subset of these more formal texts with a specific focus on economic applications.

2. Examples

The remainder of this notebook demonstrates the usage and functionality of the quantecon.ivp module by way of example. To get started, we need to import the quantecon.ivp module...

2.1 Lotka-Volterra "Predator-Prey" model

We begin with the Lotka-Volterra model, also known as the predator-prey model, which is a pair of first order, non-linear, differential equations frequently used to describe the dynamics of biological systems in which two species interact, one a predator and the other its prey. The model was proposed independently by Alfred J. Lotka in 1925 and Vito Volterra in 1926.

\begin{align} \frac{du}{dt} =& au - buv \tag{2.1.1} \\ \frac{dv}{dt} =& -cv + dbuv \tag{2.1.2} \end{align}

where $u$ is the number of preys (for example, rabbits), $v$ is the number of predators (for example, foxes) and $a, b, c, d$ are constant parameters defining the behavior of the population.

Parameter definitions are as follows:

• $a$: the natural growing rate of prey in the absence of predators.
• $b$: the natural dying rate of prey due to predation.
• $c$: the natural dying rate of predators, in teh absence of prey.
• $d$: the factor describing how many caught prey is necessary to create a new predator.

I will use $\textbf{y}=[u, v]$ to describe the state of both populations.

2.1.1 Defining an instance of the IVP class

First, we need to create an instance of the IVP class representing the Lotka-Volterra "Predator-Prey" model. To initialize an instance of the IVP class we need to define the following...

• f : Callable of the form f(t, y, *f_args). The function f is the right hand side of the system of equations defining the model. The independent variable, t, should be a scalar; y is an ndarray of dependent variables with y.shape == (n,). The function f should return a scalar, ndarray or list (but not a tuple).
• jac : Callable of the form jac(t, y, *jac_args), optional(default=None). The Jacobian of the right hand side of the system of equations defining the ODE.

$$ \mathcal{J}_{i,j} = \bigg[\frac{\partial f_i}{\partial y_j}\bigg] \tag {2.1.3}$$

Most all of this information can be found in the docstring for the ivp.IVP class.

From the docstring we see that we are required to define a function describing the right-hand side of the system of differential equations that we wish to solve. While, optional, it is always a good idea to also define a function describing the Jacobian matrix of partial derivatives.

For the Lotka-Volterra model, these two functions would look as follows...

We can go ahead and create our instance of the ivp.IVP class representing the Lotka-Volterra model using the above defined functions as follows...

2.1.2 Defining model parameters

In order to simulate the model, however, we will need to supply values for the model parameters $a,b,c,d$. First, let's define some "reasonable" values.

In order to add these parameter values to our model we need to pass them as arguments to the set_f_params and set_jac_params methods of the newly created instance of the ivp.IVP class. Check the doctrings of the methods for information on the appropriate syntax...

From the docstring we see that both the set_f_params and the set_jac_params methods take an arbitrary number of positional arguments.

...and we can inspect that values of these attributes and see that the return results are the same.

Alternatively, we could just directly set the f_params and jac_params attributes without needing to explicitly call either the set_f_params and set_jac_params methods!

2.1.3 Using ivp.IVP.solve to integrate the ODE

In order to solve a system of ODEs, the ivp.IVP.solve method provides an interface into the ODE integration routines provided in the scipy.integrate.ode module. The method takes the following parameters...

• t0 : float. Initial condition for the independent variable.
• y0 : array_like (float, shape=(n,)). Initial condition for the dependent variables.
• h : float, optional(default=1.0). Step-size for computing the solution. Can be positive or negative depending on the desired direction of integration.
• T : int, optional(default=None). Terminal value for the independent variable. One of either T or g must be specified.
• g : Callable of the form g(t, y, args), optional(default=None). Provides a stopping condition for the integration. If specified user must also specify a stopping tolerance, tol.
• tol : float, optional (default=None). Stopping tolerance for the integration. Only required if g is also specifed.
• integrator : str, optional(default='dopri5') Must be one of 'vode', 'lsoda', 'dopri5', or 'dop853'
• step : bool, optional(default=False) Allows access to internal steps for those solvers that use adaptive step size routines. Currently only 'vode', 'zvode', and 'lsoda' support step=True.
• relax : bool, optional(default=False) Currently only 'vode', 'zvode', and 'lsoda' support relax=True.
• **kwargs : dict, optional(default=None). Dictionary of integrator specific keyword arguments. See the Notes section of the docstring for scipy.integrate.ode for a complete description of solver specific keyword arguments.

... and returns:

• solution: array_like (float). Simulated solution trajectory.

The above information can be found in the doctring of the ivp.IVP.solve method.

Example usage

Using dopri5, an embedded Runge-Kutta method of order 4(5) with adaptive step size control due to Dormand and Prince, integrate the model forward from an initial condition of 10 rabbits and 5 foxes for 15 years.

Plotting the solution

Once we have computed the solution, we can plot it using the excellent matplotlib Python library.

Note that we have plotted the time paths of rabbit and fox populations as sequences of points rather than smooth curves. This is done to visually emphasize the fact that finite-difference methods used to approximate the solution return a discrete approximation to the true continuous solution.

2.1.4 Using ivp.IVP.interpolate to interpolate the solution

The IVP.interpolate method provides an interface to the parametric B-spline interpolation routines in scipy.interpolate in order to constuct a continuous approximation to the true solution. For more details on B-spline interpolation, including some additional economic applications, see chapter 6 of (Judd, 1998). The ivp.IVP.interpolate method takes the following parameters...

• traj : array_like (float) Solution trajectory providing the data points for constructing the B-spline representation.
• ti : array_like (float) Array of values for the independent variable at which to interpolate the value of the B-spline.
• k : int, optional(default=3) Degree of the desired B-spline. Degree must satisfy :math:1 \le k \le 5.
• der : int, optional(default=0) The order of derivative of the spline to compute (must be less than or equal to k).
• ext : int, optional(default=2) Controls the value of returned elements for outside the original knot sequence provided by traj. For extrapolation, set ext=0; ext=1 returns zero; ext=2 raises a ValueError.

... and returns:

• interp_traj: ndarray (float). The interpolated trajectory.

Example usage

Approximate the solution to the Lotka-Volterra model at 1000 evenly spaced points using a 5th order B-spline interpolation and no extrapolation.

Plotting the interpolated solution

We can now plot the interpolated solution using matplotlib as follows...

Note that we have plotted the time paths of rabbit and fox populations as smooth curves. This is done to visually emphasize the fact that the B-spline interpolation methods used to approximate the solution return a continuous approximation to the true continuous solution.

2.1.5 Assessing accuracy using ivp.IVP.compute_residual

After computing a continuous approximation to the solution of our IVP, it is important to verify that the computed approximation is actually a "good" approximation. To assess the accuracy of our numerical solution we first define a residual function, $R(t)$, as the difference between the derivative of the B-spline approximation of the solution trajectory, $\hat{\textbf{y}}'(t)$, and the right-hand side of the original ODE evaluated along the approximated solution trajectory.

\begin{equation} \textbf{R}(t) = \hat{\textbf{y}}'(t) - f(t, \hat{\textbf{y}}(t)) \tag{2.1.4} \end{equation}

The idea is that if our numerical approximation of the true solution is "good", then this residual function should be roughly zero everywhere within the interval of approximation. The ivp.IVP.compute_residual method takes the following parameters...

• traj : array_like (float). Solution trajectory providing the data points for constructing the B-spline representation.
• ti : array_like (float). Array of values for the independent variable at which to interpolate the value of the B-spline.
• k : int, optional(default=3). Degree of the desired B-spline. Degree must satisfy $1 \le k \le 5$.
• ext : int, optional(default=2). Controls the value of returned elements for outside the original knot sequence provided by traj. For extrapolation, set ext=0; ext=1 returns zero; ext=2 raises a ValueError.

... and returns:

• residual : array (float) Difference between the derivative of the B-spline approximation of the solution trajectory and the right-hand side of the ODE evaluated along the approximated solution trajectory.

Remember to check the docstring for more information!

Example usage

Compute the residual to the Lotka-Volterra model at 1000 evenly spaced points using a 1st order B-spline interpolation (which is equivalent to linear interpolation!).

Plotting the residual

In your introductory econometrics/statistics course, your professor likely implored you to "always plot your residuals!" This maxim of data analysis is no less true in numerical analysis. However, while patterns in residuals are generally a "bad" thing in econometrics/statistics (as they suggest model mispecification, or other related problems), patterns in a residual function, $\textbf{R}(t)$, in numerical analysis are generally OK (and in certain cases actually desireable!).

In this context, what is important is that overall size of the residuals is "small".

Understanding determinants of accuracy

We can use IPython widgets to investigate the determinants of accuracy of our approximated solution. Good candidates for exploration are...

• h: the step size used in computing the initial finite difference solution.
• atol: the absolute tolerance for the solver.
• rtol: the relative tolerance for the solver.
• k: the degree of the B-spline used in the interpolation of that finite difference solution.

Now we can make use of the @interact decorator and the various IPython widgets to create an interactive visualization of the residual plot for the Lotka-Volterra "Predator-Prey" model.

Sensitivity to parameters

Once we have computed and plotted an approximate solution (and verified that the approximation is a good one by plotting the residual function!), we can try and learn something about the dependence of the solution on model parameters.

2.2 The Lorenz equations

The Lorenz equations are a system of three coupled, first-order, non-linear differential equations which describe the trajectory of a particle through time. The system was originally derived by as a model of atmospheric convection, but the deceptive simplicity of the equations have made them an often-used example in fields beyond atmospheric physics.

The equations describe the evolution of the spatial variables x, y, and z, given the governing parameters $\sigma, \beta, \rho$, through the specification of the time-derivatives of the spatial variables:

\begin{align} \frac{dx}{dt} =& \sigma(y − x) \tag{2.2.1} \\ \frac{dy}{dt} =& x(\rho − z) − y \tag{2.2.2} \\ \frac{dz}{dt} =& xy − \beta z \tag{2.2.3} \end{align}

The resulting dynamics are entirely deterministic giving a starting point $(x_0,y_0,z_0)$. Though it looks straightforward, for certain choices of the parameters, the trajectories become chaotic, and the resulting trajectories display some surprising properties.

2.2.1 Incorporating SymPy

While deriving the Jacobian matrix by hand is trivial for most simple 2D or 3D systems, it can quickly become tedious and error prone for larger systems (or even some highly non-linear 2D systems!). In addition to being an important input to most ODE integrators/solvers, Jacobians are also useful for assessing the stability properties of equilibria.

An alternative approach for solving IVPs using the ivp module, which leverages the SymPy Python library to do the tedious computations involved in deriving the Jacobian, is as follows.

1. Define the IVP using SymPy.
2. Use SymPy routines for computing the Jacobian.
3. Wrap the symbolic expressions as callable NumPy functions.
4. Use these functions to create an instance of the ivp.IVP class.

The remainder of this notebook implements each of these steps to solve and analyze the Lorenz equations defined above.

Step 1: Defining the Lorenz equations using SymPy

We begin by defining a sp.Matrix instance containing the three Lorenz equations...

Let's take a check out our newly defined _lorenz_system and make sure it looks as expected...

Step 2: Computing the Jacobian using SymPy

Once we have defined our model as a SymPy matrix, computing the Jacobian is trivial...

Step 3: Wrap the SymPy expression to create vectorized NumPy functions

Now we wrap the SymPy matrices defining the model and the Jacobian to create vectorized NumPy functions. It is crucial that the interface for our wrapped functions matches the interface required by the f and jac parameters which we will pass to the ivp.IVP constructor to create an instance of the ivp.IVP class representing the Lorenz equations.

Recall from the ivp.IVP docstring that...

f : callable `f(t, y, *f_args)`
    Right hand side of the system of equations defining the ODE. The
    independent variable, `t`, is a `scalar`; `y` is an `ndarray`
    of dependent variables with `y.shape == (n,)`. The function `f`
    should return a `scalar`, `ndarray` or `list` (but not a
    `tuple`).
jac : callable `jac(t, y, *jac_args)`, optional(default=None)
    Jacobian of the right hand side of the system of equations defining
    the ODE.

    .. :math:

        \mathcal{J}_{i,j} = \bigg[\frac{\partial f_i}{\partial y_j}\bigg]

Thus our wrapped functions need to take a float t as the first argument, and array y as a second argument, followed by some arbitrary number of model parameters. We can handle all of this as follows...

Step 4: Use these functions to create an instance of the IVP class

First we define functions describing the right-hand side of the ODE and the Jacobian which we need to initialize the ivp.IVP class...

... next we define a tuple of model parameters...

... finally, we are ready to create the instance of the ivp.IVP class representing the Lorenz equations.

2.2.2 Solving the Lorenz equations

At this point I proceed in exactly the same fashion as in the previous Lotka-Volterra equations example:

1. Solve the model using a discretized, finite-difference approximation.
2. Use the discretized approximation in conjunction with parametric B-spline interpolation to construct a continuous approximation of the true solution.
3. Compute and analyze the residual of the approximate solution.

Step 1. Solve the model using a discretized, finite-difference approximation

Using dop853, an embedded Runge-Kutta method of order 7(8) with adaptive step size control due to Dormand and Prince, integrate the Lorenz equations forward from an initial condition of $X_0 = (1.0, 1.0, 1.0)$ from $t=0$, to $T=100$.

Plotting the solution time paths

We can use IPython widgets to construct a "poor man's" animation of the evolution of the Lorenz equations.

Step 2: Construct a continuous approximation to the solution.

Let's construct a continuous approximation to the solution of the Lorenz equations at 10000 evenly spaced points using a 5th order B-spline interpolation.

Plotting 2D projections of the solution in phase space

The underlying structure of the Lorenz system becomes more apparent when I plot the interpolated solution trajectories in phase space. We can plot 2D projections of the time paths of the solution in phase space...

Step 3: Compute the residuals to assese the accuracy of our solution

Finally, to assess the accuracy of our solution we need to compute and plot the solution residuals at 10000 evenly spaced points using a 5th order B-spline interpolation.

Again, we want to confirm that the residuals are "small" everywhere. Patterns, if they exists, are not cause for concern.