Numerical Analysis Modules in gvar

gvar.GVars can be used in many numerical algorithms, to propagates errors through the algorithm. A code that is written in pure Python is likely to work well with gvar.GVars, perhaps with minor modifications. Here we describe some sample numerical codes, included in gvar, that have been adapted to work with gvar.GVars, as well as with floats. More examples will follow with time.

The sub-modules included here are:

Cubic Splines

Module gvar.cspline implements a class for smoothing and/or interpolating one-dimensional data using cubic splines:

class gvar.cspline.CSpline(xknot, yknot, deriv=(None, None), extrap_order=3, alg='steffen', warn=False)

Cubic spline approximation to a function.

Given N values of a function yknot[i] at N points xknot[i] for i=0..N-1 (the knots of the spline), the code

from gvar.cspline import CSpline

f = CSpline(xknot, yknot)

defines a function f such that: a) f(xknot[i]) = yknot[i] for all i; and b) f(x) and its first derivative are continuous. Function f(x) is a cubic polynomial between the knots xknot[i] (determined from the values and first derivatives at the adjacent knots). Argument x in f(x) is a number or an array of numbers (any shape). Numbers can be replaced by gvar.GVar objects in x, xknot and/or yknot.

Derivatives and integrals of the spline function are also available:

f.D(x) — first derivative at x;

f.D2(x) — second derivative at x;

f.D3(x) — third derivative at x;

f.integ(x) — integral from xknot[0] to x.

The derivatives are increasingly unreliable with increasing order.

Splines can be used outside the range covered by the defining xknot values. As this is often a bad idea, keyword parameter warn can be set to True so that warnings are generated if the spline is used out of range. This keyword is set in the original constructor (CSpline(x, y, warn=True)) or in calls to a CSpline object (eg, f.D(x, warn=True)). The spline value for an out-of-range point is calculated using a polynomial whose value and derivatives match those of the spline at the knot closest to the out-of-range point. The extrapolation polynomial is cubic by default, but lower orders are specified by setting parameter extrap_order to a (non-negative) integer less than 3.

The first derivatives of f(x) at the knots are determined from the function values in the vicinity of the knot. The derivatives at the endpoints (xknot[0] and xknot[-1] if the knots are ordered) can be specified explicitly, if they are known:

f = CSpline(xknot, yknot, deriv=[dydx_i, dydx_f])

where dydx_i is the derivative at xknot[0] and dydx_f is the derivative at xknot[-1].

There are different types of spline available, selected with keyword alg:

alg='steffen'

Monotonic cubic spline that has quadratic precision where the function is monotonic between knots. This is the default algorithm.

alg='cspline'

