Two groups of intrinsics for numerical integration are provided. The first group are designed for the numerical integration of complex-valued functions along the interval [ - 1, 1]. They have been researched and implemented by Christian Neurohr. The second group consist of three basic Romberg-like numerical integration methods that are based on those in Pari.
Given n + 1 points (xi, yi), there is a unique polynomial p(x) of degree n which passes through these points. Given some t, polynomial interpolation computes the value u = p(t).
Using Neville's algorithm, interpolate the value of t with respect to a polynomial p such that p(P[i]) = V[i]. An estimate of the error is also returned.
An intrinsic is provided to compute the discrete Fourier transform of a sequence of complex numbers.
Given a sequence E of complex numbers, return their discrete Fourier transform as a sequence of complex numbers.
In this section several integration schemes that are useful for numerical integration along the interval [ - 1, 1] are described. They have been designed and implemented by Christian Neurohr. In particular, the routines are suitable for numerically approximating integrals of complex-valued functions that are analytic over the interval [-1,1] or have integrable singularities at the end points. These should not be regarded as general purpose numerical integration tools. The examples in this section merely indicate how to use these integration methods. Error bounds for them and more details on how to use them can be found in [Neu18, Chapter 3] as well as many other references on numerical analysis.
Gaussian quadratures are well-known quadrature rules that integrate exactly polynomials of degree 2N - 1 using N integration points. In particular, they are highly efficient for integrands that are well-behaved around the integration path. The following functions compute the Gauss-Jacobi integration scheme on the interval [ - 1, 1], i.e. a sequence of abscissas A and a sequence of weights W such that can be used to approximate the definite integral int - 11 (1 - x)a(1 + x)b f(x) dx approx ∑k=1N W[k] .A[k] for a well-behaved, complex-valued function f : [ - 1, 1] -> C, a, b > - 1 and N sufficiently large. The Gauss-Legendre integration scheme occurs as the special case a=b=0. For more details see also [Neu18, Section 3.2].
Construct the Gauss--Legendre integration scheme on N points and precision D.
Construct the Gauss--Jacobi integration scheme on N points and precision D.
The Clenshaw--Curtis quadrature on N + 1 is an interpolation integration scheme that integrates polynomials exactly up to degree N. The abscissas are the N - 1 extrema of the N-th Chebyshev polynomial plus the boundary points ∓ 1 and the weights are obtained by interpolation. It is best suited for integrands that are analytic around [ - 1, 1] and can be computed quite quickly using a discrete Fourier transform. As before, a sequence of abscissas A and a sequence of weights W on interval [ - 1, 1] can be used to approximate int - 11 f(x) dx approx ∑k=1N W[k] .A[k] for a well-behaved, complex-valued function f : [ - 1, 1] -> C and N chosen large enough. For more details on Clenshaw--Curtis integration and this particular implementation see [Neu18, Section 3.3].
Construct the Clenshaw--Curtis integration scheme on N + 1 points and precision D.
Double--exponential integration (or tanh--sinh quadrature) is a quite versatile integration scheme. While it is far less efficient than Gaussian or Clenshaw--Curtis quadrature for well-behaved integrands, double-exponential integration can easily handle difficult integrands with endpoint singularities. In fact, the double--exponential integration as implemented here can handle integrands similar to Gauss--Jacobi quadrature. It returns three sequences of real numbers of length 2N + 1: abscissas A, weights W1 and 'extra weights' W2. int - 11 (1 - x)a(1 + x)b f(x) dx approx ∑k=12N + 1 W1[k] .W2[k]b .(1 + A[k])b - a .f(A[k]) . for a well--behaved, complex--valued function f : [ - 1, 1] -> C and N chosen large enough and, without loss of generality, b>a.
Construct the double--exponential integration scheme with 2N + 1 points and step-length h>0. The precision of the abscissas and weights is inherited from h.
> N := 20; > D := 50; > function f(x) > return x*Log(1+x); > end function; > A_GL, W_GL := GaussLegendreIntegrationPoints(N,D); > /* Map abscissas to the interval [0,1] */ > A_GL := [ (1/2)*(a+1) : a in A_GL ]; > int_GL := (1/2) * &+[ W_GL[k] * f(A_GL[k]) : k in [1..N] ]; > Abs(int_GL - 1/4); 9.3232511485451973751923945899073760403381175208927E-33
Now, using Clenshaw--Curtis integration on N + 1=21 points
> A_CC, W_CC := ClenshawCurtisIntegrationPoints(N,D); > A_CC := [ (1/2)*(a+1) : a in A_CC ]; > int_CC := (1/2) * &+[ W_CC[k] * f(A_CC[k]) : k in [1..N+1] ]; > Abs(int_CC - 1/4); 1.4644696752310341410357330113758869025438089669910E-20and, finally, double-exponential with steplength h=6/N
> h := RealField(D)!6/N; > A_DE, W_DE := TanhSinhIntegrationPoints(N,h); > A_DE := [ (1/2)*(a+1) : a in A_DE ]; > int_DE := (1/2) * &+[ W_DE[k] * f(A_DE[k]) : k in [1..2*N+1] ]; > Abs(int_DE - 1/4); 4.4223451294316356251002317691142668875095964280941E-10
> c := 0.128866348618482136010608730745; > N := 10; > D := 30; > function f(x) > return Sin(x)*Cos(x); > end function; > A, W := GaussJacobiIntegrationPoints(N,D,-1/3,-1/5); > int_GJ := &+[ W[k] * f(A[k]) : k in [1..N] ]; > Abs(int_GJ - c); 1.14406228730447597223836768506E-20A similar result can be achieved when using double-exponential integration (or Tanh--Sinh quadrature) with steplength h = 1/8 and N = 48
> h := Real(2^-3); > N := Ceiling(6/h); > A, W1, W2 := TanhSinhIntegrationPoints(N,h); > int_DE := &+[ W1[k]*W2[k]^(1/3)*(1+A[k])^(1/3-1/5)*f(A[k]) : k in [1..2*N+1] ]; > Abs(int_DE - c); 1.32999219914892690177636985694E-20
> h := RealField(1000)!(2^-10); > N := Ceiling(7.2/h); > A_DE, W_DE := TanhSinhIntegrationPoints(N,h); > assert A_DE[1] gt -1; > assert A_DE[2*N+1] lt 1; > A_DE := [ (1/2)*(a+1) : a in A_DE ]; > function f(z) > return (Exp(1-1/z) * Cos(1/z-1)/z^2); > end function; > int_DE := (1/2) * &+[ W_DE[k] * f(A_DE[k]) : k in [1..2*N+1] ]; > Abs(int_DE - 1/2) lt 10^-300; true
A number of `Romberg-Type' integration methods are provided which are based on Pari implementations. The precision should not be made too large for this, and the range of integration (including its boundaries) must be free of singularities.
Precision: FldReElt Default: 1.0e-6
MaxSteps: RngIntElt Default: 20
K: RngIntElt Default: 5
Using Romberg's method of order 2K, approximate the integral of f from a to b. The desired accuracy may be specified by setting the Precision parameter, and the order of the algorithm by changing K. The algorithm ceases after MaxSteps iterations if the desired accuracy has not been achieved.
Using Simpson's rule on n sub-intervals, approximate the integral of f from a to b.
Using the trapezoidal rule on n sub-intervals, approximate the integral of f from a to b.
A function is provided to compute the numerical derivative of a function. This is achieved by computing enough interpolation points and performing a Taylor expansion.
Given a suitably nice function f, compute a numerical approximation to the nth derivative at the point z.
> f := func<x|Exp(2*x)>; > NumericalDerivative(f, 10, ComplexField(30)! 1.0) / f (1.0); 1024.00000000000000000000000000 > NumericalDerivative(func<x|LogGamma(x)>,1,ComplexField()!3.0); 0.922784335098467139393487909918 > Psi(3.0); // Psi is Gamma'/Gamma 0.922784335098467139393487909918