r"""
Base class for solving initial value problems (IVPs) of the form:
.. math::
\frac{dy}{dt} = f(t,y),\ y(t_0) = y_0
using finite difference methods. The `quantecon.ivp` class uses various
integrators from the `scipy.integrate.ode` module to perform the
integration (i.e., solve the ODE) and parametric B-spline interpolation
from `scipy.interpolate` to approximate the value of the solution
between grid points. The `quantecon.ivp` module also provides a method
for computing the residual of the solution which can be used for
assessing the overall accuracy of the approximated solution.
"""
import numpy as np
from scipy import integrate, interpolate
[docs]class IVP(integrate.ode):
r"""
Creates an instance of the IVP class.
Parameters
----------
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]
"""
def __init__(self, f, jac=None):
super(IVP, self).__init__(f, jac)
def _integrate_fixed_trajectory(self, h, T, step, relax):
"""Generates a solution trajectory of fixed length."""
# initialize the solution using initial condition
solution = np.hstack((self.t, self.y))
while self.successful():
self.integrate(self.t + h, step, relax)
current_step = np.hstack((self.t, self.y))
solution = np.vstack((solution, current_step))
if (h > 0) and (self.t >= T):
break
elif (h < 0) and (self.t <= T):
break
else:
continue
return solution
def _integrate_variable_trajectory(self, h, g, tol, step, relax):
"""Generates a solution trajectory of variable length."""
# initialize the solution using initial condition
solution = np.hstack((self.t, self.y))
while self.successful():
self.integrate(self.t + h, step, relax)
current_step = np.hstack((self.t, self.y))
solution = np.vstack((solution, current_step))
if g(self.t, self.y, *self.f_params) < tol:
break
else:
continue
return solution
def _initialize_integrator(self, t0, y0, integrator, **kwargs):
"""Initializes the integrator prior to integration."""
# set the initial condition
self.set_initial_value(y0, t0)
# select the integrator
self.set_integrator(integrator, **kwargs)
[docs] def compute_residual(self, traj, ti, k=3, ext=2):
r"""
The residual is the difference between the derivative of the B-spline
approximation of the solution trajectory and the right-hand side of the
original ODE evaluated along the approximated solution trajectory.
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`.
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`.
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.
"""
# B-spline approximations of the solution and its derivative
soln = self.interpolate(traj, ti, k, 0, ext)
deriv = self.interpolate(traj, ti, k, 1, ext)
# rhs of ode evaluated along approximate solution
T = ti.size
rhs_ode = np.vstack([self.f(ti[i], soln[i, 1:], *self.f_params)
for i in range(T)])
rhs_ode = np.hstack((ti[:, np.newaxis], rhs_ode))
# should be roughly zero everywhere (if approximation is any good!)
residual = deriv - rhs_ode
return residual
[docs] def solve(self, t0, y0, h=1.0, T=None, g=None, tol=None,
integrator='dopri5', step=False, relax=False, **kwargs):
r"""
Solve the IVP by integrating the ODE given some initial condition.
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 ``g(t, y, f_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.
Returns
-------
solution: ndarray (float)
Simulated solution trajectory.
"""
self._initialize_integrator(t0, y0, integrator, **kwargs)
if (g is not None) and (tol is not None):
soln = self._integrate_variable_trajectory(h, g, tol, step, relax)
elif T is not None:
soln = self._integrate_fixed_trajectory(h, T, step, relax)
else:
mesg = "Either both 'g' and 'tol', or 'T' must be specified."
raise ValueError(mesg)
return soln
[docs] def interpolate(self, traj, ti, k=3, der=0, ext=2):
r"""
Parametric B-spline interpolation in N-dimensions.
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`.
Returns
-------
interp_traj: ndarray (float)
The interpolated trajectory.
"""
# array of parameter values
u = traj[:, 0]
# build list of input arrays
n = traj.shape[1]
x = [traj[:, i] for i in range(1, n)]
# construct the B-spline representation (s=0 forces interpolation!)
tck, t = interpolate.splprep(x, u=u, k=k, s=0)
# evaluate the B-spline (returns a list)
out = interpolate.splev(ti, tck, der, ext)
# convert to a 2D array
interp_traj = np.hstack((ti[:, np.newaxis], np.array(out).T))
return interp_traj