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(xknots, yknots, deriv=(None, None), extrap_order=3, warn=True)

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) is continuous, as are its first and second derivatives. Function f(x) is a cubic polynomial between the knots xknot[i].

CSpline(xknot, yknot) creates a natural spline, which has zero second derivative at the end points, xknot[0] and xknot[-1] (assuming the knots are sorted). More generally one can specify the derivatives of f(x) at one or both of the endpoints:

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]. Replacing either (or both) of these with None results in a derivative corresponding to zero second derivative at that boundary (i.e., a natural boundary).

Derivatives and integrals of the spline function can also be evaluated:

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

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

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

Splines can be used outside the range covered by the defining xknot values. As this is often a bad idea, the CSpline methods issue a warning when called with out-of-range points. The warning can be suppressed by setting parameter warn=False. 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 can be specified by setting parameter extrap_order to a (non-negative) integer less than 3; this is often a good idea.

Examples

Typical usage is:

>>> import math
>>> 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.D2(0.7), f.integ(0.7))
0.644243383101 0.765592448296 -0.663236750777 0.234963942648

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.04458    | 0.2955  0.9553  0.04466
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.235      | 0.6442  0.7648  0.2352
0.9  0.783   0.6176  -0.7891 0.3782     | 0.7833  0.6216  0.3784
1.1  0.8902  0.452   -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.07352 -1.025  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))
gvar/cspline.py:164: UserWarning: x outside of spline range: [ 6.28318531]
1.7618635470106501

The correct answer is 0.0, of course. This is why the spline function issues a warning. Working just outside the knot region is often fine, although it is usually 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. Finally one can specify the values of the first derivatives of the function at one or the other endpoints of the spline region, if they are known. Continuing from above, for example, one would take

>>> f = gv.cspline.CSpline(xknot, yknot, deriv=[1., -1.])

since the derivatives of sin(x) at x=0 and x=3.14159265 are 1 and -1, respectively.

Parameters:
  • xknot (1-d sequence of number) – The knots of the spline, where the function values are specified. The knots are sorted (from small to large) if necessary.
  • yknot (1-d sequence of number) – Function values at the locations specified by xknot[i].
  • deriv (2-component sequence) – Derivatives at initial and final boundaries of the region specified by xknot[i]. Default value is None for each boundary.
  • 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 True; out-of-range warnings are suppressed if set to False.

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.dot(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.eigvalsh(a, eigvec=False)

Eigenvalues of Hermitian matrix a.

Parameters:
  • a – Two-dimensional, square matrix/array of numbers and/or gvar.GVars.
  • eigvec (bool) – If True, method returns a tuple of arrays (val, vec) where the val[i] are the eigenvalues. Arrays vec[:, i] are the corresponding eigenvectors of a when one ignores uncertainties (that is, they are eigenvectors of gvar.mean(a)). Only val is returned if eigvec=False (default).
Returns:

Array of eigenvalues of matrix a if parameter eigvec==False (default). where the val[i] are the eigenvalues; otherwise it returns a tuple of arrays (val, vec) where the val[i] are the eigenvalues. Arrays vec[:, i] are the corresponding eigenvectors of a when one ignores uncertainties (that is, they are eigenvectors of gvar.mean(a)).
Raises:

ValueError – If matrix is not square and 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, 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 roughly 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)
        )

A custom definition can be used to allow an Integrator to work with data types other than floats or numpy arrays of floats. All that is required of the data type is that it support ordinary arithmetic. Therefore, for example, defining delta(yerr, y, delta_y) with np.abs() instead of np.fabs() allows y to be complex valued. (Actually the default delta allows this as well.)

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).
  • 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, 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.
  • h (float) – Absolute value of initial step size. The default value equals the entire width of the integration interval.
  • hmin (float) – 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).
  • 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):

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

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.

Power series arithmetic

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

>>> from gvar.powerseries import *
>>> f = 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 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 = 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]

Numerical evaluation of power series

The power series can also be evaluated for a particular numerical value of x: continuing the example,

>>> x = 0.01
>>> print f(x)          # should be exp(0.01)-1 approximately
0.0100501670833

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

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 locally analytic (or at least 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:

>>> from gvar.powerseries import *
>>> def g(x):           # an example of a Python function
...     return 0.5/sqrt(1+x) + 0.5/sqrt(1-x)
...
>>> x = 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 = exp(x) # Taylor series for exp(x) about x=0
>>> print exp_taylor
[ 1.          1.          0.5         0.16666667  0.04166667  0.00833333]
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 (list or array) – Power series coefficients (optional if parameter order specified).
  • order (integer) – Highest power in power series (optional if parameter c specified).
deriv(n=1)

Compute n-th derivative of self.

Parameters:n (positive integer) – Number of derivatives.
Returns:n-th derivative of self.
integ(n=1, x0=None)

Compute n-th indefinite integral of self.

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

Parameters:
  • n (integer) – Number of integrations.
  • x0 – Starting point for definite integral (optional).
Returns:

n-th integral of self.
order

Highest power in 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(fcn, 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(fcn, 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.