|
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.
The easy integral
int 01 x log(1 + x) dx = (1 /4) .
is chosen to illustrate how these integration schemes may be used.
The computation is performed using precision D=50.
> 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-20
and, 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
One could, for example, want to approximate the integral
int - 11 (cos(x)sin(x) /(1 - x) 1 /3(1 + x) 1 /3) dx approx 0.128866348618482136010608730745 .
using N=10 integration points in precision D=30.
> 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-20
A 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
The double-exponential integration can also be used to approximate tough integrals such as
int 0∞ e - xcos(x) dx = int 01 (e^(1 - (1 /z))cos((1 /z) - 1) /z 2) dz = (1 /2) .
> 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
[Next][Prev] [Right] [Left] [Up] [Index] [Root]
|