Module scid.calculus

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.