src/numericalnim/integrate

Search:
Group by:
Source   Edit  

Integration

This module implements various integration routines. It provides:

Integrate discrete data:

  • trapz, simpson: works for any spacing between points.
  • romberg: requires equally spaced points and the number of points must be of the form 2^k + 1 ie 3, 5, 9, 17, 33, 65, 129 etc.

Example:

import src/numericalnim/integrate
import numericalnim, std/[sequtils]
let x = linspace(0.0, 5.0, 5)
let y = x.mapIt(it*it/2)

let integral = simpson(y, x)

Integrate continous functions:

  • adaptiveGauss (recommended): Adaptive Gaussian quadrature. This is often the best choice as it is both fast and accurate. It also handles infinite integration limits.
  • gaussQuad: Fixed step size Gaussian quadrature.
  • romberg: Adaptive method based on Richardson Extrapolation.
  • adaptiveSimpson: Adaptive step size.
  • simpson: Fixed step size.
  • trapz: Fixed step size.

Example:

import src/numericalnim/integrate
import numericalnim, std/math

proc f(x: float, ctx: NumContext[float, float]): float =
    exp(x)

let a = 0.0
let b = Inf
let integral = adaptiveGauss(f, a, b)

Cumulative integration:

  • Discrete: cumtrapz, cumsimpson.
  • Continous:
    • cumGauss: adaptive step size.
    • cumGaussSpline: returns a HermiteSpline constructed from the cumulative integral.
    • cumsimpson, cumtrapz: fixed step size.

Example:

import src/numericalnim/integrate
import numericalnim, std/[sequtils]
# Discrete cumulative
let x = linspace(0.0, 5.0, 5)
let y = x.mapIt(it*it/2)

let cumintegral = cumsimpson(y, x)

Example:

import src/numericalnim/integrate
import numericalnim, std/math
# Continous cumulative
proc f(x: float, ctx: NumContext[float, float]): float =
    exp(x)

# The values we want the cumulative integral at:
let xspan = linspace(0.0, 5.0, 10)
let cumintegral = cumGauss(f, xspan)

Procs

proc adaptiveGauss[T; U](f_in: NumContextProc[T, U]; xStart_in, xEnd_in: U;
                         tol = 1e-8; initialPoints: openArray[U] = @[];
                         maxintervals: int = 10000; ctx: NumContext[T, U] = nil): T

Calculate the integral of f using an globally adaptive Gauss-Kronrod Quadrature. Inf and -Inf can be used as integration limits.

Input:

  • f: the function that is integrated.
  • xStart: The start of the integration interval.
  • xEnd: The end of the integration interval.
  • tol: The error tolerance that must be satisfied on every subinterval.
  • maxintervals: maximum numbers of intervals to divide integral in before stopping.
  • initialPoints: A list of known difficult points (integrable singularities, discontinouities etc) that will be used as the inital interval boundaries.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • The value of the integral of f from xStart to xEnd calculated using an adaptive Gauss-Kronrod Quadrature.
Source   Edit   
proc adaptiveGaussLocal[T](f: NumContextProc[T, float]; xStart, xEnd: float;
                           tol = 1e-8; ctx: NumContext[T, float] = nil): T

Calculate the integral of f using an locally adaptive Gauss-Kronrod Quadrature.

Input:

  • f: the function that is integrated.
  • xStart: The start of the integration interval.
  • xEnd: The end of the integration interval.
  • tol: The error tolerance that must be satisfied on every subinterval.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • The value of the integral of f from xStart to xEnd calculated using an adaptive Gauss-Kronrod Quadrature.
Source   Edit   
proc adaptiveSimpson[T](f: NumContextProc[T, float]; xStart, xEnd: float;
                        tol = 1e-8; ctx: NumContext[T, float] = nil): T

Calculate the integral of f using an adaptive Simpson's rule.

Input:

  • f: the function that is integrated.
  • xStart: The start of the integration interval.
  • xEnd: The end of the integration interval.
  • tol: The error tolerance that must be satisfied on every subinterval.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • The value of the integral of f from xStart to xEnd calculated using an adaptive Simpson's rule.
Source   Edit   
proc adaptiveSimpson2[T](f: NumContextProc[T, float]; xStart, xEnd: float;
                         tol = 1e-8; ctx: NumContext[T, float] = nil): T

Calculate the integral of f using an adaptive Simpson's rule.