Cubic spline with continuous second derivatives f''(x) (in addition to f'(x) and f(x)). Sets f''(x) to zero at the endpoints (natural boundary conditions) if the first derivatives have not been specified there (using keyword deriv).

alg='pchip'

Monotonic cubic spline used in the scipy function scipy.interpolate.PchipInterpolator. This algorithm is usually less accurate than 'steffen' where the function is smooth (and monotonic).

The 'cspline' algorithm gives the smoothest interpolations and the most accuracy for very smooth functions, but it tends to overreact to abrupt jumps in the function, leading to unrealistic oscillations around the jump. Monotonic splines are guaranteed to be monotonic between knots and so are unable to overshoot the input data in this way. (The monotonic splines may not be monotonic in the first or last intervals if derivatives are supplied for the endpoints, using the deriv keyword. Also third derivatives (f.D3(x)) are unreliable for monotonic splines.)

The algorithm is irrelevant when there are only two knots. In that case the spline is linear if no derivatives are specified (using deriv), quadratic if one or the other derivative is specified, or cubic if both derivatives are specified.

Examples

Typical usage is:

>>> import gvar as gv
>>> xknot = [0., 0.78539816, 1.57079633, 2.35619449, 3.14159265]
>>> yknot = [0., 0.70710678, 1.0, 0.70710678, 0.]
>>> f = gv.cspline.CSpline(xknot, yknot)
>>> print(f(0.7), f.D(0.7), f.integ(0.7))
0.644243383100598 0.7655924482958195 0.23496394264839268

Here the yknot values were obtained by taking sin(xknot). Tabulating results from the spline together with the exact results shows that this 5-knot spline gives a pretty good approximation of the function sin(x), as well as its derivatives and integral:

x    f(x)    f.D(x)  f.D2(x) f.integ(x) | sin(x)  cos(x)  1-cos(x)
------------------------------------------------------------------
0.3  0.2951  0.9551  -0.2842 0.0446     | 0.2955  0.9553  0.0447
0.5  0.4791  0.8793  -0.4737 0.1222     | 0.4794  0.8776  0.1224
0.7  0.6442  0.7656  -0.6632 0.2350     | 0.6442  0.7648  0.2352
0.9  0.7830  0.6176  -0.7891 0.3782     | 0.7833  0.6216  0.3784
1.1  0.8902  0.4520  -0.8676 0.5461     | 0.8912  0.4536  0.5464
1.3  0.9627  0.2706  -0.9461 0.7319     | 0.9636  0.2675  0.7325
1.5  0.9974  0.0735  -1.0246 0.9286     | 0.9975  0.07074 0.9293

Using the spline outside the range covered by the knots is less good:

>>> print(f(2 * math.pi))
1.7618635470106572

The correct answer is 0.0, of course. Working just outside the knot region is often fine, although it might be a good idea to limit the order of the polynomial used in such regions: for example, setting

>>> f = gv.cspline.CSpline(xknot, yknot, extrap_order=2)

implies that quadratic polynomials are used outside the spline range.

Parameters:

  • xknot (1-d sequence) – Location of the spline’s knots, where the function values are specified. The locations can be numbers of gvar.GVar objects.

  • yknot (1-d sequence) – Function values at the locations specified by xknot[i]. The values can be numbers or gvar.GVar objects.

  • deriv (2-component sequence) – Derivatives at initial and final boundaries of the region specified by xknot[i]. The derivatives can be numbers or gvar.GVar objects. Default value is None for each boundary, which implies natural boundary conditions (vanishing second derivative).

  • extrap_order (int) – Order of polynomial used for extrapolations outside of the spline range. The polynomial is constructed from the spline’s value and derivatives at the (nearest) knot of the spline. The allowed range is 0 <= extrap_order <= 3. The default value is 3 although it is common practice to use smaller values.

  • warn (bool) – If True, warnings are generated when the spline function is called for x values that fall outside of the original range of xknots used to define the spline. Default value is False, which means out-of-range warnings are suppressed.

  • alg (str) – Spline algorithm used, which is one of: 'steffen' (default), 'cspline', and 'pchip'. The first and last of these give monotonic splines.

Linear Algebra

Module gvar.linalg implements several methods for doing basic linear algebra with matrices whose elements can be either numbers or gvar.GVars:

linalg.det(a)

Determinant of matrix a.

Parameters:

a – Two-dimensional, square matrix/array of numbers and/or gvar.GVars.

Returns:

Deterimant of the matrix.

Raises:

ValueError – If matrix is not square and two-dimensional.

linalg.slogdet(a)

Sign and logarithm of determinant of matrix a.

Parameters:

a – Two-dimensional, square matrix/array of numbers and/or gvar.GVars.

Returns:

Tuple (s, logdet) where the determinant of matrix a is

s * exp(logdet).

Raises:

ValueError – If matrix is not square and two-dimensional.

linalg.inv(a)

Inverse of matrix a.

Parameters:

a – Two-dimensional, square matrix/array of numbers and/or gvar.GVars.

Returns:

The inverse of matrix a.

Raises:

ValueError – If matrix is not square and two-dimensional.

linalg.solve(a, b)

Find x such that a @ x = b for matrix a.

Parameters:

  • a – Two-dimensional, square matrix/array of numbers and/or gvar.GVars.

  • b – One-dimensional vector/array of numbers and/or gvar.GVars, or an array of such vectors. Requires b.shape[0] == a.shape[1].

Returns:

The solution x of a.dot(x) = b, which is equivalent to inv(a).dot(b).

Raises:

  • ValueError – If a is not square and two-dimensional.

  • ValueError – If shape of b does not match that of a (that is b.shape[0] != a.shape[1]).

linalg.lstsq(a, b, rcond=None, weighted=False, extrainfo=False)

Least-squares solution x to a @ x = b for gvar.GVars.

Here x is defined to be the solution that minimizes ||b - a @ x||. If b has a covariance matrix, another option is to weight the norm with the inverse covariance matrix: i.e., minimize || isig @ b - isig @ a @ x|| where isig is the square root of the inverse of b’s covariance matrix. Set parameter weighted=True to obtain the weighted-least-squares solution.

Parameters:

  • a – Matrix/array of shape (M,N) containing numbers and/or gvar.GVars.

  • b – Vector/array of shape (M,) containing numbers and/or gvar.GVars.

  • rcond (float) – Cutoff for singular values of a. Singular values smaller than rcond times the maximum eigenvalue are ignored. Default (rcond=None) is max(M,N) times machine precision.

  • weighted (bool) – If True, use weighted least squares; otherwise use unweighted least squares.

  • extrainfo (bool) – If False (default) only x is returned; otherwise (x, residual, rank, s) is returned.

Returns:

Array x of shape (N,) that minimizes || b - a @ x|| if extrainfo==False (default); otherwise returns a tuple (x, residual, rank, s) where residual is the sum of the squares of b - a @ x, rank is the rank of matrix a, and s is an array containing the singular values.

linalg.eigvalsh(a, eigvec=False)

Eigenvalues of Hermitian matrix a.

Parameters:

  • a – Two-dimensional, square Hermitian matrix/array of numbers and/or gvar.GVars. Array elements must be real-valued if gvar.GVars are involved (i.e., symmetric matrix).

  • eigvec (bool) – If True, method returns a tuple of arrays (val, vec) where val[i] are the eigenvalues of a, and vec[:, i] are the mean values of the corresponding eigenvectors. Only val is returned if eigvec=False (default).

Returns:

Array val of eigenvalues of matrix a if parameter eigvec==False (default); otherwise a tuple of arrays (val, vec) where val[i] are the eigenvalues (in ascending order) and vec[:, i] are the mean values of the corresponding eigenvectors.

Raises:

ValueError – If matrix is not square and two-dimensional.

linalg.eigh(a, eigvec=True, rcond=None)

Eigenvalues and eigenvectors of symmetric matrix a.

Parameters:

  • a – Two-dimensional, square Hermitian matrix/array of numbers and/or gvar.GVars. Array elements must be real-valued if gvar.GVars are involved (i.e., symmetric matrix).

  • eigvec (bool) – If True (default), method returns a tuple of arrays (val, vec) where val[i] are the eigenvalues of a (in ascending order), and vec[:, i] are the corresponding eigenvectors of a. Only val is returned if eigvec=False.

  • rcond (float) – Eigenvalues whose difference is smaller than rcond times their sum are assumed to be degenerate (and ignored) when computing variances for the eigvectors. Default (rcond=None) is N times machine precision, where N is the dimension of the array.

Returns:

Tuple (val,vec) of eigenvalues and eigenvectors of matrix a if parameter eigvec==True (default). The eigenvalues val[i] are in ascending order and vec[:, i] are the corresponding eigenvalues. Only the eigenvalues val are returned if eigvec=False.

Raises:

ValueError – If matrix is not square and two-dimensional.

linalg.svd(a, b, compute_uv=True, rcond=None)

svd decomposition of matrix a containing gvar.GVars.

Parameters:

  • a – Two-dimensional matrix/array of numbers and/or gvar.GVars.

  • compute_uv (bool) – It True (default), returns tuple (u,s,vT) where matrix a = u @ np.diag(s) @ vT where matrices u and vT satisfy u.T @ u = 1 and vT @ vT.T = 1, and s is the list of singular values. Only s is returned if compute_uv=False.

  • rcond (float) – Singular values whose difference is smaller than rcond times their sum are assumed to be degenerate for calculating variances for u and vT. Default (rcond=None) is max(M,N) times machine precision.

Returns:

Tuple (u,s,vT) where matrix a = u @ np.diag(s) @ vT where matrices u and vT satisfy u.T @ u = 1 and vT @ vT.T = 1, and s is the list of singular values. If a.shape=(N,M), then u.shape=(N,K) and vT.shape=(K,M) where K is the number of nonzero singular values (len(s)==K). If compute_uv==False only s is returned.

Raises:

ValueError – If matrix is not two-dimensional.

Ordinary Differential Equations

Module gvar.ode implements two classes for integrating systems of first-order differential equations using an adaptive Runge-Kutta algorithm. One integrates scalar- or array-valued equations, while the other integrates dictionary-valued equations:

class gvar.ode.Integrator(deriv, tol=1e-05, h=None, hmin=None, hmax=None, maxstep=None, delta=None, analyzer=None)

Integrate dy/dx = deriv(x,y).

An Integrator object odeint integrates dy/dx = f(x,y) to obtain y(x1) from y0 = y(x0). y and f(x,y) can be scalars or numpy arrays. Typical usage is illustrated by the following code for integrating dy/dx = y:

from gvar.ode import Integrator

def f(x, y):
    return y

odeint = Integrator(deriv=f,  tol=1e-8)
y0 = 1.
y1 = odeint(y0, interval=(0, 1.))
y2 = odeint(y1, interval=(1., 2.))
...

Here the first call to odeint integrates the differential equation from x=0 to x=1 starting with y=y0 at x=0; the result is y1=exp(1), of course. Similarly the second call to odeint continues the integration from x=1 to x=2, giving y2=exp(2).

If the interval is a list with more than two entries, then odeint(y0, interval=[x0, x1, x2 ...]) in the example above returns an array of solutions for points x1, x2 .... So the example above could have been written equivalently as

...

odeint = Integrator(deriv=f,  tol=1e-8)
y0 = 1.
y1, y2 ... = odeint(y0, interval=[0, 1., 2. ...])

An alternative interface creates a new function which is the solution of the differential equation for specific initial conditions. The code above could be rewritten:

x0 = 0.         # initial conditions
y0 = 1.
y = Integrator(deriv=f, tol=1e-8).solution(x0, y0)
y1 = y(1)
y2 = y(2)
...

Here method Integrator.solution() returns a function y(x) where: a) y(x0) = y0; and b) y(x) uses the integator to integrate the differential equation to point x starting from the last point at which y was evaluated (or from x0 for the first call to y(x)). The function can also be called with an array of x values, in which case an array containing the corresponding y values is returned.

