Module scid.calculus

Functions for numerical integration (quadrature) and differentiation.

It is recommended to read the following introduction before using any of the functions in this module, as it discusses features and issues that are common to several of them.

Integration

This module contains several different integration methods, suitable for a wide range of problems. The choice of which method to use depends on several things, such as whether the integration interval is finite or infinite, whether the integrand oscillates, or whether it exhibits other difficult behaviour such as singularities or discontinuities in the integration interval. Choosing the right method can both speed up the calculation and make it more accurate. Therefore, each method below has a short description of which types of integral it is suitable for.

When you just want to "get things done" without any fuss, use the general-purpose integrate() function. It tries to select the method which is most likely to succeed for the given interval.

All the integration functions have a set of input parameters in common: The function to integrate, the limits of the integral, and some accuracy requirements. All of these are described in the documentation for integrate().

Most of the quadrature routines in SciD are ported from from QUADPACK which is a well-known FORTRAN package for numerical evaluation of integrals using Gaussian quadrature. Gaussian quadrature is considered to be the most accurate method available for numerical integration of smooth functions. The QUADPACK integrators are named integrateQ*() below.

In addition, SciD has support for double exponential integration. The DE algorithms are ported from Takuya Ooura's DE-Quadrature package, and are available in this module through the integrateDE*() functions.

Differentiation

The simplest, but fastest, form of numerical differentiation is using finite differences. The functions in this module use three different finite-difference formulas: Forward differences,

df   f(x+h) - f(x)
-- = -------------
dx         h

backward differences,

df   f(x) - f(x-h)
-- = -------------
dx         h

and central differences,

df   f(x+h) - f(x-h)
-- = ---------------
dx         2 h

Of these three, the latter is the most accurate, but in the cases where f(x) is already known it requires one more function evaluation compared to the forward-/backward-difference formulas. This is of particular importance when one needs to evaluate several derivatives, as is the case with e.g. gradients and Jacobians.

Some functions in this module, like diff(), use another method by Ridders, which extrapolates the results of finite-difference formulas like the above to make the approximation more accurate – usually a lot more so. This requires several function evaluations and is therefore also a lot slower than simple finite differences.

All the differentiation methods in this module take an optional scale parameter. This is a "characteristic scale" for the function, i.e. a scale over which the function changes significantly. It is worth experimenting a bit with this parameter, as it can have a drastic impact on accuracy. (For an extreme example, see the "Example" section for the diff() function below.) For practical reasons, 1.0 is chosen as the default scale in this module. Another, probably more common, choice that sometimes works well is to set scale = x, where x is the point at which the derivative is taken.

Functions

NameDescription
checkQuadpackStatus(ier, result, abserr, file, line)
diff(f, x, scale, tableauSize)Calculate the derivative of a function.
diff2(f, x, fx, scale)Calculate the derivative of a function using a two-point formula, i.e. a forward- or backward-difference formula.
diff3(f, x, scale)Calculate the derivative of a function using a three-point formula, a.k.a. a central difference formula.
gradient(f, x, scale, buffer)Calculate the gradient of a function of several variables.
gradientR(f, x, scale, buffer)Calculate the gradient of a function of several variables using Ridders' method.
hessian(f, x, scale, buffer)Calculate the Hessian matrix of a function of several variables using a central-difference approximation.
integrate(f, a, b, epsRel, epsAbs)Integrate a function from a to b.
integrateDE(f, a, b, epsRel)Calculate the integral of f(x) over the finite interval (a,b) using double exponential integration.
integrateDEI(f, a, epsRel)Calculate the integral of f(x) over the infinite interval (a,∞) using double exponential integration.
integrateDEO(f, a, omega, epsRel)Calculate the integral of an oscillating function f(x) over the infinite interval (a,∞) using double exponential integration.
integrateQAG(f, a, b, rule, epsRel, epsAbs)Calculate the integral of f(x) over the finite interval (a,b) using a simple globally adaptive integrator.
integrateQAGI(f, epsRel, epsAbs)Calculate the integral of f(x) over an infinite interval.
integrateQAGP(f, a, b, trouble, epsRel, epsAbs)Calculate the integral of f(x) over the finite interval (a,b), taking into account known points of special difficulty inside the interval.
integrateQAGS(f, a, b, epsRel, epsAbs)Calculate the integral of f(x) over the finite interval (a,b) using a general-purpose integration algorithm.
integrateQAWC(f, a, b, c, epsRel, epsAbs)Calculate a Cauchy principal value integral.
integrateQAWF(f, a, omega, weight, epsAbs)Calculates a Fourier transform integral.
integrateQAWO(f, a, b, omega, weight, epsRel, epsAbs)Calculate the integral of an oscillatory function over the finite interval (a,b).
integrateQAWS(f, a, b, alpha, beta, weight, epsRel, epsAbs)Calculate an integral over the finite interval (a,b), where the integrand has algebraic and/or logarithmic endpoint singularities of a known type.
integrateQNG(f, a, b, epsRel, epsAbs)Calculate the integral of f(x) over the finite interval (a,b) using a simple non-adaptive automatic integrator, based on a sequence of rules with increasing degree of algebraic precision.
jacobian(f, m, x, scale, buffer)Calculate the Jacobian matrix associated with a set of m functions in n variables using a central-difference approximation to the Jacobian.
jacobian2(f, m, x, scale, fx, buffer)Calculate the Jacobian using forward-/backward-difference methods, also known as 2-point formulas.

Classes

NameDescription
IntegrationExceptionException thrown when an integration routine fails.

Enums

NameDescription
GaussKronrod
Infinite
Oscillation
Weight