Input:

  • f: the function that is integrated.
  • xStart: The start of the integration interval.
  • xEnd: The end of the integration interval.
  • tol: The error tolerance that must be satisfied on every subinterval.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • The value of the integral of f from xStart to xEnd calculated using an adaptive Simpson's rule.
Source   Edit   
proc cmpInterval[T; U; V](interval1, interval2: IntervalType[T, U, V]): int
Source   Edit   
proc cumGauss[T](f_in: NumContextProc[T, float]; X: openArray[float];
                 tol = 1e-8; initialPoints: openArray[float] = @[];
                 maxintervals: int = 10000; ctx: NumContext[T, float] = nil): seq[
    T]

Calculate the cumulative integral of f using an globally adaptive Gauss-Kronrod Quadrature. Returns a sequence of values which is the cumulative integral of f at the points defined in X. Important: because of the much higher order of the Gauss-Kronrod quadrature (order 21) compared to the interpolating Hermite spline (order 3) you have to give it a large amount of initialPoints. Otherwise it may only use a couple of points which gives quite a bad interpolant. By default if no initial points are given, 100 uniformly spaced points are used. The more points the better interpolant but the longer it will take to run the integration.

Input:

  • f: the function that is integrated.
  • X: the points the cumulative integral should be evaluated at.
  • tol: The error tolerance that must be satisfied on every subinterval.
  • maxintervals: maximum numbers of intervals to divide integral in before stopping.
  • initialPoints: A list of known difficult points (integrable singularities, discontinouities etc) that will be used as the inital interval boundaries. This is also the minimum number of points it will evaluate f at, if the function is too smooth it may be evaluated in too few points to give a good enough interpolation.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • A sequence of values which is the cumulative integral of f at the points defined in X.
Source   Edit   
proc cumGaussSpline[T; U](f_in: NumContextProc[T, U]; xStart_in, xEnd_in: U;
                          tol = 1e-8; initialPoints: openArray[U] = @[];
                          maxintervals: int = 10000; ctx: NumContext[T, U] = nil): InterpolatorType[
    T]

Calculate the cumulative integral of f using an globally adaptive Gauss-Kronrod Quadrature. Inf and -Inf can be used as integration limits. Returns a Hermite spline that can be evaluated at any point between xStart and xEnd. Important: because of the much higher order of the Gauss-Kronrod quadrature (order 21) compared to the interpolating Hermite spline (order 3) you have to give it a large amount of initialPoints. Otherwise it may only use a couple of points which gives quite a bad interpolant. By default if no initial points are given, 100 uniformly spaced points are used. The more points the better interpolant but the longer it will take to run the integration.

Input:

  • f: the function that is integrated.
  • xStart: The start of the integration interval.
  • xEnd: The end of the integration interval.
  • tol: The error tolerance that must be satisfied on every subinterval.
  • maxintervals: maximum numbers of intervals to divide integral in before stopping.
  • initialPoints: A list of known difficult points (integrable singularities, discontinouities etc) that will be used as the inital interval boundaries. This is also the minimum number of points it will evaluate f at, if the function is too smooth it may be evaluated in too few points to give a good enough interpolation.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • A 1D spline which represents the cumulative integral of f from xStart to xEnd.
Source   Edit   
proc cumsimpson[T](f: NumContextProc[T, float]; X: openArray[float];
                   ctx: NumContext[T, float] = nil; dx = 0.00001): seq[T]

Calculate the cumulative integral of f using Simpson's rule.

Input:

  • f: the function that is integrated.
  • X: The x-values of the returned values.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.
  • dx: The step length to use when integrating.

Returns:

  • The value of the integral of f from the smallest to the largest value of X calculated using Simpson's rule.
Source   Edit   
proc cumsimpson[T](Y: openArray[T]; X: openArray[float]): seq[T]

Calculate the cumulative integral of f using Simpson's rule from a set of values.

Input:

  • Y: A seq of values of the integrand.
  • X: A seq with the corresponding x-values.

Returns:

  • The cumulative integral evaluated from the smallest to the largest value in X calculated using Simpson's rule.
Source   Edit   
proc cumtrapz[T](f: NumContextProc[T, float]; X: openArray[float];
                 ctx: NumContext[T, float] = nil; dx = 0.00001): seq[T]

Calculate the cumulative integral of f using the trapezoidal rule at the points in X.

Input:

  • f: the function that is integrated.
  • X: The x-values of the returned values.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.
  • dx: The step length to use when integrating.

