Numerical Analysis Modules in gvar
¶
gvar.GVar
s 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.GVar
s, perhaps with minor modifications.
Here we describe some sample numerical codes, included in
gvar
, that have been adapted
to work with gvar.GVar
s, as well as with float
s.
More examples will follow with time.
The sub-modules included here are:
gvar.cspline
— cubic splines for 1-d data.
gvar.linalg
— basic linear algebra.
gvar.ode
— integration of systems of ordinary differential equations; one-dimensional integrals.
gvar.powerseries
— power series representation of functions.
gvar.root
— root-finding for one-dimensional functions.
See Case Study: Creating an Integrator for an example
illustrating how existing code can be repurposed to work with gvar.GVar
s.
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 functionyknot[i]
atN
pointsxknot[i]
fori=0..N-1
(the knots of the spline), the codefrom gvar.cspline import CSpline f = CSpline(xknot, yknot)
defines a function
f
such that: a)f(xknot[i]) = yknot[i]
for alli
; and b)f(x)
and its first derivative are continuous. Functionf(x)
is a cubic polynomial between the knotsxknot[i]
(determined from the values and first derivatives at the adjacent knots). Argumentx
inf(x)
is a number or an array of numbers (any shape). Numbers can be replaced bygvar.GVar
objects inx
,xknot
and/oryknot
.Derivatives and integrals of the spline function are also available:
f.D(x)
— first derivative atx
;f.D2(x)
— second derivative atx
;f.D3(x)
— third derivative atx
;f.integ(x)
— integral fromxknot[0]
tox
.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 parameterwarn
can be set toTrue
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 aCSpline
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 parameterextrap_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]
andxknot[-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 atxknot[0]
anddydx_f
is the derivative atxknot[-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 tof'(x)
andf(x)
). Setsf''(x)
to zero at the endpoints (natural boundary conditions) if the first derivatives have not been specified there (using keywordderiv
).alg='pchip'
Monotonic cubic spline used in the
scipy
functionscipy.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 thederiv
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 takingsin(xknot)
. Tabulating results from the spline together with the exact results shows that this 5-knot spline gives a pretty good approximation of the functionsin(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 orgvar.GVar
objects.deriv (2-component sequence) – Derivatives at initial and final boundaries of the region specified by
xknot[i]
. The derivatives can be numbers orgvar.GVar
objects. Default value isNone
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 forx
values that fall outside of the original range ofxknot
s used to define the spline. Default value isFalse
, 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.GVar
s:
- linalg.det(a)¶
Determinant of matrix
a
.- Parameters:
a – Two-dimensional, square matrix/array of numbers and/or
gvar.GVar
s.- 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.GVar
s.- Returns:
- Tuple
(s, logdet)
where the determinant of matrixa
is s * exp(logdet)
.
- Tuple
- 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.GVar
s.- Returns:
The inverse of matrix
a
.- Raises:
ValueError – If matrix is not square and two-dimensional.
- linalg.solve(a, b)¶
Find
x
such thata @ x = b
for matrixa
.- Parameters:
- Returns:
The solution
x
ofa.dot(x) = b
, which is equivalent toinv(a).dot(b)
.- Raises:
ValueError – If
a
is not square and two-dimensional.ValueError – If shape of
b
does not match that ofa
(that isb.shape[0] != a.shape[1]
).
- linalg.lstsq(a, b, rcond=None, weighted=False, extrainfo=False)¶
Least-squares solution
x
toa @ x = b
forgvar.GVar
s.Here
x
is defined to be the solution that minimizes||b - a @ x||
. Ifb
has a covariance matrix, another option is to weight the norm with the inverse covariance matrix: i.e., minimize|| isig @ b - isig @ a @ x||
whereisig
is the square root of the inverse ofb
’s covariance matrix. Set parameterweighted=True
to obtain the weighted-least-squares solution.- Parameters:
a – Matrix/array of shape
(M,N)
containing numbers and/orgvar.GVar
s.b – Vector/array of shape
(M,)
containing numbers and/orgvar.GVar
s.rcond (float) – Cutoff for singular values of
a
. Singular values smaller thanrcond
times the maximum eigenvalue are ignored. Default (rcond=None
) ismax(M,N)
times machine precision.weighted (bool) – If
True
, use weighted least squares; otherwise use unweighted least squares.extrainfo (bool) – If
False
(default) onlyx
is returned; otherwise(x, residual, rank, s)
is returned.
- Returns:
Array
x
of shape(N,)
that minimizes|| b - a @ x||
ifextrainfo==False
(default); otherwise returns a tuple(x, residual, rank, s)
whereresidual
is the sum of the squares ofb - a @ x
,rank
is the rank of matrixa
, ands
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.GVar
s. 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)
whereval[i]
are the eigenvalues ofa
, andvec[:, i]
are the mean values of the corresponding eigenvectors. Onlyval
is returned ifeigvec=False
(default).
- Returns:
Array
val
of eigenvalues of matrixa
if parametereigvec==False
(default); otherwise a tuple of arrays(val, vec)
whereval[i]
are the eigenvalues (in ascending order) andvec[:, 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.GVar
s. 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)
whereval[i]
are the eigenvalues ofa
(in ascending order), andvec[:, i]
are the corresponding eigenvectors ofa
. Onlyval
is returned ifeigvec=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
) isN
times machine precision, whereN
is the dimension of the array.
- Returns:
Tuple
(val,vec)
of eigenvalues and eigenvectors of matrixa
if parametereigvec==True
(default). The eigenvaluesval[i]
are in ascending order andvec[:, i]
are the corresponding eigenvalues. Only the eigenvaluesval
are returned ifeigvec=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
containinggvar.GVar
s.- Parameters:
a – Two-dimensional matrix/array of numbers and/or
gvar.GVar
s.compute_uv (bool) – It
True
(default), returns tuple(u,s,vT)
where matrixa = u @ np.diag(s) @ vT
where matricesu
andvT
satisfyu.T @ u = 1
andvT @ vT.T = 1
, ands
is the list of singular values. Onlys
is returned ifcompute_uv=False
.rcond (float) – Singular values whose difference is smaller than
rcond
times their sum are assumed to be degenerate for calculating variances foru
andvT
. Default (rcond=None
) ismax(M,N)
times machine precision.
- Returns:
Tuple
(u,s,vT)
where matrixa = u @ np.diag(s) @ vT
where matricesu
andvT
satisfyu.T @ u = 1
andvT @ vT.T = 1
, ands
is the list of singular values. Ifa.shape=(N,M)
, thenu.shape=(N,K)
andvT.shape=(K,M)
whereK
is the number of nonzero singular values (len(s)==K
). Ifcompute_uv==False
onlys
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
objectodeint
integratesdy/dx = f(x,y)
to obtainy(x1)
fromy0 = y(x0)
.y
andf(x,y)
can be scalars ornumpy
arrays. Typical usage is illustrated by the following code for integratingdy/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 fromx=0
tox=1
starting withy=y0
atx=0
; the result isy1=exp(1)
, of course. Similarly the second call toodeint
continues the integration fromx=1
tox=2
, givingy2=exp(2)
.If the
interval
is a list with more than two entries, thenodeint(y0, interval=[x0, x1, x2 ...])
in the example above returns an array of solutions for pointsx1, 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 functiony(x)
where: a)y(x0) = y0
; and b)y(x)
uses the integator to integrate the differential equation to pointx
starting from the last point at whichy
was evaluated (or fromx0
for the first call toy(x)
). The function can also be called with an array ofx
values, in which case an array containing the correspondingy
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 theIntegrator
by specifying parameterh
. A minimum step sizehmin
can also be specified; theIntegrator
raises an exception if the step size becomes smaller thanhmin
. TheIntegrator
keeps track of the number of good steps, whereh
is increased, and the number of bad steps, whereh
is decreased and the step is repeated:odeint.ngood
andodeint.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 erroryerr
after a given step, the proposed value fory
, and the proposed changedelta_y
iny
— that returns a number to compare with tolerancetol
. The step size is decreased and the step repeated ifdelta(yerr, y, delta_y) > tol
; otherwise the step is accepted and the step size increased. The default definition ofdelta
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 parameteranalyzer
. This function is called after every full step of the integration, with the current values ofx
andy
. Objects of typegvar.ode.Solution
are examples of (simple) analyzers.- Parameters:
deriv – Function of
x
andy
that returnsdy/dx
. The return value should have the same shape asy
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 toNone
.delta – Function
delta(yerr, y, delta_y)
that returns a number to be compared withtol
at each integration step: if it is larger thantol
, 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. Hereyerr
is an estimate of the error iny
on the last step;y
is the proposed value; anddelta_y
is the change iny
over the last step.analyzer – Function of
x
andy
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)
wherey
is a dictionary.An
DictIntegrator
objectodeint
integratesdy/dx = f(x,y)
to obtainy(x1)
fromy0 = y(x0)
.y
andf(x,y)
are dictionary types having the same keys, and containing scalars and/ornumpy
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 fromx=x0
tox=x1
, returningy1=y(x1)
. The second call continues the integration tox=x2
, returningy2=y(x2)
. Multiple integration points can be specified ininterval
, in which case a list of the correspondingy
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 theDictIntegrator
by specifying parameterh
. A minimum ste psizehmin
can also be specified; theIntegrator
raises an exception if the step size becomes smaller thanhmin
. TheDictIntegrator
keeps track of the number of good steps, whereh
is increased, and the number of bad steps, whereh
is decreases and the step is repeated:odeint.ngood
andodeint.nbad
, respectively.An analyzer
analyzer(x,y)
can be specified using parameteranalyzer
. This function is called after every full step of the integration with the current values ofx
andy
. Objects of typegvar.ode.Solution
are examples of (simple) analyzers.- Parameters:
deriv – Function of
x
andy
that returnsdy/dx
. The return value should be a dictionary with the same keys asy
, and values that have the same shape as the corresponding values iny
.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 toNone
.delta – Function
delta(yerr, y, delta_y)
that returns a number to be compared withtol
at each integration step: if it is larger thantol
, 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. Hereyerr
is an estimate of the error iny
on the last step;y
is the proposed value; anddelta_y
is the change iny
over the last step.analyzer – Function of
x
andy
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, andsoln.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 callresult = integral(fcn, interval=(x0, x1))
calculates the integral of
fcn(x)
fromx0
tox1
. 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. Settingfshape=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
forGVar
s.The
order=(m,n)
Pade approximant to a series given bysum_i f[i] * x**i
is the ratio of polynomials of orderm
(numerator) andn
(denominator) whose Taylor expansion agrees with that of the original series up to orderm+n
.A
Pade
objectpade
createsgvar.powerseries.PowerSeries
objects for the numerator (pade.num
) and the denominator (pade.den
) of the Pade approximant corresponding to the input seriesf[i]
. The approximant can be evaluated for arbitraryx
usingpade(x)
. The coefficients used in the numerator and denominator are given bypade.num.c
andpade.den.c
, respectively.Elements in the series
f[i]
may be numbers orgvar.GVar
s. When the latter appear, the code uses an SVD algorithm (seepade_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 bypade.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 fori=0...n+m
.order (tuple) –
order=(m,n)
specifies that the numerator and denominator of the Pade approximant have maximum orderm
andn
, 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 thef[i]
. Setrtol
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 inf[i]
are ignored. The SVD analysis is skipped ifrtol=None
.
- gvar.pade.pade_gvar(f, m, n, rtol='gavg')¶
(m,n)
Pade approximant tosum_i f[i] x**i
forGVar
s.The
(m,n)
Pade approximant to a series given bysum_i f[i] * x**i
is the ratio of polynomials of orderm
(numerator) andn
(denominator) whose Taylor expansion agrees with that of the original series up to orderm+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 fori=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 thef[i]
. Setrtol
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 inf[i]
are ignored. The SVD analysis is skipped ifrtol=None
.
- Returns:
Tuple of power series coefficients
(p, q)
such thatsum_i p[i] x**i
is the numerator of the approximant, andsum_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 tosum_i f[i] x**i
.The
(m,n)
Pade approximant to a series given bysum_i f[i] * x**i
is the ratio of polynomials of orderm
(numerator) andn
(denominator) whose Taylor expansion agrees with that of the original series up to orderm+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 coefficientsf[i]
. It automatically reduces the order of the approximant if the extraction of Pade coefficients is too unstable given tolerancertol
.- Parameters:
f – Array
f[i]
of power series coefficients fori=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 thatsum_i p[i] x**i
is the numerator of the approximant, andsum_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 theorder
parameter. The list ofc[i]
s is then padded with zeros ifc
is too short, or truncated if it is too long. Omittingc
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 ofPowerSeries p
isp[3]
. The order ofp
isp.order
.PowerSeries
should work for coefficients of any data type that supports ordinary arithmetic.Arithmetic expressions of
PowerSeries
variables yield newPowerSeries
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 arrayc
is too small for the specifiedorder
, 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 ofself
.
- 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 ofself
.
- 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. Ind
dimensions,c[i1,i2,...,id]
is the coefficient ofx1**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
. IfNone
, 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 adim
-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 offcn(x)=0
. It examines pointsx[j + 1] = fac * x[j] + incr
where
x[0]=x0
andj=0...maxit-1
, looking for a pair of successive points wherefcn(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 usingroot.refine()
.Example
The following code seeks to bracket the first zero of
sin(x)
withx>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 pointx
that is examined. This can be used, for example, to monitor the search while debugging. Default isNone
.
- Returns:
Tuple
(a, b)
wherefcn(a) * fcn(b) <= 0
, which implies that a root occurs betweena
andb
(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 functionfcn
on an interval.This method finds a root
x
offcn(x)=0
inside aninterval=(a,b)
that brackets the root, withfcn(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 interval1 <= 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, withfcn(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 setsrtol
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 pointx
examined by the algorithm. This can be used, for example, to monitor convergence while debugging. Default isNone
.
- Returns:
The root, which is either a
float
or agvar.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, whereb
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.