The integrator uses an adaptive Runge-Kutta algorithm that adjusts the integrator’s step size to obtain relative accuracy tol in the solution. An initial step size can be set in the Integrator by specifying parameter h. A minimum step size hmin can also be specified; the Integrator raises an exception if the step size becomes smaller than hmin. The Integrator keeps track of the number of good steps, where h is increased, and the number of bad steps, where h is decreased and the step is repeated: odeint.ngood and odeint.nbad, respectively.

A custom criterion for step-size changes can be implemented by specifying a function for parameter delta. This is a function delta(yerr, y, delta_y) — of the estimated error yerr after a given step, the proposed value for y, and the proposed change delta_y in y — that returns a number to compare with tolerance tol. The step size is decreased and the step repeated if delta(yerr, y, delta_y) > tol; otherwise the step is accepted and the step size increased. The default definition of delta is equivalent to:

import numpy as np
import gvar as gv

def delta(yerr, y, delta_y):
    return np.max(
        np.fabs(yerr) /
        (np.fabs(y) + np.fabs(delta_y) + gv.ode.TINY)
        )

An analyzer analyzer(x,y) can be specified using parameter analyzer. This function is called after every full step of the integration, with the current values of x and y. Objects of type gvar.ode.Solution are examples of (simple) analyzers.

