DAMath: Numerical Quadrature

DAMath: Numerical Quadrature#

For a general overview, see: Wikipedia [1187], Wikipedia [1194], MathWorld [874], Wikipedia [1174].

General error codes for Amath integration functions, type 1#

Amath integration functions which are NOT using the double exponential transformation return an error code ier, which has the following meaning:

ier = 0: Normal and reliable termination of the routine. It is assumed that the requested accuracy has been achieved.

ier > 0: Abnormal termination of the routine the estimates for integral and error are less reliable. It is assumed that the requested accuracy has not been achieved.

ier = 1: Maximum number of subdivisions allowed has been achieved. One can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). However, if this yields no improvement it is advised to analyze the integrand in order to determine the integration difficulties. If the position of a local difficulty can be determined (e.g. singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling the integrator on the subranges. If possible, an appropriate special-purpose integrator should be used, which is designed for handling the type of difficulty involved.

ier = 2: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be under-estimated.

ier = 3: Extremely bad integrand behaviour occurs at some points of the integration interval.

ier = 4: The algorithm does not converge. Roundoff error is detected in the extrapolation table. It is presumed that the requested tolerance cannot be achieved, and that the returned result is the best which can be obtained.

ier = 5: The integral is probably divergent, or slowly convergent. It must be noted that divergence can occur with any other value of ier.

ier = 6: The input is invalid, because epsabs <= 0 and epsrel < 50*eps_d. result, abserr, last are set to zero.

ier = 7: The input is invalid, limit < 0 or limit > QMAXLIM. result, abserr, last are set to zero.

ier = 8: Dynamic list vectors cannot be allocated. result, abserr, last are set to zero.

ier = 9: At least one limit a or b is NaN, or infinite a=b. result, abserr, last are set to NaN_d.

General error codes for Amath integration functions, type 2#

Amath integration functions which are using the double exponential transformation return an error code ier, which has the following meaning:

ier = 0: Normal and reliable termination of the routine. It is assumed that the requested accuracy has been achieved.

ier > 0: Abnormal termination of the routine the estimates for integral and error are less reliable. It is assumed that the requested accuracy has not been achieved.

ier = 1: eps < 8*eps_x.

ier = 2: roundoff problems.

ier = 3: max. iterations, result/abserr have values.

Global adaptive quadrature by Forsythe, Malcolm, Moler (quanc8)#

ctx.quanc8(f, a, b, epsabs, epsrel, limit)#

where ctx is math53. Performs a global adaptive quadrature of f over (a,b) based on a Fortran subroutine by Forsythe, Malcolm, Moler.

Parameters:

  • f – the function for which the integral is determined.

  • a – the left border of the integration interval.

  • b – the right border of the integration interval.

  • abserr – the absolute accuracy requested.

  • relerr – the relative accuracy requested.

Returns:

a tuple (integral, err, flag, neval), where integral is an approximation to the integral; abserr is an estimate of of the modulus of the absolute error; neval is the number of integrand evaluations; and ier is the error code (success=0, error > 0: see General error codes for Amath integration functions, type 1).

Include reference to Forsythe, Malcolm, Moler.

The routine in Python:

>>> from xlcalcnet import xreal
>>> a = 0.0; b = 2.0; epsabs = 1E-8; epsrel = 0.0; limit = 0
>>> f = lambda x: xreal.Exp(-x * x / 2)
>>> result, abserr, neval, ier = xreal.Qags(f, a, b, epsabs, epsrel, limit)
>>> print("result: ", result, "abserr: ", abserr)
>>> print("neval: ", neval, "ier: ", ier)

The same routine in Visual Basic:

Function XF4(x As xreal) As xreal
    Dim y As New xreal
    y = xreal.StudentTPdf(10, x)
    'Console.WriteLine("x : {0}, y: {1}", x, y)
    Return y
End Function

Sub DemoQags()
    Dim a, b, epsabs, epsrel As New xreal, limit As Int32=0
    a = "0" : b = "2" : epsabs = "1E-15" : epsrel = "1E-15"
    Dim OutputFormat = "Qags: f = F4, a = {0}, b = {1}, epsabs = {2}, epsrel = {3}, limit = {4}"
    Console.WriteLine(OutputFormat, a, b, epsabs, epsrel, limit)
    Dim Res2 = xreal.Qags(AddressOf XF4, a, b, epsabs, epsrel)
    Console.WriteLine("Res2:(result, abserr, neval, ier) {0}", Res2)
    Console.WriteLine()
End Sub

21-point Gauss-Kronrod rule, finite interval (qags)#

ctx.qags(f, a, b, epsabs, epsrel, limit)#

where ctx is math53. Performs a global adaptive quadrature of f over (a,b) based on 21-point Gauss-Kronrod rule for the subintervals, with acceleration by Wynn’s epsilon algorithm.

Parameters:

  • f – the function for which the integral is determined.

  • a – the left border of the integration interval.

  • b – the right border of the integration interval.

  • epsabs – the absolute accuracy requested.

  • epsrel – the relative accuracy requested.

  • limit – upperbound on the no. of subintervals, 0: use DefLimit.

Returns:

a tuple (integral, abserr, neval, ier), where integral is an approximation to the integral; abserr is an estimate of of the modulus of the absolute error; neval is the number of integrand evaluations; and ier is the error code (success=0, error > 0: see General error codes for Amath integration functions, type 1).

See also: Piessens et al. [495], Wikipedia [1188], Wikipedia [1202].

The routine in Python:

>>> from xlcalcnet import xreal
>>> a = 0.0; b = 2.0; epsabs = 1E-8; epsrel = 0.0; limit = 0
>>> f = lambda x: xreal.Exp(-x * x / 2)
>>> result, abserr, neval, ier = xreal.Qags(f, a, b, epsabs, epsrel, limit)
>>> print("result: ", result, "abserr: ", abserr)
>>> print("neval: ", neval, "ier: ", ier)

The same routine in Visual Basic:

Function XF4(x As xreal) As xreal
    Dim y As New xreal
    y = xreal.StudentTPdf(10, x)
    'Console.WriteLine("x : {0}, y: {1}", x, y)
    Return y
End Function

Sub DemoQags()
    Dim a, b, epsabs, epsrel As New xreal, limit As Int32=0
    a = "0" : b = "2" : epsabs = "1E-15" : epsrel = "1E-15"
    Dim OutputFormat = "Qags: f = F4, a = {0}, b = {1}, epsabs = {2}, epsrel = {3}, limit = {4}"
    Console.WriteLine(OutputFormat, a, b, epsabs, epsrel, limit)
    Dim Res2 = xreal.Qags(AddressOf XF4, a, b, epsabs, epsrel)
    Console.WriteLine("Res2:(result, abserr, neval, ier) {0}", Res2)
    Console.WriteLine()
End Sub

15-point Gauss-Kronrod rule, infinite interval (qagi)#

math53.qagi(f, bound, inf, epsabs, epsrel, limit)#

where ctx is math53. Performs a global adaptive quadrature of f over an infinite interval based on a transformed 15-point Gauss-Kronrod for the subintervals, with acceleration by Wynn’s epsilon algorithm.

Parameters:

  • f – the function for which the integral is determined.

  • bound – finite bound of integration range (if any).

  • inf – indicating the kind of integration range involved: 1 corresponds to (bound, +infinity), -1 to (-infinity, bound), 2 to (-infinity, +infinity).

  • epsabs – the absolute accuracy requested.

  • epsrel – the relative accuracy requested.

  • limit – upperbound on the no. of subintervals, 0: use DefLimit.

Returns:

a tuple (integral, abserr, neval, ier), where integral is an approximation to the integral; abserr is an estimate of of the modulus of the absolute error; neval is the number of integrand evaluations; and ier is the error code (success=0, error > 0: see General error codes for Amath integration functions, type 1).

See also: Piessens et al. [495], Wikipedia [1188].

The routine in Python:

>>> from xlcalcnet import xreal
>>> bound = 0.0; infcode = 2; epsabs = 1E-8; epsrel = 0.0; limit = 0
>>> f = lambda x: xreal.Exp(-x * x / 2)
>>> result, abserr, neval, ier = xreal.Qagi(f, bound, infcode, epsabs, epsrel, limit)
>>> print("result: ", result, "abserr: ", abserr)
>>> print("neval: ", neval, "ier: ", ier)

The same routine in Visual Basic:

Function XF4(x As xreal) As xreal
    Dim y As New xreal
    y = xreal.StudentTPdf(10, x)
    'Console.WriteLine("x : {0}, y: {1}", x, y)
    Return y
End Function

Sub DemoQagi()
    Dim bound, epsabs, epsrel As New xreal, inf, limit As Short
    bound = "0": inf = 1: epsabs = "1E-15" : epsrel = "1E-15": limit = 0
    Dim OutputFormat = "Qagi: f = XF4, a = {0}, b = {1}, epsabs = {2}, epsrel = {3}, limit = {4}"
    Console.WriteLine(OutputFormat, bound, inf, epsabs, epsrel, limit)
    Dim Res3 = xreal.Qagi(AddressOf XF4, bound, inf, epsabs, epsrel)
    Console.WriteLine("Res3:(result, abserr, neval, ier) {0}", Res3)
    Console.WriteLine()
End Sub

Cauchy principal value, finite interval (qawc)#

math53.qawc(f, a, b, c, epsabs, epsrel, limit)#

where ctx is math53. Performs an adaptive quadrature of \(f(x)/(x-c)\) over the finite interval \((a,b)\) with the singularity at \(c\) and \(c\) not equal \(a\) or \(b\). The routine calculates an approximation result to the Cauchy principal value.

Parameters:

  • f – the function for which the integral is determined.

  • a – the left border of the integration interval.

  • b – the right border of the integration interval.

  • c – singularity, ier=6 if \(c=a\) or \(c=b\).

  • epsabs – the absolute accuracy requested.

  • epsrel – the relative accuracy requested.

  • limit – upperbound on the no. of subintervals, 0: use DefLimit.

Returns:

a tuple (integral, abserr, neval, ier), where integral is an approximation to the integral; abserr is an estimate of of the modulus of the absolute error; neval is the number of integrand evaluations; and ier is the error code (success=0, error > 0: see General error codes for Amath integration functions, type 1).

See also: https://maxima.sourceforge.io/docs/manual/maxima_100.html

See also: Piessens et al. [495], Wikipedia [1188], Wikipedia [1197].

The routine in Python:

# Qawc Example is missing
bound = 0.0; infcode = 2; epsabs = 1E-8; epsrel = 0.0; limit = 0
f = lambda x: xreal.Exp(-x * x / 2)
result, abserr, neval, ier = xreal.Qagi(f, bound, infcode, epsabs, epsrel, limit)
print("xreal.Qagi")
print("result: ", result, "abserr: ", abserr)
print("neval: ", neval, "ier: ", ier)
print()

The same routine in Visual Basic:

Function XF5(x As xreal) As xreal
    Dim y As New xreal
    y = 1.0/(1 * (5*x*x*x + 6))
    'Console.WriteLine("x : {0}, y: {1}", x, y)
    Return y
End Function

Sub DemoQawc()
    Dim a, b, c, epsabs, epsrel As New xreal, limit As Int32=0
    a = "-1" : b = "5" : c = "0": epsabs = "1E-5" : epsrel = "1E-5"
    Dim OutputFormat = "Qawc: f = XF5, a = {0}, b = {1}, c = {2}, epsabs = {3}, epsrel = {4}, limit = {5}"
    Console.WriteLine(OutputFormat, a, b, c, epsabs, epsrel, limit)
    Dim Res4 = xreal.Qawc(AddressOf XF5, a, b, c, epsabs, epsrel)
    Console.WriteLine("Res4:(result, abserr, neval, ier) {0}", Res4)
End Sub

Double Exponential (DE) transformation, finite interval (intde)#

math53.intde(f, a, b, eps)#

where ctx is math53. Performs an automatic quadrature of f(x) over the finite interval (a,b)} using Double Exponential (DE) transformation.