Returns:

  • The value of the integral of f from the smallest to the largest value of X calculated using the trapezoidal rule.
Source   Edit   
proc cumtrapz[T](Y: openArray[T]; X: openArray[float]): seq[T]

Calculate the cumulative integral of f using the trapezoidal rule from a set of values.

Input:

  • Y: A seq of values of the integrand.
  • X: A seq with the corresponding x-values.

Returns:

  • The cumulative integral evaluated from the smallest to the largest value in X calculated using the trapezoidal rule.
Source   Edit   
proc gaussQuad[T](f: NumContextProc[T, float]; xStart, xEnd: float; N = 100;
                  nPoints = 7; ctx: NumContext[T, float] = nil): T

Calculate the integral of f using Gaussian Quadrature. Has 20 different sets of weights, ranging from 1 to 20 function evaluations per subinterval.

Input:

  • f: the function that is integrated.
  • xStart: The start of the integration interval.
  • xEnd: The end of the integration interval.
  • N: The number of subintervals to divide the integration interval into.
  • nPoints: The number of points to evaluate f at per interval. Choose between 1 to 20 with increasing accuracy.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • The integral evaluated from the xStart to xEnd calculated using Gaussian Quadrature.
Source   Edit   
proc insert[T; U; V](intervalList: var IntervalList[T, U, V];
                     el: IntervalType[T, U, V]) {.inline.}
Source   Edit   
proc pop[T; U; V](intervalList: var IntervalList[T, U, V]): IntervalType[T, U, V] {.
    inline.}
Source   Edit   
proc romberg[T](f: NumContextProc[T, float]; xStart, xEnd: float; depth = 8;
                tol = 1e-8; ctx: NumContext[T, float] = nil): T

Calculate the integral of f using Romberg Integration.

Input:

  • f: the function that is integrated.
  • xStart: The start of the integration interval.
  • xEnd: The end of the integration interval.
  • depth: The maximum depth of the Richardson Extrapolation.
  • tol: The error tolerance that must be satisfied.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • The value of the integral of f from xStart to xEnd calculated using Romberg integration.
Source   Edit   
proc romberg[T](Y: openArray[T]; X: openArray[float]): T

Calculate the integral of f using Romberg Integration from a set of values.

Input:

  • Y: A seq of values of the integrand.
  • X: A seq with the corresponding x-values.

Returns:

  • The integral evaluated from the smallest to the largest value in X calculated using Romberg Integration.
Source   Edit   
proc simpson[T](f: NumContextProc[T, float]; xStart, xEnd: float; N = 500;
                ctx: NumContext[T, float] = nil): T

Calculate the integral of f using Simpson's rule.

Input:

  • f: the function that is integrated.
  • xStart: The start of the integration interval.
  • xEnd: The end of the integration interval.
  • N: The number of subintervals to divide the integration interval into. Must be 2 or greater.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • The value of the integral of f from xStart to xEnd calculated using Simpson's rule.
Source   Edit   
proc simpson[T](Y: openArray[T]; X: openArray[float]): T

Calculate the integral of f using Simpson's rule from a set of values.

Input:

  • Y: A seq of values of the integrand.
  • X: A seq with the corresponding x-values.

Returns:

  • The integral evaluated from the smallest to the largest value in X calculated using Simpson's rule.
Source   Edit   
proc trapz[T](f: NumContextProc[T, float]; xStart, xEnd: float; N = 500;
              ctx: NumContext[T, float] = nil): T

Calculate the integral of f using the trapezoidal rule.

Input:

  • f: the function that is integrated.
  • xStart: The start of the integration interval.
  • xEnd: The end of the integration interval.
  • N: The number of subintervals to divide the integration interval into.
  • ctx: A context variable that can be accessed and modified in f. It is a ref type so IT IS MUTABLE. It can be used to save extra information during the solving for example, or to pass in big Tensors.

Returns:

  • The value of the integral of f from xStart to xEnd calculated using the trapezoidal rule.
Source   Edit   
proc trapz[T](Y: openArray[T]; X: openArray[float]): T

Calculate the integral of f using the trapezoidal rule from a set values.

Input:

  • Y: A seq of values of the integrand.
  • X: A seq with the corresponding x-values.

Returns:

  • The integral evaluated from the smallest to the largest value in X calculated using the trapezoidal rule.
Source   Edit  

Made with Nim. Generated: 2024-03-29 09:01:52 UTC