Parameters:

  • deriv – Function of x and y that returns dy/dx. The return value should have the same shape as y if arrays are used.

  • tol (float) – Relative accuracy in y relative to |y| + h|dy/dx| for each step in the integration. Any integration step that achieves less precision is repeated with a smaller step size. The step size is increased if precision is higher than needed. Default is 1e-5.

  • h (float or None) – Absolute value of initial step size. The default value equals the entire width of the integration interval.

  • hmin (float or None) – Smallest step size allowed. A warning is raised if a smaller step size is requested, and the step size is not decreased. This prevents infinite loops at singular points, but the solution may not be reliable when a warning has been issued. The default value is None (which does not prevent infinite loops).

  • hmax (float or None) – Largest step allowed. Ignored if set to None.

  • maxstep (int or None) – Maximum number of integration steps allowed, after which a RuntimeError exception is raised. Ignored if set to None.

  • delta – Function delta(yerr, y, delta_y) that returns a number to be compared with tol at each integration step: if it is larger than tol, the step is repeated with a smaller step size; if it is smaller the step is accepted and a larger step size used for the subsequent step. Here yerr is an estimate of the error in y on the last step; y is the proposed value; and delta_y is the change in y over the last step.

  • analyzer – Function of x and y that is called after each step of the integration. This can be used to analyze intermediate results.

class gvar.ode.DictIntegrator(deriv, tol=1e-05, h=None, hmin=None, hmax=None, maxstep=None, delta=None, analyzer=None)

Integrate dy/dx = deriv(x,y) where y is a dictionary.

An DictIntegrator object odeint integrates dy/dx = f(x,y) to obtain y(x1) from y0 = y(x0). y and f(x,y) are dictionary types having the same keys, and containing scalars and/or numpy arrays as values. Typical usage is:

from gvar.ode import DictIntegrator

def f(x, y):
    ...

odeint = DictIntegrator(deriv=f,  tol=1e-8)
y1 = odeint(y0, interval=(x0, x1))
y2 = odeint(y1, interval=(x1, x2))
...

The first call to odeint integrates from x=x0 to x=x1, returning y1=y(x1). The second call continues the integration to x=x2, returning y2=y(x2). Multiple integration points can be specified in interval, in which case a list of the corresponding y values is returned: for example,

odeint = DictIntegrator(deriv=f,  tol=1e-8)
y1, y2 ... = odeint(y0, interval=[x0, x1, x2 ...])

The integrator uses an adaptive Runge-Kutta algorithm that adjusts the integrator’s step size to obtain relative accuracy tol in the solution. An initial step size can be set in the DictIntegrator by specifying parameter h. A minimum ste psize hmin can also be specified; the Integrator raises an exception if the step size becomes smaller than hmin. The DictIntegrator keeps track of the number of good steps, where h is increased, and the number of bad steps, where h is decreases and the step is repeated: odeint.ngood and odeint.nbad, respectively.

An analyzer analyzer(x,y) can be specified using parameter analyzer. This function is called after every full step of the integration with the current values of x and y. Objects of type gvar.ode.Solution are examples of (simple) analyzers.

Parameters:

  • deriv – Function of x and y that returns dy/dx. The return value should be a dictionary with the same keys as y, and values that have the same shape as the corresponding values in y.

  • tol (float) – Relative accuracy in y relative to |y| + h|dy/dx| for each step in the integration. Any integration step that achieves less precision is repeated with a smaller step size. The step size is increased if precision is higher than needed. Default is 1e-5.

  • h (float or None) – Absolute value of initial step size. The default value equals the entire width of the integration interval.

  • hmin (float or None) – Smallest step size allowed. A warning is raised if a smaller step size is requested, and the step size is not decreased. This prevents infinite loops at singular points, but the solution may not be reliable when a warning has been issued. The default value is None (which does not prevent infinite loops).

  • hmax (float or None) – Largest step allowed. Ignored if set to None.

  • maxstep (int or None) – Maximum number of integration steps allowed, after which a RuntimeError exception is raised. Ignored if set to None.

  • delta – Function delta(yerr, y, delta_y) that returns a number to be compared with tol at each integration step: if it is larger than tol, the step is repeated with a smaller step size; if it is smaller the step is accepted and a larger step size used for the subsequent step. Here yerr is an estimate of the error in y on the last step; y is the proposed value; and delta_y is the change in y over the last step.

  • analyzer – Function of x and y that is called after each step of the integration. This can be used to analyze intermediate results.

A simple analyzer class is:

class gvar.ode.Solution

ODE analyzer for storing intermediate values.

Usage: eg, given

odeint = Integrator(...)
soln = Solution()
y0 = ...
y = odeint(y0, interval=(x0, x), analyzer=soln)

then the soln.x[i] are the points at which the integrator evaluated the solution, and soln.y[i] is the solution of the differential equation at that point.

One-Dimensional Integration

Module gvar.ode also provides a method for evaluating one-dimensional integrals (using its adaptive Runge-Kutta algorithm):

gvar.ode.integral(fcn, interval, fcnshape=None, tol=1e-08, hmin=None)

Compute integral of fcn(x) on interval.

Given a function fcn(x) the call

result = integral(fcn, interval=(x0, x1))

calculates the integral of fcn(x) from x0 to x1. For example:

>>> def fcn(x):
...    return math.sin(x) ** 2 / math.pi
>>> result = integral(fcn, (0, math.pi))
>>> print(result)
0.500000002834