Parameters:

  • f – the function for which the integral is determined.

  • a – the left border of the integration interval.

  • b – the right border of the integration interval.

  • eps – the relative accuracy requested.

Returns:

a tuple (integral, abserr, neval, ier), where integral is an approximation to the integral, if ier=0 or 3; abserr is an estimate of of the modulus of the absolute error, if ier=0 or 3; neval is the number of integrand evaluations; and ier is the error code (success=0, error > 0: see General error codes for Amath integration functions, type 2).

See also: Wikipedia [1211], Takahasi and Mori [537], Ooura and M. [472], Mori [438].

The routine in Python:

>>> from xlcalcnet import xreal
>>> a = 0.0; b = 2.0; eps = 1E-8
>>> f = lambda x: xreal.Exp(-x * x / 2)
>>> result, abserr, neval, ier = xreal.Intde(f, a, b, eps)
>>> print("result: ", result, "abserr: ", abserr)
>>> print("neval: ", neval, "ier: ", ier)

The same routine in Visual Basic:

Function XF4(x As xreal) As xreal
    Dim y As New xreal
    y = xreal.StudentTPdf(10, x)
    'Console.WriteLine("x : {0}, y: {1}", x, y)
    Return y
End Function

Sub DemoIntde()
    Dim a, b, epsabs As New xreal, limit As Int32=0
    a = "0" : b = "2" :epsabs = "1E-12"
    Console.WriteLine("Intde: f = XF4, a = {0}, b = {1}, epsabs = {2}", a, b, epsabs)
    Dim Res5 = xreal.Intde(AddressOf XF4, a, b, epsabs)
    Console.WriteLine("Res5:(result, abserr, neval, ier) {0}", Res5)
