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

## 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 `xknot`s 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.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 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.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 that `a @ x = b` for matrix `a`.

  Parameters:

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

    - b – One-dimensional vector/array of numbers and/or `gvar.GVar`s, 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.GVar`s.

  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.GVar`s.

    - b – Vector/array of shape `(M,)` containing numbers and/or `gvar.GVar`s.

    - 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.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)` 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.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)` 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.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 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 `GVar`s.

  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.GVar`s. 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 `GVar`s.

  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.