Function fcn(x) can return a scalar or an array (any shape): for example,

>>> def fcn(x):
...    return np.array([1., x, x**3])

>>> result = integral(fcn, (0,1))
>>> print(result)
[1. 0.5 0.25]

The function can also return dictionaries whose values are scalars or arrays: for example,

>>> def fcn(x):
...    return dict(x=x, x3=x**3)
>>> result = integral(fcn, (0,1))
>>> print(result)
{'x': 0.5,'x3': 0.25}
Parameters:

  • fcn – Function of scalar variable x that returns the integrand. The return value should be either a scalar or an array, or a dictionary whose values are scalars and/or arrays.

  • interval – Contains the interval (x0,x1) over which the integral is computed.

  • fcnshape – Contains the shape of the array returned by f(x) or () if the function returns a scalar. Setting fshape=None (the default) results in an extra function evaluation to determine the shape.

  • tol – Relative accuracy of result.

  • hmin – Smallest step size allowed in adaptive integral. A warning is raised if a smaller step size is requested, and the step size is not decreased. This prevents infinite loops at singular points, but the integral may not be accurate when a warning has been issued. The default value is None (which does not prevent infinite loops).

Pade Approximants

Module gvar.pade provides a class to represent Pade approximants of functions:

class gvar.pade.Pade(f, order, rtol='gavg')

Pade approximant to sum_i f[i] x**i for GVars.

The order=(m,n) Pade approximant to a series given by sum_i f[i] * x**i is the ratio of polynomials of order m (numerator) and n (denominator) whose Taylor expansion agrees with that of the original series up to order m+n.

A Pade object pade creates gvar.powerseries.PowerSeries objects for the numerator (pade.num) and the denominator (pade.den) of the Pade approximant corresponding to the input series f[i]. The approximant can be evaluated for arbitrary x using pade(x). The coefficients used in the numerator and denominator are given by pade.num.c and pade.den.c, respectively.

Elements in the series f[i] may be numbers or gvar.GVars. When the latter appear, the code uses an SVD algorithm (see pade_svd()) to deal with the imprecision in the input data. It automatically reduces the order of the approximant if the extraction of Pade coefficients is too unstable given the noise in the input data. The actual order used in the approximant is given by pade.order.

Examples

Typical usage is:

>>> import gvar as gv
>>> c = gv.gvar(['1(0)', '1.0(1)', '0.500(10)', '0.1667(100)','.04167(100)'])
>>> pade = gv.pade.Pade(c, order=(2,2))
>>> print(pade.num.c, pade.den.c)               # num and den coefficients
[1(0) 0.50(14) 0.08(11)] [1(0) -0.50(14) 0.083(55)]
>>> print(pade(1))                              # evaluate pade at x=1
2.714(99)
>>> print(pade.num(1), pade.den(1), pade.order) # evaluate num and den at x=1
1.58(21) 0.583(82) (2, 2)

When errors are larger (see c[0]), the order may be automatically reduced (here to (1,2)):

>>> c = gv.gvar(['1.0(1)', '1.0(1)', '0.500(10)', '0.1667(100)','.04167(100)'])
>>> pade = gv.pade.Pade(c, order=(2,2))
>>> print(pade.num.c, pade.den.c)               # num and den coefficients
[1.00(10) 0.33(29)] [1(0) -0.67(17) 0.17(11)]
>>> print(pade(1))                              # evaluate pade at x=1
2.67(20)
>>> print(pade.num(1), pade.den(1), pade.order) # evaluate num and den at x=1
1.33(27) 0.500(69) (1, 2)
Parameters:

  • f – Array f[i] of power series coefficients for i=0...n+m.

  • order (tuple) – order=(m,n) specifies that the numerator and denominator of the Pade approximant have maximum order m and n, respectively.

  • rtol (float or str) – Relative uncertainty in the coefficients, for use by the SVD algorithm. If rtol is a string, it determines how the relative tolerance is determined from the relative uncertainties in the f[i]. Set rtol equal to: 'gavg' for the geometric mean (default); 'avg' for the average; 'min' for the minimum; or 'max' for the maximum. Otherwise a number can be specified, in which case the uncertainties in f[i] are ignored. The SVD analysis is skipped if rtol=None.

gvar.pade.pade_gvar(f, m, n, rtol='gavg')

(m,n) Pade approximant to sum_i f[i] x**i for GVars.

The (m,n) Pade approximant to a series given by sum_i f[i] * x**i is the ratio of polynomials of order m (numerator) and n (denominator) whose Taylor expansion agrees with that of the original series up to order m+n.

This code uses an SVD algorithm (see pade_svd()) to deal with imprecision in the input data. It automatically reduces the order of the approximant if the extraction of Pade coefficients is too unstable given noise in the input data.

Parameters:

  • f – Array f[i] of power series coefficients for i=0...n+m.

  • m – Maximum order of polynomial in numerator of Pade approximant (m>=0).

  • n – Maximum order of polynomial in denominator of Pade approximant (m>=0).

  • rtol (float or str) – If rtol is a string, it determines how the relative tolerance is determined from the relative uncertainties in the f[i]. Set rtol equal to: 'gavg' for the geometric mean (default); 'avg' for the average; 'min' for the minimum; or 'max' for the maximum. Otherwise a number can be specified, in which case the uncertainties in f[i] are ignored. The SVD analysis is skipped if rtol=None.