End Sub

DE transformation, infinite interval, no oscillatory factor (intdei)#

math53.intdei(f, a, eps)#

where ctx is math53. Performs an automatic quadrature of f(x) over (a,INF) using Double Exponential transformation when f(x) has no oscillatory factor.

Parameters:

  • f – the function for which the integral is determined.

  • a – the left border of the integration interval.

  • eps – the relative accuracy requested.

Returns:

a tuple (integral, abserr, neval, ier), where integral is an approximation to the integral, if ier=0 or 3; abserr is an estimate of of the modulus of the absolute error, if ier=0 or 3; neval is the number of integrand evaluations; and ier is the error code (success=0, error > 0: see General error codes for Amath integration functions, type 2).

See also: Wikipedia [1211], Takahasi and Mori [537], Ooura and M. [472], Mori [438].

The routine in Python:

>>> from xlcalcnet import xreal
>>> a = 0.0; eps = 1E-8
>>> f = lambda x: xreal.Exp(-x * x / 2)
>>> result, abserr, neval, ier = xreal.Intdei(f, a, eps)
>>> print("result: ", result, "abserr: ", abserr)
>>> print("neval: ", neval, "ier: ", ier)

The same routine in Visual Basic:

Function XF4(x As xreal) As xreal
    Dim y As New xreal
    y = xreal.StudentTPdf(10, x)
    'Console.WriteLine("x : {0}, y: {1}", x, y)
    Return y
