Numerical Functions

This section contains some functions for numerical analysis, taken from Pari.

Subsections

• Summation of Infinite Series
• Integration

Summation of Infinite Series

There are three functions for evaluating infinite sums of real numbers. The sum should be specified as a map m from the integers to the free real field, such that m(n) is the n^( th) term of the sum. The summation begins at term i. The precision of the result will be the default precision of the free real field.

InfiniteSum(m, i) : Map, RngIntElt -> FldPrElt

An approximation to the infinite sum m(i) + m(i + 1) + m(i + 2) + ... . This function also works for maps to the free complex field.

PositiveSum(m, i) : Map, RngIntElt -> FldPrElt

An approximation to the infinite sum m(i) + m(i + 1) + m(i + 2) + ... . Designed for series in which every term is positive, it uses van Wijngaarden's trick for converting the series into an alternating one. Due to the stopping criterion, terms equal to 0 will create problems and should be removed.

AlternatingSum(m, i) : Map, RngIntElt -> FldPrElt

    Al: MonStgElt                       Default: "Villegas"

An approximation to the infinite sum m(i) + m(i + 1) + m(i + 2) + ... . Designed for series in which the terms alternate in sign. The optional argument Al can be used to specify the algorithm used. The possible values are "Villegas" (the default), and "EulerVanWijngaarden". Due to the stopping criterion, terms equal to 0 will create problems and should be removed.

Integration

A number of `Romberg-like' integration methods have been taken from Pari. The precision should not be made too large for this, and singularities are not allowed in the interval of integration (including its boundaries).

Interpolation(P, V, x) : [FldPrElt], [FldPrElt], FldPrElt -> FldPrElt, FldPrElt

Using Neville's algorithm, interpolate the value of x under a polynomial p such that p(P[i]) = V[i]. An estimate of the error is also returned.

Integral(m, a, b) : Map, FldPrElt, FldPRElt -> FldPrElt

    Al: MonStgElt                       Default: "Open"

Given a function f, (at least) continuous everywhere on the open interval (a, b) and defined on [a, b], for free real variables a and b, this function numerically integrates f from a to b. The function f must be defined as a map with as its domain the free real field and as its codomain either the free real or the free complex field. The precision should be not much higher than 28. For improper integrals, where one or both of the limits of integration are plus or minus infinity, the integrand must decrease sufficiently rapidly at infinity (this can often be achieved by integration by parts).

Using the optional argument Al it is possible to choose a special algorithm. With Al := "Infinite" a method tailored for use when either a or b is `infinite' is used; in that case -a or b could be chosen `large' (how large depends on how rapidly the function tends to 0 at +- Infinity), but it must be the case that a and b have the same sign (and it is suggested that a and b have magnitude at least 1). When neither -a nor b are very large, f is smooth and it can be evaluated everywhere on the closed interval, the method chosen with Al := "Simple" is the fastest. The case Al := "Open" is best suited for cases in which f is undefined (but continuous) at one of the end-points.

The default algorithm depends on the limits of integration. It always uses the "Open" or "Infinite" algorithm (or some combination of these) and so the function need not be defined at a and b.

RombergQuadrature(f, a, b: parameters) : Program, FldPrElt, FldPrElt -> FldPrElt

    Precision: FldPrElt                 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.

SimpsonQuadrature(f, a, b, n) : Program, FldPrElt, FldPrElt, RngIntElt -> FldPrElt

Using Simpson's rule on n sub-intervals, approximate the integral of f from a to b.

TrapezoidalQuadrature(f, a, b, n) : Program, FldPrElt, FldPrElt, RngIntElt -> FldPrElt

Using the trapezoidal rule on n sub-intervals, approximate the integral of f from a to b.

Example FldRe_Integral (H41E9)

> R := RealField();
> f := map< R -> R | x :-> Exp(-x^2) >;
> Integral(f, 0, 10);
0.88622692545275801364908374188
> f := map< R -> R | x :-> x^2 >;
> Integral(f, -1, 2);
2.99999999999999999999999999888

> f := map< R -> R | x :-> 1/(x+4)^2 >;
> Integral(f, 1, 10^4000: Al := "Infinite");
0.199999999999999999999999989688
> Integral(f, 1, 10^40: Al := "Infinite");
0.199999999999999999999999989688

> R := RealField();
> C<i> := ComplexField();
> f := map< R -> C | x :-> i*x^2 + 5*x>;
> Integral(f, 0, 10: Al := "Simple");
250.000000000000000000000000 + 333.333333333333333333333333335487*i

> f := map< R -> R | x :-> Sin(x)/x >;
> Integral(f, 0, 1: Al := "Open");
0.94608307036718301494135330891