Returns:

Tuple of power series coefficients (p, q) such that sum_i p[i] x**i is the numerator of the approximant, and sum_i q[i] x**i is the denominator. q[0] is normalized to 1.

gvar.pade.pade_svd(f, m, n, rtol=1e-14)

(m,n) Pade approximant to sum_i f[i] x**i.

The (m,n) Pade approximant to a series given by sum_i f[i] * x**i is the ratio of polynomials of order m (numerator) and n (denominator) whose Taylor expansion agrees with that of the original series up to order m+n.

This code is adapted from P. Gonnet, S. Guttel, L. N. Trefethen, SIAM Review Vol 55, No. 1, 101 (2013). It uses an SVD algorithm to deal with imprecision in the input data, here specified by the relative tolerance rtol for the input coefficients f[i]. It automatically reduces the order of the approximant if the extraction of Pade coefficients is too unstable given tolerance rtol.

Parameters:

  • f – Array f[i] of power series coefficients for i=0...n+m.

  • m – Maximum order of polynomial in numerator of Pade approximant (m>=0).

  • n – Maximum order of polynomial in denominator of Pade approximant (m>=0).

  • rtol – Relative accuracy of input coefficients. (Default is 1e-14.)

Returns:

Tuple of power series coefficients (p, q) such that sum_i p[i] x**i is the numerator of the approximant, and sum_i q[i] x**i is the denominator. q[0] is normalized to 1.

Power Series

Module gvar.powerseries provides tools for manipulating power series approximations of functions. A function’s power series is specified by the coefficients in its Taylor expansion with respect to an independent variable, say x:

f(x) = f(0) + f'(0) * x + (f''(0)/2) * x**2 + (f'''(0)/6) * x**3 + ...
     = f0 + f1 * x + f2 * x**2 + f3 * x**3 + ...

In practice a power series is different from a polynomial because power series, while infinite order in principle, are truncated at some finite order in numerical applications. The order of a power series is the highest power of x that is retained in the approximation; coefficients for still higher-order terms are assumed to be unknown (as opposed to zero).

Taylor’s theorem can be used to generate power series for functions of power series:

g(f(x)) = g(f0) + g'(f0) * (f(x)-f0) + (g''(f0)/2) * (f(x)-f0)**2 + ...
        = g0 + g1 * x + g2 * x**2 + ...

This allows us to define a full calculus for power series, where arithmetic expressions and (sufficiently differentiable) functions of power series return new power series.

Using power series

Class PowerSeries provides a numerical implementation of the power series calculus. PowerSeries([f0,f1,f2,f3...]) is a numerical representation of a power series with coefficients f0, f1, f2, f3... (as in f(x) above). Thus, for example, we can define a 4th-order power series approximation f to exp(x)=1+x+x**2/2+... using

>>> import gvar as gv
>>> import gvar.powerseries as ps
>>> f = ps.PowerSeries([1., 1., 1/2., 1/6., 1/24.])
>>> print(f)            # print the coefficients
[ 1.          1.          0.5         0.16666667  0.04166667]

Arithmetic expressions involving instances of class PowerSeries are themselves PowerSeries as in, for example,

>>> print(1/f)              # power series for exp(-x)
[ 1.         -1.          0.5        -0.16666667  0.04166667]
>>> print(gv.log(f))        # power series for x
[ 0.  1.  0. -0.  0.]
>>> print(f / f)            # power series for 1
[ 1.  0.  0.  0.  0.]

The standard arithmetic operators (+,-,*,/,=,**) are supported, as are the usual elementary functions (exp, log, sin, cos, tan ...). Different PowerSeries can be combined arithmetically to create new PowerSeries; the order of the result is that of the operand with the lowest order.

PowerSeries can be differentiated and integrated:

>>> print(f.deriv())    # derivative of exp(x)
[ 1.          1.          0.5         0.16666667]
>>> print(f.integ())    # integral of exp(x) (from x=0)
[ 0.          1.          0.5         0.16666667  0.04166667  0.00833333]

Each PowerSeries represents a function. The PowerSeries for a function of a function is easily obtained. For example, assume f represents function f(x)=exp(x), as above, and g represents g(x)=log(1+x):

>>> g = ps.PowerSeries([0, 1, -1/2., 1/3., -1/4.])

Then f(g) gives the PowerSeries for exp(log(1+x)) = 1 + x:

>>> print(f(g))
[  1.0000e+00   1.0000e+00   0.0000e+00  -2.7755e-17 -7.6327e-17]

Individual coefficients from the powerseries can be accessed using array-element notation: for example,

>>> print(f[0], f[1], f[2], f[3])
1.0 1.0 0.5 0.166666666667
>>> f[0] = f[0] - 1.
>>> print(f)            # f is now the power series for exp(x)-1
[ 0.          1.          0.5         0.16666667  0.04166667]

Finally, a power series can be evaluated for a particular numerical value of x:

>>> x = 0.01
>>> print(f(x))             # should be exp(0.01)-1 approximately
0.0100501670833
>>> print(gv.exp(x)-1)      # verify that it is
0.0100501670842

The independent variable x could be of any arithmetic type (it need not be a float).

Multivariate power series

The coefficients in a PowerSeries object can themselves by PowerSeries objects. The is used to represent multivariate power series such as:

f(x,y) = f00 + f10 * x + f01 * y + f20 * x**2 + f11 * x*y + f02 * y**2 + ...

One way to construct a PowerSeries object f representing this series, through order=2, is from an array c containing the coefficients:

c = [[f00, f01, f02], [f10, f11, 0], [f20, 0, 0]]
f = ps.multiseries(c, order=2)

Here entries for c[1,2], c[2,1], and c[2,2] are ignored because they correspond to order=3 or higher. The individual coefficients c[i,j] are accessed using f[i,j], and the function is evaluate at point (x,y) using f(x,y). Similarly the first-order partial derivative with respect to x and y, for example, is given by PowerSeries object f.deriv(1,1), while first-order integrals with respect to x and y are given by f.integ(1,1).

Taylor expansions of Python functions

PowerSeries can be used to compute Taylor series for more-or-less arbitrary pure-Python functions provided the functions are sufficiently differentiable. To compute the N-th order expansion of a Python function g(x), first create a N-th order PowerSeries variable that represents the expansion parameter: say, x = PowerSeries([0.,1.],order=N). The Taylor series for function g is then given by g_taylor = g(x) which is a PowerSeries instance. For example, consider:

>>> def g(x):              # an example of a Python function
...     return 0.5/sqrt(1+x) + 0.5/sqrt(1-x)
...
>>> x = ps.PowerSeries([0.,1.],order=5)    # Taylor series for x
>>> print(x)
[ 0.  1.  0.  0.  0.  0.]
>>> g_taylor = g(x)        # Taylor series for g(x) about x=0
>>> print(g_taylor)
[ 1.         0.         0.375      0.         0.2734375  0.       ]
>>> exp_taylor = gv.exp(x) # Taylor series for exp(x) about x=0
>>> print(exp_taylor)
[ 1.          1.          0.5         0.16666667  0.04166667  0.00833333]

This generalizes easily to multivariate expansions. For example, one can calculate the Taylor expansion coefficients for exp(x+y) using:

>>> x,y = ps.multivar(dim=2, order=3)
>>> exp_taylor = gv.exp(x + y)
>>> print(exp_taylor)
[[1.         1.         0.5        0.16666667], [1.  1.  0.5], [0.5 0.5], [0.16666667]]
>>> print(exp_taylor(.1,.2), gv.exp(.1 + .2))
1.3495000000000001 1.3498588075760032

Here function multivar() creates PowerSeries objects corresponding to the expansion variables through a given order.

class gvar.powerseries.PowerSeries(c=None, order=None)

Power series representation of a function.

The power series created by PowerSeries(c) corresponds to:

c[0] + c[1]*x + c[2]*x**2 + ... .

The order of the power series is normally determined by the length of the input list c. This can be overridden by specifying the order of the power series using the order parameter. The list of c[i]s is then padded with zeros if c is too short, or truncated if it is too long. Omitting c altogether results in a power series all of whose coefficients are zero. Individual series coefficients are accessed using array/list notation: for example, the 3rd-order coefficient of PowerSeries p is p[3]. The order of p is p.order. PowerSeries should work for coefficients of any data type that supports ordinary arithmetic.

Arithmetic expressions of PowerSeries variables yield new PowerSeries results that represent the power series expansion of the expression. Expressions can include the standard mathematical functions (log, exp, sqrt, sin, cos, tan...). PowerSeries can also be differentiated (p.deriv()) and integrated (p.integ()).

Parameters:

  • c (array) – Power series coefficients.

  • order (int or None) – Highest power in power series. If None, the order is inferred from the array of coefficients. If array c is too small for the specified order, the array is padded with zeros at the end.

coeff

Copy of power series coefficients (numpy.array).

deriv(*n)

Compute n-th derivative (or partial derivative) of self.

Parameters:

n (array) – Number of derivatives in each direction. Default is n=[1].

Returns:

PowerSeries object representing the n-th derivative or partial derivative of self.

integ(*n, x0=0.)

Compute n-th indefinite integral of self.

If x0 is specified, then the definite integral, integrating from point x0, is returned.

Parameters:

  • n (array) – Number of integrations in each direction. Default is n=[1].

  • x0 (array or float) – Starting point for definite integral in each direction (default is 0).

Returns:

PowerSeries object representing the n-th integral of self.

gvar.powerseries.multiseries(c, order=None)

Create multivariate power series from coefficients in array c.

Parameters:

  • c (array) – numpy-like array containing the power series coefficients. In d dimensions, c[i1,i2,...,id] is the coefficient of x1**i1 * x2**i2 * ... * xd**id.

  • order (int or None) – Highest power in power series, where the power associated with term x1**i1 * x2**i2 * ... * xd**id is the sum of the exponents: i1 + i2 + ... + id. If None, the order is inferred from the array of coefficients.

Returns:

PowerSeries object representing the multivariate power series.

gvar.powerseries.multivar(dim, order)

Create PowerSeries objects representing the expansion variables.

Parameters:

  • dim (int) – The dimensionality of the multivariate space.

  • order (int) – Highest power in the power series, where the power associated with term x1**i1 * x2**i2 * ... * xd**id is the sum of the exponents: i1 + i2 + ... + id.

Returns:

An array of dim PowerSeries objects corresponding to each of the expansion variables in a dim-dimensional multivariate power series.

Root Finding