End Function

Sub DemoIntdei()
    Dim a, epsabs As New xreal, limit As Int32=0
    a = "0" : epsabs = "1E-15"
    Console.WriteLine("Intdei: f = XF4, a = {0}, epsabs = {1}", a, epsabs)
    Dim Res6 = xreal.Intdei(AddressOf XF4, a, epsabs)
    Console.WriteLine("Res6:(result, abserr, neval, ier) {0}", Res6)
    Console.WriteLine()
End Sub

DE transformation, infinite interval, oscillatory factor (intdeo)#

math53.intdeo(f, a, omega, eps)#

where ctx is math53. Performs an automatic quadrature of f(x) over (a,INF) using Double Exponential transformation when f(x) has an oscillatory factor.

Parameters:

  • f – the function for which the integral is determined.

  • a – the left border of the integration interval.

  • omega – oscillatory factor.

  • eps – the relative accuracy requested.

Returns:

a tuple (integral, abserr, neval, ier), where integral is an approximation to the integral, if ier=0 or 3; abserr is an estimate of of the modulus of the absolute error, if ier=0 or 3; neval is the number of integrand evaluations; and ier is the error code (success=0, error > 0: see General error codes for Amath integration functions, type 2).

See also: Wikipedia [1211], Takahasi and Mori [537], Ooura and M. [472], Mori [438].

The routine in Python:

# Intdeo example is missing
a = 0.0; eps = 1E-8
f = lambda x: xreal.Exp(-x * x / 2)
result, abserr, neval, ier = xreal.Intdei(f, a, eps)
print("xreal.Intdei")
print("result: ", result, "abserr: ", abserr)
print("neval: ", neval, "ier: ", ier)
return "0"

The same routine in Visual Basic:

Function XF6(x As xreal) As xreal
    Dim y, alpha As New xreal
    alpha = "2.0"
    y = xreal.Sin(x*alpha)/xreal.Sqrt(x)
    'Console.WriteLine("x : {0}, y: {1}", x, y)
    Return y
End Function

Sub DemoIntdeo()
    Dim a, alpha, epsabs As New xreal, limit As Int32=0
    a = "0" : alpha = "2" :epsabs = "1E-16"
    Console.WriteLine("Intdeo: f = XF6, a = {0}, alpha = {1}, epsabs = {2}", a, alpha, epsabs)
    Dim Res7 = xreal.Intdeo(AddressOf XF6, a, alpha, epsabs)
    Console.WriteLine("Res7:(result, abserr, neval, ier) {0}", Res7)
    Console.WriteLine("Analytic: {0}", Math.Sqrt(0.5*Math.PI/2))
    Console.WriteLine()
End Sub
Contents