API

riccati.chebutils module

riccati.chebutils.cheb(n)

Returns a differentiation matrix D of size (n+1, n+1) and (n+1) Chebyshev nodes x for the standard 1D interval [-1, 1]. The matrix multiplies a vector of function values at these nodes to give an approximation to the vector of derivative values. Nodes are output in descending order from 1 to -1. The nodes are given by

\[x_p = \cos \left( \frac{\pi n}{p} \right), \quad n = 0, 1, \ldots p.\]

Parameters:

n: int

Number of Chebyshev nodes - 1.

Returns:

D: numpy.ndarray [float]

Array of size (n+1, n+1) specifying the differentiation matrix.

x: numpy.ndarray [float]

Array of size (n+1,) containing the Chebyshev nodes.

riccati.chebutils.coeffs2vals(coeffs)

Convert the Chebyshev coefficient representation of a set of polynomials P_j to their values at Chebyshev nodes of the second kind (ordered from +1 to -1). This function returns a matrix V such that for an input coefficient matrix C,

\[V_{ij} = P_j(x_i) = \sum_{k=0}^{n} C_{kj}T_k(x_i).\]

Taken from the coeff2vals function in the chebfun package.

Parameters:

coeffs: numpy.ndarray [float (real)]

An array of size (n+1, m), with the (i, j)th element representing the projection of the jth input polynomial onto the ith Chebyshev polynomial.

Returns:

values: np.ndarray [float (real)]

An array of size (n+1, m), with the (i, j)th element representing the jth input polynomial evaluated at the ith Chebyshev node.

riccati.chebutils.integrationm(n)

Chebyshev integration matrix. It maps function values at n Chebyshev nodes of the second kind, ordered from +1 to -1, to values of the integral of the interpolating polynomial at those points, with the last value (start of the interval) being zero. Taken from the cumsummat function in the chebfun package.

Parameters:

n: int

Number of Chebyshev nodes the integrand is evaluated at. Note the nodes are always ordered from +1 to -1.

Returns:

Q: np.ndarray [float (real)]

Integration matrix of size (n, n) that maps values of the integrand at the n Chebyshev nodes to values of the definite integral on the interval, up to each of the Chebyshev nodes (the last value being zero by definition).

riccati.chebutils.interp(s, t)

Creates interpolation matrix from an array of source nodes s to target nodes t. Taken from here .

Parameters:

s: numpy.ndarray [float]

Array specifying the source nodes, at which the function values are known.

t: numpy.ndarray [float]

Array specifying the target nodes, at which the function values are to be interpolated.

Returns:

L: numpy.ndarray [float]

Array defining the inteprolation matrix L, which takes function values at the source points s and yields the function evaluated at target points t. If s has size (p,) and t has size (q,), then L has size (q, p).

riccati.chebutils.quadwts(n)

Clenshaw-Curtis quadrature weights mapping function evaluations at (n+1) Chebyshev nodes of the second kind, ordered from +1 to -1, to value of the definite integral of the interpolating function on the same interval. Taken from [1] Ch 12, clencurt.m

Parameters:

n: int

The (number of Chebyshev nodes - 1) for which the quadrature weights are to be computed.

Returns:

w: np.ndarray [float (real)]

Array of size (n+1,) containing the quadrature weights.

References

[1]

Trefethen, Lloyd N. Spectral methods in MATLAB. Society for industrial and applied mathematics, 2000.

riccati.chebutils.spectral_cheb(info, x0, h, y0, dy0, niter)

Utility function to apply a spectral collocation method based on Chebyshev nodes from x = x0 to x = x0+h, starting from the initial conditions y(x0) = y0, y’(x0) = dy0. In each spectral collocation/Chebyshev step, the solver iterates over how many nodes to perform the step with, starting from info.nini, doubling in every iteration until info.nmax is reached or the required tolerance is satisfied, whichever comes sooner. The niter parameter keeps track of how many iterations have been done and is used to retrieve the relevant pre-computed differentiation matrix and vector of nodes from the input info object.

Parameters:

info: Solverinfo object

x0: float (real)

Value of the independent variable to take a spectral step from.

h: float (real)

Stepsize.

y0, dy0: complex

Initial conditions (value of the solution and its derivative) at x0.

niter: int

Integer to keep track of how many iterations of the spectral collocation step have been done.

Returns:

y1, dy1: complex

Numerical estimate of the solution and its derivative at the end of the step, at x0+h.

xscaled: numpy.ndarray [float]

Chebyshev nodes used for the current iteration of the spectral collocation method scaled to lie between [x0, x0+h].

