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 generalpurpose 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 wellknown 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
DEQuadrature 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 finitedifference formulas: Forward differences,
df f(x+h)  f(x)
 = 
dx h
backward differences,
df f(x)  f(xh)
 = 
dx h
and central differences,
df f(x+h)  f(xh)
 = 
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/backwarddifference 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 finitedifference 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
Name  Description 

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 twopoint formula, i.e. a forward or backwarddifference formula. 
diff3(f, x, scale)  Calculate the derivative of a function using a threepoint 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 centraldifference 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 generalpurpose 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
nonadaptive 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 centraldifference approximation to the Jacobian. 
jacobian2(f, m, x, scale, fx, buffer)  Calculate the Jacobian using forward/backwarddifference methods, also known as 2point formulas. 
Classes
Name  Description 

IntegrationException  Exception thrown when an integration routine fails. 
Enums
Name  Description 

GaussKronrod  
Infinite  
Oscillation  
Weight 