Module gvar.root contains methods for finding the roots of of one-dimensional functions: that is, finding x such that fcn(x)=0 for a given function fcn. Typical usage is:

>>> import math
>>> import gvar as gv
>>> interval = gv.root.search(math.sin, 1.)     # bracket root
>>> print(interval)
(3.1384283767210035, 3.4522712143931042)
>>> root = gv.root.refine(math.sin, interval)   # refine root
>>> print(root)
3.14159265359

This code finds the first root of sin(x)=0 larger than 1. The first setp is a search to find an interval containing a root. Here gvar.root.search() examines sin(x) for a sequence of points 1. * 1.1 ** n for n=0,1,2..., stopping when the function changes sign. The last two points in the sequence then bracket a root since sin(x) is continuous; they are returned as a tuple to interval. The final root is found by refining the interval, using gvar.root.refine. By default, the root is refined iteratively to machine precision, but this requires only a small number (4) of iterations:

>>> print(root.nit)                             # number of iterations
4

The most challenging situations are ones where the function is extremely flat in the vicinity of the root — that is, two or more of its leading derivatives vanish there. For example:

>>> import gvar as gv
>>> def f(x):
...     return (x + 1) ** 3 * (x - 0.5) ** 11
>>> root = gv.root.refine(f, (0, 2))
>>> print(root)
0.5
>>> print(root.nit)                             # number of iterations
142

This routine also works with variables of type gvar.GVar: for example,

>>> import gvar as gv
>>> def f(x, w=gv.gvar(1, 0.1)):
...     return gv.sin(w * x)
>>> root = gv.root.refine(f, (1, 4))
>>> print(root)
3.14(31)

returns a root with a 10% uncertainty, reflecting the uncertainty in parameter w.

Descriptions of the two methods follow.

root.search(x0, incr=0, fac=1.1, maxit=100, analyzer=None)

Search for and bracket root of one-dimensional function fcn(x).

This method searches for an interval in x that brackets a root of fcn(x)=0. It examines points

x[j + 1] = fac * x[j] + incr

where x[0]=x0 and j=0...maxit-1, looking for a pair of successive points where fcn(x[j]) changes sign. These points bracket a root (assuming the function is continuous), providing a coarse estimate of the root. That estimate can be refined using root.refine().

Example

The following code seeks to bracket the first zero of sin(x) with x>0.1:

>>> import math
>>> import gvar as gv
>>> interval = gv.root.search(math.sin, 0.1)
>>> print(interval)
(3.0912680532870755, 3.4003948586157833)

The resulting interval correctly brackets the root at pi.

Parameters:

  • fcn – One dimenionsal function whose root is sought.

  • x0 (float) – Starting point for search.

  • incr (float, optional) – Increment used for linear searches. Default value is 0.

  • fac (float, optional) – Rescaling factor for exponential searches. Default value is 1.1.

  • maxit (int, optional) – Maximum number of steps allowed for search. An exception is raised if a root is not found in time. Default value is 100.

  • analyzer – Optional function f(x, fcn(x)) that is called for each point x that is examined. This can be used, for example, to monitor the search while debugging. Default is None.

Returns:

Tuple (a, b) where fcn(a) * fcn(b) <= 0, which implies that a root occurs between a and b (provided the function is continuous). The tuple has extra attributes that provide additional information about the search:

  • nit — Number of iterations used to find interval (a,b).

  • fcnval — Tuple containing the function values at (a,b).

Raises:

RuntimeError – If unable to find a root in maxit steps.

root.refine(interval, rtol=None, maxit=1000, analyzer=None)

Find root x of one-dimensional function fcn on an interval.

This method finds a root x of fcn(x)=0 inside an interval=(a,b) that brackets the root, with fcn(a) * fcn(b) <= 0.

This method is a pure Python adaptation of an algorithm from Richard Brent’s book “Algorithms for Minimization without Derivatives” (1973). Being pure Python it works with gvar.GVar-valued functions and variables.

Example

The following code finds a root of sin(x) in the interval 1 <= x <= 4, using 7 iterative refinements of the initial interval:

>>> import math
>>> import gvar as gv
>>> root = gv.root.refine(math.sin, (1, 4))
>>> print(root)
3.14159265359
>>> print(root.nit)
7
Parameters:

  • fcn – One-dimensional function whose zero/root is sought.

  • interval – Tuple (a,b) specifying an interval containing the root, with fcn(a) * fcn(b) <= 0. The search for a root is confined to this interval.

  • rtol (float, optional) – Relative tolerance for the root. The default value is None, which sets rtol equal to machine precision (sys.float_info.epsilon). A larger value usually leads to less precision but is faster.

  • maxit (int, optional) – Maximum number of iterations used to find a root with the given tolerance. A warning is issued if the algorithm does not converge in time. (Default value is 1000.)

  • analyzer – Optional function f(x, fcn(x)) that is called for each point x examined by the algorithm. This can be used, for example, to monitor convergence while debugging. Default is None.

Returns:

The root, which is either a float or a gvar.GVar but with extra attributes that provide additional information about the root:

  • nit — Number of iterations used to find the root.

  • interval — Smallest interval (b,c) found containing the root, where b is the root returned by the method.

  • fcnval — Value of fcn(x) at the root.

Raises:

  • ValueError – If fcn(a) * fcn(b) > 0 for initial interval (a,b).

  • UserWarning – If the algorithm fails to converge after maxit iterations.