riccati.chebutils.vals2coeffs(values)

Convert a matrix of values of m polynomials evaluated at n+1 Chebyshev nodes of the second kind, ordered from +1 to -1, to their interpolating Chebyshev coefficients. This function returns the coefficient matrix C such that for an input matrix of values V,

\[F_j(x) = \sum_{k=0}^{n} C_{kj}T_k(x) \]

interpolates the values \([V_{0j}, V_{1j}, \ldots, V_{nj}]\) for \(j = 0..(m-1)\).

Taken from the vals2coeffs function in the chebfun package.

Parameters:

values: np.ndarray [float (real)]

An array of size (n+1, m), with the (i, j)th element being the jth polynomial evaluated at the ith Chebyshev node.

Returns:

coeffs: np.ndarray [float (real)]

An array of size (n+1, m) with the (i, j)th element being the coefficient multiplying the ith Chebyshev polynomial for interpolating the jth input polynomial.

riccati.solversetup module

class riccati.solversetup.Solverinfo(w, g, h, nini, nmax, n, p)

Bases: object

Class to store information of an ODE solve run. When initialized, it computes Chebyshev nodes and differentiation matrices neede for the run once and for all.

Attributes:

w, g: callable

Frequency and friction function associated with the ODE, defined as \(\omega\) and \(\gamma\) in \(u'' + 2\gamma u' + \omega^2 u = 0\).

h0: float

Initial interval length which will be used to estimate the initial derivatives of w, g.

nini, nmax: int

Minimum and maximum (number of Chebyshev nodes - 1) to use inside Chebyshev collocation steps. The step will use nmax nodes or the minimum number of nodes necessary to achieve the required local error, whichever is smaller. If nmax > 2`nini`, collocation steps will be attempted with \(2^i\) nini nodes at the ith iteration.

n: int

(Number of Chebyshev nodes - 1) to use for computing Riccati steps.

p: int

(Number of Chebyshev nodes - 1) to use for estimating Riccati stepsizes.

denseout: bool

Defines whether or not dense output is required for the current run.

h: float

Current stepsize.

y: np.ndarray [complex]

Current state vector of size (2,), containing the numerical solution and its derivative.

wn, gn: np.ndarray [complex]

The frequency and friction function evaluated at n + 1 Chebyshev nodes over the interval [x, x + h ], where x is the value of the independent variable at the start of the current step and h is the current stepsize.

Ds: list [np.ndarray [float]]

List containing differentiation matrices of sizes \((2^i n_{\mathrm{ini}} + 1,2^i n_{\mathrm{ini}} + 1)\) for \(i = 0, 1, \ldots, \lfloor \log_2\frac{n_{\mathrm{max}}}{n_{\mathrm{ini}}} \rfloor.\)

nodes: list [np.ndarray [float]]

List containing vectors of Chebyshev nodes over the standard interval, ordered from +1 to -1, of sizes \((2^i n_{\mathrm{ini}} + 1,)\) for \(i = 0, 1, \ldots, \lfloor \log_2 \frac{n_{\mathrm{max}}}{n_{\mathrm{ini}}} \rfloor.\)

ns: np.ndarray [int]

Vector of lengths of node vectors stored in nodes, i.e. the integers \(2^i n_{\mathrm{ini}} + 1\) for \(i = 0, 1, \ldots, \lfloor \log_2 \frac{n_{\mathrm{max}}}{n_{\mathrm{ini}}} \rfloor.\)

xn: np.ndarray [float]

Values of the independent variable evaluated at (n + 1) Chebyshev nodes over the interval [x, x + h], where x is the value of the independent variable at the start of the current step and h is the current stepsize.

xp: np.ndarray [float]

Values of the independent variable evaluated at (p + 1) Chebyshev nodes over the interval [x, x + h], where x is the value of the independent variable at the start of the current step and h is the current stepsize.

xpinterp: np.ndarray [float]

Values of the independent variable evaluated at p points over the interval [x, x + h] lying in between Chebyshev nodes, where x is the value of the independent variable at the start of the current step and h is the current stepsize. The in-between points \(\tilde{x}_p\) are defined by

L: np.ndarray [float]

Interpolation matrix of size (p`+1, `p), used for interpolating a function between the nodes xp and xpinterp (for computing Riccati stepsizes).

quadwts: np.ndarray [float]

Vector of size (n + 1,) containing Clenshaw-Curtis quadrature weights.

n_chebnodes: int

Number of times Chebyhev nodes have been calculated, i.e. riccati.chebutils.cheb has been called.

n_chebstep: int

Number of Chebyshev steps attempted.

n_chebits: int

Number of times an iteration of the Chebyshev-grid-based collocation method has been performed (note that if nmax >= 4`nini` then a single Chebyshev step may include multiple iterations!).

n_LS: int

Number of times a linear system has been solved.

n_riccstep: int

Number of Riccati steps attempted.

Methods

increase:	Increase the various counters (attributes starting with n_) by given values.
output:	Return the state of the counters as a dictionary.

increase(chebnodes=0, chebstep=0, chebits=0, LS=0, riccstep=0)

Increases the relevant attribute of the class (a counter for a specific arithmetic operation) by a given number. Used for generating performance statistics.

Parameters:

chebnodes: int

Count by which to increase riccati.solversetup.Solverinfo.n_chebnodes.

chebstep: int

Count by which to increase riccati.solversetup.Solverinfo.n_chebstep.

chebits: int

Count by which to increase riccati.solversetup.Solverinfo.n_chebits.

LS: int

Count by which to increase riccati.solversetup.Solverinfo.n_LS.

riccstep: int

Count by which to increase riccati.solversetup.Solverinfo.n_riccstep.

Returns:

None

output(steptypes)

Creates a dictionary of the counter-like attributes of Solverinfo, namely: the number of attempted Chebyshev steps n_chebsteps, the number of Chebyshev iterations n_chebits, the number of attempted Riccati steps n_riccstep, the number of linear solver n_LS, and the number of times Chebyshev nodes have been computed, n_chebnodes. It also logs the number of steps broken down by steptype: in the relevant fields, (n, m) means there were n attempted steps out of which m were successful.

Parameters:

steptypes: list [int]

List of steptypes (of successful steps) produced by riccati.evolve.solve(), each element being 0 (Chebyshev step) or 1 (Riccati step).

Returns:

statdict: dict

Dictionary with the following keywords:

cheb steps: tuple [int]

(n, m), where n is the total number of attempted Chebyshev steps out of which m were successful.

cheb iterations: int

Total number of iterations of the Chebyshev collocation method.

ricc steps: tuple [int]

(n, m), where n is the total number of attempted Chebyshev steps out of which m were successful.

linear solves: int

Total number of times a linear solve has been performed.

cheb nodes: int

Total number of times a call to compute Chebyshev nodes has been made.

riccati.solversetup.solversetup(w, g, h0=0.1, nini=16, nmax=32, n=16, p=16)

Sets up the solver by generating differentiation matrices based on an increasing number of Chebyshev gridpoints (see riccati.solversetup.Solverinfo). Needs to be called before the first time the solver is ran or if nini or nmax are changed.

Parameters:

w: callable(t)

Frequency function in y’’ + 2g(t)y’ + w^2(t)y = 0.

g: callable(t)

Damping function in y’’ + 2g(t)y’ + w^2(t)y = 0.

h0: float

Initial stepsize.

nini: int

Lowest (number of nodes - 1) to start the Chebyshev spectral steps with.

nmax: int

Maximum (number of nodes - 1) for the spectral Chebyshev steps to try before decreasing the stepsize and reattempting the step with (nini + 1) nodes.

n: int

(Fixed) number of Chebyshev nodes to use for interpolation, integration, and differentiation in Riccati steps.

p: int

(Fixed) number of Chebyshev nodes to use for interpolation when determining the stepsize for Riccati steps.

Returns:

info: Solverinfo object

Solverinfo object created with attributes set as per the input parameters.

riccati.stepsize module

riccati.stepsize.choose_nonosc_stepsize(info, x0, h, epsh=0.2)

Chooses the stepsize for spectral Chebyshev steps, based on the variation of 1/w, the approximate timescale over which the solution changes. If over the suggested interval h 1/w changes by a fraction of \(\pm``epsh\) or more, the interval is halved, otherwise it’s accepted.

Parameters:

info: Solverinfo object

Solverinfo object which is used to retrieve Solverinfo.xp, the (p+1) Chebyshev nodes used for interpolation to determine the stepsize.

x0: float

Current value of the independent variable.

h: float

Initial estimate of the stepsize.

epsh: float

Tolerance parameter defining how much 1/w(x) is allowed to change over the course of the step.

Returns:

h: float

Refined stepsize over which 1/w(x) does not change by more than epsh/w(x).

riccati.stepsize.choose_osc_stepsize(info, x0, h, epsh=1e-12)

Chooses the stepsize h over which the functions w(x), g(x) can be represented sufficiently accurately. Evaluations of w(x) and g(x) at (p+1) Chebyshev nodes (given by info.xp) between [x0, x0+h] are used to infer w(x), g(x) values at p points halfway between the (p+1) nodes on the half circle,

\[x_{p,\mathrm{interp}} = \cos\left( \frac{\pi(2n+1)}{2p} \right), \quad n = 0, 1, \ldots, p-1.\]

These target nodes given by info.xpinterp. The same (p+1) points will then be used to compute a Riccati step. The interpolated values are then compared to the actual values of w(x), g(x) at info.xpinterp. If the largest relative error in w, g exceeds epsh, h is halved.

Parameters:

info: Solverinfo object

Solverinfo object used for two purposes: info.xpinterp, info.xp are looked up to determine the source and target points to determine whether Chebyshev interpolation of w(x), g(x) is accurate enough over the step; and info.wn, info.gn may be used to read off the w, g values at the source points.

x0: float

Current value of the independent variable.

h: float

Initial estimate of the stepsize.

epsh: float

Tolerance parameter defining the maximum relative error Chebyshev interpolation of w, g is allowed to have over the course of the proposed step [x0, x0+h].

Returns:

h: float

Refined stepsize over which Chebyshev interpolation of w, g has a relative error no larger than epsh.

riccati.step module

riccati.step.nonosc_step(info, x0, h, y0, dy0, epsres=1e-12)

A single Chebyshev step to be called from the solve() function. Advances the solution from x0 by h, starting from the initial conditions y(x0) = y0, y’(x0) = dy0. The function uses a Chebyshev spectral method with an adaptive number of nodes. Initially, info.nini nodes are used, which is doubled in each iteration until epsres relative accuracy is reached or the number of nodes would exceed info.nmax. The relative error is measured as the difference between the predicted value of the dependent variable at the end of the step obtained in the current iteration and in the previous iteration (with half as many nodes). If the desired relative accuracy cannot be reached with info.nmax nodes, it is advised to decrease the stepsize h, increase info.nmax, or use a different approach.

Parameters:

info: `Solverinfo` object

Solverinfo object used to read off various matrices required for numerical differentiation, and info.wn, info.gn for evaluations of w(x), g(x) over [x0, x0+h].

x0: float

Starting value of the independent variable.

h: float

Stepsize.

y0, dy0: complex

Value of the dependent variable and its derivative at x = x0.

epsres: float

Tolerance for the relative accuracy of Chebyshev steps.

Returns:

y[0], dy[0]: complex

Value of the dependent variable and its derivative at the end of the step, x = x0 + h.

maxerr: float

(Absolute) value of the relative difference of the dependent variable at the end of the step as predicted in the last and the previous iteration.

success: int

Takes the value 1 if the asymptotic series has reached epsres residual, 0 otherwise.

riccati.step.osc_step(info, x0, h, y0, dy0, epsres=1e-12, plotting=False, k=0)

A single Riccati step to be called from within the solve() function. Advances the solution from x0 by h, starting from the initial conditions y(x0) = y0, y’(x0) = dy0. The function will increase the order of the asymptotic series used for the Riccati equation until a residual of epsres is reached or the residual stops decreasing. In the latter case, the asymptotic series cannot approximate the solution of the Riccati equation with the required accuracy over the given interval; the interval (h) should be reduced or another approximation should be used instead.

Parameters:

info: `Solverinfo` object

Solverinfo object used to read off various matrices required for numerical differentiation, and info.wn, info.gn for evaluations of w(x), g(x) over [x0, x0+h].

x0: float

Starting value of the independent variable.

h: float

Stepsize.

y0, dy0: complex

Value of the dependent variable and its derivative at x = x0.

epsres: float

Tolerance for the relative accuracy of Riccati steps.

Returns:

y1[0], dy1[0]: complex

Value of the dependent variable and its derivative at the end of the step, x = x0 + h.

maxerr: float

Maximum value of the residual (after the final iteration of the asymptotic approximation) over the Chebyshev nodes across the interval.

success: int

Takes the value 1 if the asymptotic series has reached epsres residual, 0 otherwise.

phase: complex

Total phase change (not mod :math: 2pi!) of the dependent variable over the step.

Warning

This function relies on info.wn, info.gn being set correctly, as appropriate for a step of size h. If solve() is calling this funciont, that is automatically taken care of, but otherwise needs to be done manually.

riccati.evolve module

riccati.evolve.nonosc_evolve(info, x0, x1, h, y0, epsres=1e-12, epsh=0.2)

Allows continuous evolution between the independent variable values x0 and x1. Starting from x0 and y = [y0, yd0] takes a Chebyshev step of size h or until x1 is reached, and updates the attributes of info such that it can be called inside a while loop.

Parameters:

info: `Solverinfo` object

Solverinfo object used to read off info.Dn for differentiation, and info.wn, info.gn for evaluations of w(x), g(x) over [x0, x0+h].

x0: float

Starting value of the independent variable.

x1: float

Maximum value of the independent variable. This may not be reached in a single step (i.e. a single call of osc_evolve), but will not be exceeded.

h: float

Stepsize. If the step would result in the independent variable exceeding x1, this will be adjusted.

y0: np.ndarray [complex]

Value of the state vector at x0.

epsres: float

Tolerance for the relative accuracy of Riccati steps.

epsh: float

Tolerance for choosing the stepsize for Riccati steps.

Returns:

success: int

0 if the step failed, 1 if successful.

Warning

The user will need to set info.wn, info.gn correctly, as appropriate for a step of size h, before this function is called. This can be done by e.g. calling choose_osc_stepsize.

Notes

Note that epsh is defined differently for Riccati and Chebyshev steps; the same value may therefore not be appropriate for both osc_evolve() and nonosc_evolve().

riccati.evolve.osc_evolve(info, x0, x1, h, y0, epsres=1e-12, epsh=1e-12)

Allows continuous evolution between the independent variable values x0 and x1. Starting from x0 and y = [y0, yd0], this function takes a Riccati step of size h or until x1 is reached, and updates the attributes of info such that it can be called inside a while loop.

Parameters:

info: `Solverinfo` object

Solverinfo object used to read off info.Dn for differentiation, and info.wn, info.gn for evaluations of w(x), g(x) over [x0, x0+h].

x0: float

Starting value of the independent variable.

x1: float

Maximum value of the independent variable. This may not be reached in a single step (i.e. a single call of osc_evolve), but will not be exceeded.

h: float

Stepsize. If the step would result in the independent variable exceeding x1, this will be adjusted.

y0: np.ndarray [complex]

Value of the state vector at x0.

epsres: float

Tolerance for the relative accuracy of Riccati steps.

epsh: float

Tolerance for choosing the stepsize for Riccati steps.

Returns:

success: int

0 if the step failed, 1 if successful.

Warning

The user will need to set info.wn, info.gn correctly, as appropriate for a step of size h, before this function is called. This can be done by e.g. calling choose_osc_stepsize.

riccati.evolve.solve(info, xi, xf, yi, dyi, eps=1e-12, epsh=1e-12, xeval=array([], dtype=float64), hard_stop=False, warn=False)

Solves y’’ + 2gy’ + w^2y = 0 on the interval (xi, xf), starting from the initial conditions y(xi) = yi, y’(xi) = dyi. Keeps the residual of the ODE below eps, and returns an interpolated solution (dense output) at points xeval.

Parameters:

info: Solverinfo object

Objects containing differentiation matrices, etc.

xi, xf: float

Solution range.

yi, dyi: complex

Initial conditions, value of the dependent variable and its derivative at xi.

eps: float

Relative tolerance for the local error of both Riccati and Chebyshev type steps.

epsh: float

Relative tolerance for choosing the stepsize of Riccati steps.

xeval: list

List of x-values where the solution is to be interpolated (dense output) and returned.

hard_stop: bool

Whether to force the solver to have a potentially smaller last stepsize, in order to stop exactly at xf (rather than allowing the solver to step over it and get the value of the solution by interpolation).

warn: bool

Whether to display warnings, e.g. RuntimeWarnings, during a run. Due to the high level of adaptivity in this algorithm, it may throw several RuntimeWarnings even in a standard setup as it chooses the type of step, stepsize, and other parameters. For this reason, all warnings are silenced by default (warn = False). Set to True if you wish to see the warnings.

Returns:

xs: list [float]

Values of the independent variable at the internal steps of the solver.

ys, dys: list [complex]

Values of the dependent variable and its derivative at the internal steps of the solver.

successes: list [int]

Has elements 1 and 0: 1 denoting each successful step, 0 denoting unsuccessful steps.

phases: list [complex]

Complex phase of the solution accumulated during each successful Riccati step.

steptypes: list [int]

Types of successful steps taken: 1 for Riccati and 0 for Chebyshev.

yeval: numpy.array [complex]

Dense output, i.e. values of the solution at the requested independent-variable values specified in xeval. If xeval was not given, then it is an empty numpy array of shape (0,).