Numerical integration/Gauss-Legendre Quadrature
You are encouraged to solve this task according to the task description, using any language you may know.
| In a general Gaussian quadrature rule, an definite integral of is first approximated over the interval by a polynomial approximable function and a known weighting function . | |
| Those are then approximated by a sum of function values at specified points multiplied by some weights : | |
| In the case of Gauss-Legendre quadrature, the weighting function , so we can approximate an integral of with: |
![{\displaystyle [-1,1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/51e3b7f14a6f70e614728c583409a0b9a8b9de01)
For this, we first need to calculate the nodes and the weights, but after we have them, we can reuse them for numerious integral evaluations, which greatly speeds up the calculation compared to more simple numerical integration methods.
| The evaluation points for a n-point rule, also called "nodes", are roots of n-th order Legendre Polynomials . Legendre polynomials are defined by the following recursive rule: |
|
| There is also a recursive equation for their derivative: | |
| The roots of those polynomials are in general not analytically solvable, so they have to be approximated numerically, for example by Newton-Raphson iteration: | |
| The first guess for the -th root of a -order polynomial can be given by | |
| After we get the nodes , we compute the appropriate weights by: | |
| After we have the nodes and the weights for a n-point quadrature rule, we can approximate an integral over any interval by |
![{\displaystyle w_i = \frac{2}{\left( 1-x_i^2 \right) [P'_n(x_i)]^2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1f1e447dea4bf0b3477dc54698ec8bbd50bb62c9)
![{\displaystyle [a,b]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9c4b788fc5c637e26ee98b45f89a5c08c85f7935)
Task description
Similar to the task Numerical Integration, the task here is to calculate the definite integral of a function , but by applying an n-point Gauss-Legendre quadrature rule, as described here, for example. The input values should be an function f to integrate, the bounds of the integration interval a and b, and the number of gaussian evaluation points n. An reference implementation in Common Lisp is provided for comparison.
To demonstrate the calculation, compute the weights and nodes for an 5-point quadrature rule and then use them to compute
C
<lang C>#include <stdio.h>
- include <math.h>
- define N 5
double Pi; double lroots[N]; double weight[N]; double lcoef[N + 1][N + 1] = Template:0;
void lege_coef() { int n, i; lcoef[0][0] = lcoef[1][1] = 1; for (n = 2; n <= N; n++) { lcoef[n][0] = -(n - 1) * lcoef[n - 2][0] / n; for (i = 1; i <= n; i++) lcoef[n][i] = ((2 * n - 1) * lcoef[n - 1][i - 1] - (n - 1) * lcoef[n - 2][i] ) / n; } }
double lege_eval(int n, double x) { int i; double s = lcoef[n][n]; for (i = n; i; i--) s = s * x + lcoef[n][i - 1]; return s; }
double lege_diff(int n, double x) { return n * (x * lege_eval(n, x) - lege_eval(n - 1, x)) / (x * x - 1); }
void lege_roots() { int i; double x, x1; for (i = 1; i <= N; i++) { x = cos(Pi * (i - .25) / (N + .5)); do { x1 = x; x -= lege_eval(N, x) / lege_diff(N, x); } while (x != x1); /* x != x1 is normally a no-no, but this task happens to be * well behaved. */ lroots[i - 1] = x;
x1 = lege_diff(N, x); weight[i - 1] = 2 / ((1 - x * x) * x1 * x1); } }
double lege_inte(double (*f)(double), double a, double b) { double c1 = (b - a) / 2, c2 = (b + a) / 2, sum = 0; int i; for (i = 0; i < N; i++) sum += weight[i] * f(c1 * lroots[i] + c2); return c1 * sum; }
int main() { int i; Pi = atan2(1, 1) * 4;
lege_coef(); lege_roots();
printf("Roots: "); for (i = 0; i < N; i++) printf(" %g", lroots[i]);
printf("\nWeight:"); for (i = 0; i < N; i++) printf(" %g", weight[i]);
printf("\nintegrating Exp(x) over [-3, 3]:\n\t%10.8f,\n" "compred to actual\n\t%10.8f\n", lege_inte(exp, -3, 3), exp(3) - exp(-3)); return 0; }</lang>
- Output:
Roots: 0.90618 0.538469 0 -0.538469 -0.90618
Weight: 0.236927 0.478629 0.568889 0.478629 0.236927
integrating Exp(x) over [-3, 3]:
20.03557772,
compred to actual
20.03574985
Common Lisp
<lang lisp>;; Computes the initial guess for the root i of a n-order Legendre polynomial. (defun guess (n i)
(cos (* pi
(/ (- i 0.25d0)
(+ n 0.5d0)))))
- Computes and evaluates the n-order Legendre polynomial at the point x.
(defun legpoly (n x)
(let ((pa 1.0d0)
(pb x)
(pn))
(cond ((= n 0) pa)
((= n 1) pb)
(t (loop for i from 2 to n do
(setf pn (- (* (/ (- (* 2 i) 1) i) x pb)
(* (/ (- i 1) i) pa)))
(setf pa pb)
(setf pb pn)
finally (return pn))))))
- Computes and evaluates the derivative of an n-order Legendre polynomial at point x.
(defun legdiff (n x)
(* (/ n (- (* x x) 1))
(- (* x (legpoly n x))
(legpoly (- n 1) x))))
- Computes the n nodes for an n-point quadrature rule. (i.e. n roots of a n-order polynomial)
(defun nodes (n)
(let ((x (make-array n :initial-element 0.0d0)))
(loop for i from 0 to (- n 1) do
(let ((val (guess n (+ i 1))) ;Nullstellen-Schätzwert.
(itermax 5))
(dotimes (j itermax)
(setf val (- val
(/ (legpoly n val)
(legdiff n val)))))
(setf (aref x i) val)))
x))
- Computes the weight for an n-order polynomial at the point (node) x.
(defun legwts (n x)
(/ 2
(* (- 1 (* x x))
(expt (legdiff n x) 2))))
- Takes a array of nodes x and computes an array of corresponding weights w.
(defun weights (x)
(let* ((n (car (array-dimensions x)))
(w (make-array n :initial-element 0.0d0)))
(loop for i from 0 to (- n 1) do
(setf (aref w i) (legwts n (aref x i))))
w))
- Integrates a function f with a n-point Gauss-Legendre quadrature rule over the interval [a,b].
(defun int (f n a b)
(let* ((x (nodes n))
(w (weights x)))
(* (/ (- b a) 2.0d0)
(loop for i from 0 to (- n 1)
sum (* (aref w i)
(funcall f (+ (* (/ (- b a) 2.0d0)
(aref x i))
(/ (+ a b) 2.0d0))))))))</lang>
- Example:
<lang lisp>(nodes 5)
- (0.906179845938664d0 0.5384693101056831d0 2.996272867003007d-95 -0.5384693101056831d0 -0.906179845938664d0)
(weights (nodes 5))
- (0.23692688505618917d0 0.47862867049936647d0 0.5688888888888889d0 0.47862867049936647d0 0.23692688505618917d0)
(int #'exp 5 -3 3) 20.035577718385568d0</lang> Comparison of the 5-point rule with simpler, but more costly methods from the task Numerical Integration: <lang lisp>(int #'(lambda (x) (expt x 3)) 5 0 1) 0.24999999999999997d0
(int #'(lambda (x) (/ 1 x)) 5 1 100) 4.059147508941519d0
(int #'(lambda (x) x) 5 0 5000) 1.25d7
(int #'(lambda (x) x) 5 0 6000) 1.8000000000000004d7</lang>
Delphi
<lang Delphi>program Legendre;
{$APPTYPE CONSOLE}
const Order = 5;
Epsilon = 1E-14;
var Roots : array[0..Order-1] of double;
Weight : array[0..Order-1] of double; LegCoef : array [0..Order,0..Order] of double;
function F(X:double) : double; begin
Result := Exp(X);
end;
procedure PrepCoef; var I, N : integer; begin
for I:=0 to Order do
for N := 0 to Order do
LegCoef[I,N] := 0;
LegCoef[0,0] := 1;
LegCoef[1,1] := 1;
For N:=2 to Order do
begin
LegCoef[N,0] := -(N-1) * LegCoef[N-2,0] / N;
For I := 1 to Order do
LegCoef[N,I] := ((2*N-1) * LegCoef[N-1,I-1] - (N-1)*LegCoef[N-2,I]) / N;
end;
end;
function LegEval(N:integer; X:double) : double; var I : integer; begin
Result := LegCoef[n][n]; for I := N-1 downto 0 do Result := Result * X + LegCoef[N][I];
end;
function LegDiff(N:integer; X:double) : double; begin
Result := N * (X * LegEval(N,X) - LegEval(N-1,X)) / (X*X-1);
end;
procedure LegRoots; var I : integer;
X, X1 : double;
begin
for I := 1 to Order do
begin
X := Cos(Pi * (I-0.25) / (Order+0.5));
repeat
X1 := X;
X := X - LegEval(Order,X) / LegDiff(Order, X);
until Abs (X-X1) < Epsilon;
Roots[I-1] := X;
X1 := LegDiff(Order,X);
Weight[I-1] := 2 / ((1-X*X) * X1*X1);
end;
end;
function LegInt(A,B:double) : double; var I : integer;
C1, C2 : double;
begin
C1 := (B-A)/2; C2 := (B+A)/2; Result := 0; For I := 0 to Order-1 do Result := Result + Weight[I] * F(C1*Roots[I] + C2); Result := C1 * Result;
end;
var I : integer;
begin
PrepCoef; LegRoots;
Write('Roots: ');
for I := 0 to Order-1 do
Write (' ',Roots[I]:13:10);
Writeln;
Write('Weight: ');
for I := 0 to Order-1 do
Write (' ', Weight[I]:13:10);
writeln;
Writeln('Integrating Exp(x) over [-3, 3]: ',LegInt(-3,3):13:10);
Writeln('Actual value: ',Exp(3)-Exp(-3):13:10);
Readln;
end.</lang>
Roots: 0.9061798459 0.5384693101 0.0000000000 -0.5384693101 -0.9061798459 Weight: 0.2369268851 0.4786286705 0.5688888889 0.4786286705 0.2369268851 Integrating Exp(X) over [-3, 3]: 20.0355777184 Actual value: 20.0357498548
D
<lang d>import std.stdio, std.math;
immutable struct GaussLegendreQuadrature(size_t N, FP=double,
size_t NBITS=50) {
immutable static double[N] lroots, weight;
alias FP[N + 1][N + 1] CoefMat;
pure nothrow static this() {
static FP legendreEval(in ref FP[N + 1][N + 1] lcoef,
in int n, in FP x) pure nothrow {
FP s = lcoef[n][n];
foreach_reverse (immutable i; 1 .. n+1)
s = s * x + lcoef[n][i - 1];
return s;
}
static FP legendreDiff(in ref CoefMat lcoef,
in int n, in FP x) pure nothrow {
return n * (x * legendreEval(lcoef, n, x) -
legendreEval(lcoef, n - 1, x)) /
(x ^^ 2 - 1);
}
static void legendreRoots(in ref CoefMat lcoef) pure nothrow {
foreach (immutable i; 1 .. N + 1) {
FP x = cos(PI * (i - 0.25) / (N + 0.5));
FP x1;
do {
x1 = x;
x -= legendreEval(lcoef, N, x) /
legendreDiff(lcoef, N, x);
} while (feqrel(x, x1) < NBITS);
lroots[i - 1] = x;
x1 = legendreDiff(lcoef, N, x);
weight[i - 1] = 2 / ((1 - x ^^ 2) * (x1 ^^ 2));
}
}
CoefMat lcoef = 0.0;
legendreCoefInit(/*ref*/ lcoef);
legendreRoots(lcoef);
}
static private void legendreCoefInit(ref CoefMat lcoef)
pure nothrow {
lcoef[0][0] = lcoef[1][1] = 1;
foreach (immutable int n; 2 .. N + 1) { // n must be signed.
lcoef[n][0] = -(n - 1) * lcoef[n - 2][0] / n;
foreach (immutable i; 1 .. n + 1)
lcoef[n][i] = ((2 * n - 1) * lcoef[n - 1][i - 1] -
(n - 1) * lcoef[n - 2][i]) / n;
}
}
static public FP integrate(in FP function(FP x) f,
in FP a, in FP b) {
immutable FP c1 = (b - a) / 2;
immutable FP c2 = (b + a) / 2;
FP sum = 0.0;
foreach (immutable i; 0 .. N)
sum += weight[i] * f(c1 * lroots[i] + c2);
return c1 * sum;
}
}
void main() {
GaussLegendreQuadrature!(5, real) glq;
writeln("Roots: ", glq.lroots);
writeln("Weight: ", glq.weight);
writefln("Integrating exp(x) over [-3, 3]: %10.12f",
glq.integrate(&exp, -3, 3));
writefln("Compred to actual: %10.12f",
3.0.exp - -3.0.exp);
}</lang>
- Output:
Roots: [0.90618, 0.538469, 0, -0.538469, -0.90618] Weight: [0.236927, 0.478629, 0.568889, 0.478629, 0.236927] Integrating exp(x) over [-3, 3]: 20.035577718386 Compred to actual: 20.035749854820
Fortran
<lang Fortran>! Works with gfortran but needs the option ! -assume realloc_lhs ! when compiled with Intel Fortran.
program gauss
implicit none integer, parameter :: p = 16 ! quadruple precision integer :: n = 10, k real(kind=p), allocatable :: r(:,:) real(kind=p) :: z, a, b, exact do n = 1,20 a = -3; b = 3 r = gaussquad(n) z = (b-a)/2*dot_product(r(2,:),exp((a+b)/2+r(1,:)*(b-a)/2)) exact = exp(3.0_p)-exp(-3.0_p) print "(i0,1x,g0,1x,g10.2)",n, z, z-exact end do contains
function gaussquad(n) result(r)
integer :: n
real(kind=p), parameter :: pi = 4*atan(1._p)
real(kind=p) :: r(2, n), x, f, df, dx
integer :: i, iter
real(kind = p), allocatable :: p0(:), p1(:), tmp(:)
p0 = [1._p]
p1 = [1._p, 0._p]
do k = 2, n
tmp = ((2*k-1)*[p1,0._p]-(k-1)*[0._p, 0._p,p0])/k
p0 = p1; p1 = tmp
end do
do i = 1, n
x = cos(pi*(i-0.25_p)/(n+0.5_p))
do iter = 1, 10
f = p1(1); df = 0._p
do k = 2, size(p1)
df = f + x*df
f = p1(k) + x * f
end do
dx = f / df
x = x - dx
if (abs(dx)<10*epsilon(dx)) exit
end do
r(1,i) = x
r(2,i) = 2/((1-x**2)*df**2)
end do
end function
end program </lang>
n numerical integral error -------------------------------------------------- 1 6.00000000000000000000000000000000 -14. 2 17.4874646410555689643606840462449 -2.5 3 19.8536919968055821921309108927158 -.18 4 20.0286883952907008527738054439858 -.71E-02 5 20.0355777183855621539285357252751 -.17E-03 6 20.0357469750923438830654575585499 -.29E-05 7 20.0357498197266007755718729372892 -.35E-07 8 20.0357498544945172882260918041684 -.33E-09 9 20.0357498548174338368864419454859 -.24E-11 10 20.0357498548197898711175766908548 -.14E-13 11 20.0357498548198037305529147159695 -.67E-16 12 20.0357498548198037976759531014464 -.27E-18 13 20.0357498548198037979482458119095 -.94E-21 14 20.0357498548198037979491844483597 -.28E-23 15 20.0357498548198037979491872317190 -.72E-26 16 20.0357498548198037979491872388913 -.40E-28 17 20.0357498548198037979491872389166 -.15E-28 18 20.0357498548198037979491872389259 -.58E-29 19 20.0357498548198037979491872388910 -.41E-28 20 20.0357498548198037979491872388495 -.82E-28
Go
Implementation pretty much by the methods given in the task description. <lang go>package main
import (
"fmt" "math"
)
// cFunc for continuous function. A type definition for convenience. type cFunc func(float64) float64
func main() {
fmt.Println("integral:", glq(math.Exp, -3, 3, 5))
}
// glq integrates f from a to b by Guass-Legendre quadrature using n nodes. // For the task, it also shows the intermediate values determining the nodes: // the n roots of the order n Legendre polynomal and the corresponding n // weights used for the integration. func glq(f cFunc, a, b float64, n int) float64 {
x, w := glqNodes(n, f)
show := func(label string, vs []float64) {
fmt.Printf("%8s: ", label)
for _, v := range vs {
fmt.Printf("%8.5f ", v)
}
fmt.Println()
}
show("nodes", x)
show("weights", w)
var sum float64
bma2 := (b - a) * .5
bpa2 := (b + a) * .5
for i, xi := range x {
sum += w[i] * f(bma2*xi+bpa2)
}
return bma2 * sum
}
// glqNodes computes both nodes and weights for a Gauss-Legendre // Quadrature integration. Parameters are n, the number of nodes // to compute and f, a continuous function to integrate. Return // values have len n. func glqNodes(n int, f cFunc) (node []float64, weight []float64) {
p := legendrePoly(n)
pn := p[n]
n64 := float64(n)
dn := func(x float64) float64 {
return (x*pn(x) - p[n-1](x)) * n64 / (x*x - 1)
}
node = make([]float64, n)
for i := range node {
x0 := math.Cos(math.Pi * (float64(i+1) - .25) / (n64 + .5))
node[i] = newtonRaphson(pn, dn, x0)
}
weight = make([]float64, n)
for i, x := range node {
dnx := dn(x)
weight[i] = 2 / ((1 - x*x) * dnx * dnx)
}
return
}
// legendrePoly constructs functions that implement Lengendre polynomials. // This is done by function composition by recurrence relation (Bonnet's.) // For given n, n+1 functions are returned, computing P0 through Pn. func legendrePoly(n int) []cFunc {
r := make([]cFunc, n+1)
r[0] = func(float64) float64 { return 1 }
r[1] = func(x float64) float64 { return x }
for i := 2; i <= n; i++ {
i2m1 := float64(i*2 - 1)
im1 := float64(i - 1)
rm1 := r[i-1]
rm2 := r[i-2]
invi := 1 / float64(i)
r[i] = func(x float64) float64 {
return (i2m1*x*rm1(x) - im1*rm2(x)) * invi
}
}
return r
}
// newtonRaphson is general purpose, although totally primitive, simply // panicking after a fixed number of iterations without convergence to // a fixed error. Parameter f must be a continuous function, // df its derivative, x0 an initial guess. func newtonRaphson(f, df cFunc, x0 float64) float64 {
for i := 0; i < 30; i++ {
x1 := x0 - f(x0)/df(x0)
if math.Abs(x1-x0) <= math.Abs(x0*1e-15) {
return x1
}
x0 = x1
}
panic("no convergence")
}</lang>
- Output:
nodes: 0.90618 0.53847 0.00000 -0.53847 -0.90618 weights: 0.23693 0.47863 0.56889 0.47863 0.23693 integral: 20.035577718385564
J
Solution: <lang j>P =: 3 :0 NB. list of coefficients for yth Legendre polynomial
if. y<:1 do. 1{.~->:y return. end.
y%~ (<:(,~+:)y) -/@:* (0,P<:y),:(P y-2)
)
getpoints =: 3 :0 NB. points,:weights for y points
x=. 1{:: p. p=.P y
w=. 2% (-.*:x)**:(p..p)p.x
x,:w
)
GaussLegendre =: 1 :0 NB. npoints function GaussLegendre (a,b)
'x w'=.getpoints x -:(-~/y)* +/w* u -:((+/,-~/)y)p.x
)</lang>
- Example use:
<lang j> 5 ^ GaussLegendre _3 3 20.0356</lang>
Julia
Julia provides a built-in function quadgk for adaptive Gauss–Kronrod quadrature, which is built on a function gauss(T, n) to compute the points and weights for an n-point Gauss–Legendre quadrature rule on (-1,1), which can be easily rescaled for any interval
<lang julia>x,w = Base.gauss(Float64, 5)
sum(exp(3x) .* 3w)</lang>
- Output:
20.035577718385564
Mathematica
code assumes function to be integrated has attribute Listable which is true of most built in Mathematica functions <lang Mathematica>gaussLegendreQuadrature[func_, {a_, b_}, degree_: 5] := Block[{nodes, x, weights},
nodes = Cases[NSolve[LegendreP[degree, x] == 0, x], _?NumericQ, Infinity]; weights = 2 (1 - nodes^2)/(degree LegendreP[degree - 1, nodes])^2; (b - a)/2 weights.func[(b - a)/2 nodes + (b + a)/2]]
gaussLegendreQuadrature[Exp, {-3, 3}]</lang>
- Output:
20.0356
Maxima
<lang maxima>gauss_coeff(n) := block([p, q, v, w],
p: expand(legendre_p(n, x)), q: expand(n/2*diff(p, x)*legendre_p(n - 1, x)), v: map(rhs, bfallroots(p)), w: map(lambda([z], 1/subst([x = z], q)), v), [map(bfloat, v), map(bfloat, w)])$
gauss_int(f, a, b, n) := block([u, x, w, c, h],
u: gauss_coeff(n), x: u[1], w: u[2], c: bfloat((a + b)/2), h: bfloat((b - a)/2), h*sum(w[i]*bfloat(f(c + x[i]*h)), i, 1, n))$
fpprec: 40$
gauss_int(lambda([x], 4/(1 + x^2)), 0, 1, 20);
/* 3.141592653589793238462643379852215927697b0 */
% - bfloat(%pi); /* -3.427286956499858315999116083264403489053b-27 */
gauss_int(exp, -3, 3, 5);
/* 2.003557771838556215392853572527509393154b1 */
% - bfloat(integrate(exp(x), x, -3, 3)); /* -1.721364342416440206515136565621888185351b-4 */</lang>
OCaml
<lang OCaml>let rec leg n x = match n with (* Evaluate Legendre polynomial *)
| 0 -> 1.0
| 1 -> x
| k -> let u = 1.0 -. 1.0 /. float k in
(1.0+.u)*.x*.(leg (k-1) x) -. u*.(leg (k-2) x);;
let leg' n x = match n with (* derivative *)
| 0 -> 0.0 | 1 -> 1.0 | _ -> ((leg (n-1) x) -. x*.(leg n x)) *. (float n)/.(1.0-.x*.x);;
let approx_root k n = (* Reversed Francesco Tricomi: 1 <= k <= n *)
let pi = acos (-1.0) and s = float(2*n) and t = 1.0 +. float(1-4*k)/.float(4*n+2) in (1.0 -. (float (n-1))/.(s*.s*.s))*.cos(pi*.t);;
let rec refine r n = (* Newton-Raphson *)
let r1 = r -. (leg n r)/.(leg' n r) in if abs_float (r-.r1) < 2e-16 then r1 else refine r1 n;;
let root k n = refine (approx_root k n) n;;
let node k n = (* Abscissa and weight *)
let r = root k n in let deriv = leg' n r in let w = 2.0/.((1.0-.r*.r)*.(deriv*.deriv)) in (r,w);;
let nodes n =
let rec aux k = if k > n then [] else node k n :: aux (k+1) in aux 1;;
let quadrature n f a b =
let f1 x = f ((x*.(b-.a) +. a +. b)*.0.5) in let eval s (x,w) = s +. w*.(f1 x) in 0.5*.(b-.a)*.(List.fold_left eval 0.0 (nodes n));;</lang>
which can be used in: <lang OCaml>let calc n =
Printf.printf
"Gauss-Legendre %2d-point quadrature for exp over [-3..3] = %.16f\n"
n (quadrature n exp (-3.0) 3.0);;
calc 5;; calc 10;; calc 15;; calc 20;;</lang>
- Output:
Gauss-Legendre 5-point quadrature for exp over [-3..3] = 20.0355777183855608 Gauss-Legendre 10-point quadrature for exp over [-3..3] = 20.0357498548197839 Gauss-Legendre 15-point quadrature for exp over [-3..3] = 20.0357498548198052 Gauss-Legendre 20-point quadrature for exp over [-3..3] = 20.0357498548198052
This shows convergence to the correct double-precision value of the integral <lang Ocaml>Printf.printf "%.16f\n" ((exp 3.0) -.(exp (-3.0)));; 20.0357498548198052</lang> although going beyond 20 points starts reducing the accuracy, due to accumulated rounding errors.
PARI/GP
This task is easy in GP thanks to built-in support for Legendre polynomials and efficient (Schonhage-Gourdon) polynomial root finding. <lang parigp>GLq(f,a,b,n)={
my(P=pollegendre(n),Pp=P',x=polroots(P)); (b-a)*sum(i=1,n,f((b-a)*x[i]/2+(a+b)/2)/(1-x[i]^2)/subst(Pp,'x,x[i])^2)
};
- \\ Turn on timer
GLq(x->exp(x), -3, 3, 5) \\ As of version 2.4.4, this can be written GLq(exp, -3, 3, 5)</lang>
- Output:
time = 0 ms. %1 = 20.035577718385562153928535725275093932 + 0.E-37*I
ooRexx
<lang oorexx>/*---------------------------------------------------------------------
- 31.10.2013 Walter Pachl Translation from REXX (from PL/I)
- using ooRexx' rxmath package
- which limits the precision to 16 digits
- --------------------------------------------------------------------*/
prec=60 Numeric Digits prec epsilon=1/10**prec pi=3.141592653589793238462643383279502884197169399375105820974944592307 exact = RxCalcExp(3,prec)-RxCalcExp(-3,prec) Do n = 1 To 20
a = -3; b = 3 r.=0 call gaussquad sum=0 Do j=1 To n sum=sum + r.2.j * RxCalcExp((a+b)/2+r.1.j*(b-a)/2,prec) End z = (b-a)/2 * sum Say right(n,2) format(z,2,40) format(z-exact,2,4,,0) End Say ' ' exact '(exact)' Exit
gaussquad:
p0.0=1; p0.1=1
p1.0=2; p1.1=1; p1.2=0
Do k = 2 To n
tmp.0=p1.0+1
Do L = 1 To p1.0
tmp.l = p1.l
End
tmp.l=0
tmp2.0=p0.0+2
tmp2.1=0
tmp2.2=0
Do L = 1 To p0.0
l2=l+2
tmp2.l2=p0.l
End
Do j=1 To tmp.0
tmp.j = ((2*k-1)*tmp.j - (k-1)*tmp2.j)/k
End
p0.0=p1.0
Do j=1 To p0.0
p0.j = p1.j
End
p1.0=tmp.0
Do j=1 To p1.0
p1.j=tmp.j
End
End
Do i = 1 To n
x = RxCalcCos(pi*(i-0.25)/(n+0.5),prec,'R')
Do iter = 1 To 10
f = p1.1; df = 0
Do k = 2 To p1.0
df = f + x*df
f = p1.k + x * f
End
dx = f / df
x = x - dx
If abs(dx) < epsilon Then Leave
End
r.1.i = x
r.2.i = 2/((1-x**2)*df**2)
End
Return
- requires 'rxmath' LIBRARY</lang>
Output:
1 6.0000000000000000000000000000000000000000 -1.4036E+1 2 17.4874646410555686000000000000000000000000 -2.5483 3 19.8536919968055914500000000000000000000000 -1.8206E-1 4 20.0286883952907032246391703165575495371776 -7.0615E-3 5 20.0355777183855623345965085871972344078167 -1.7214E-4 6 20.0357469750923433031000982816859525440756 -2.8797E-6 7 20.0357498197266007450081506439422093510041 -3.5093E-8 8 20.0357498544945192648654062025059252571210 -3.2529E-10 9 20.0357498548174362426073138353882519240177 -2.3698E-12 10 20.0357498548197905075149387536361754813374 -1.5552E-14 11 20.0357498548198049052166074059523608613749 -1.1548E-15 12 20.0357498548198068119347633275378821700762 7.5193E-16 13 20.0357498548198063256375618073806663013152 2.6564E-16 14 20.0357498548198035202546245888922276792447 -2.5397E-15 15 20.0357498548198027919824444452012138941729 -3.2680E-15 16 20.0357498548198037471314715729442546019171 -2.3129E-15 17 20.0357498548198067452377635761033686644343 6.8524E-16 18 20.0357498548198042026084719530842757694873 -1.8574E-15 19 20.0357498548198042304714191024916472961732 -1.8295E-15 20 20.0357498548198034525095801113268011014944 -2.6075E-15 20.03574985481980606 (exact)
Pascal
See Delphi
Perl 6
A free translation of the OCaml solution. We save half the effort to calculate the nodes by exploiting the (skew-)symmetry of the Legendre Polynomials. The evaluation of Pn(x) is kept linear in n by also passing Pn-1(x) in the recursion.
The quadrature function allows passing in a precalculated list of nodes for repeated integrations.
Note: The calculations of Pn(x) and P'n(x) could be combined to further reduce duplicated effort. We also could cache P'n(x) from the last Newton-Raphson step for the weight calculation.
<lang perl6>multi legendre-pair( 1 , $x) { $x, 1 } multi legendre-pair(Int $n, $x) {
my ($m1, $m2) = legendre-pair($n - 1, $x); my \u = 1 - 1 / $n; (1 + u) * $x * $m1 - u * $m2, $m1;
}
multi legendre( 0 , $ ) { 1 } multi legendre(Int $n, $x) { legendre-pair($n, $x)[0] }
multi legendre-prime( 0 , $ ) { 0 } multi legendre-prime( 1 , $ ) { 1 } multi legendre-prime(Int $n, $x) {
my ($m0, $m1) = legendre-pair($n, $x); ($m1 - $x * $m0) * $n / (1 - $x**2);
}
sub approximate-legendre-root(Int $n, Int $k) {
# Approximation due to Francesco Tricomi my \t = (4*$k - 1) / (4*$n + 2); (1 - ($n - 1) / (8 * $n**3)) * cos(pi * t);
}
sub newton-raphson(&f, &f-prime, $r is copy, :$eps = 2e-16) {
while abs(my \dr = - f($r) / f-prime($r)) >= $eps {
$r += dr;
}
$r;
}
sub legendre-root(Int $n, Int $k) {
newton-raphson(&legendre.assuming($n), &legendre-prime.assuming($n),
approximate-legendre-root($n, $k));
}
sub weight(Int $n, $r) { 2 / ((1 - $r**2) * legendre-prime($n, $r)**2) }
sub nodes(Int $n) {
gather {
take 0 => weight($n, 0) if $n !%% 2;
for 1 .. $n div 2 {
my $r = legendre-root($n, $_);
my $w = weight($n, $r);
take $r => $w, -$r => $w;
}
}
}
sub quadrature(Int $n, &f, $a, $b, :@nodes = nodes($n)) {
sub scale($x) { ($x * ($b - $a) + $a + $b) / 2 }
($b - $a) / 2 * [+] @nodes.map: { .value * f(scale(.key)) }
}
say "Gauss-Legendre $_.fmt('%2d')-point quadrature ∫₋₃⁺³ exp(x) dx ≈ ",
quadrature($_, &exp, -3, +3) for 5 .. 10, 20;</lang>
- Output:
Gauss-Legendre 5-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0355777183856 Gauss-Legendre 6-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0357469750923 Gauss-Legendre 7-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0357498197266 Gauss-Legendre 8-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0357498544945 Gauss-Legendre 9-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0357498548174 Gauss-Legendre 10-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0357498548198 Gauss-Legendre 20-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0357498548198
PL/I
Translated from Fortran. <lang PL/I>(subscriptrange, size, fofl): Integration_Gauss: procedure options (main);
declare (n, k) fixed binary; declare r(*,*) float (18) controlled; declare (z, a, b, exact) float (18);
do n = 1 to 20; a = -3; b = 3; if allocation(r) > 0 then free r; allocate r(2, n); r = 0; call gaussquad(n, r); z = (b-a)/2 * sum(r(2,*) * exp((a+b)/2+r(1,*)*(b-a)/2)); exact = exp(3.0q0)-exp(-3.0q0); put skip edit (n, z, z-exact) (f(5), f(25,16), e(15,2)); end;
gaussquad: procedure(n, r); /*declare n fixed binary, r(2, n) float (18);*/
declare n fixed binary, r(2, *) float (18);/* corrected */
declare pi float (18) value (4*atan(1.0q0));
declare (x, f, df, dx) float (18);
declare (i, iter, L) fixed binary;
declare (p0(*), p1(*), tmp(*), tmp2(*)) float (18) controlled;
allocate p0(1) initial (1);
allocate p1(2) initial (1, 0);
do k = 2 to n;
allocate tmp(hbound(p1)+1); do L = 1 to hbound(p1); tmp(L) = p1(L); end; tmp(L) = 0;
allocate tmp2(hbound(p0)+2); tmp2(1), tmp2(2) = 0;
do L = 1 to hbound(p0); tmp2(L+2) = p0(L); end;
tmp = ((2*k-1)*tmp - (k-1)*tmp2)/k;
free p0; allocate p0(hbound(p1)); p0 = p1;
free p1; allocate p1(hbound(tmp)); p1 = tmp;
free tmp, tmp2;
end;
do i = 1 to n;
x = cos(pi*(i-0.25q0)/(n+0.5q0));
do iter = 1 to 10;
f = p1(1); df = 0;
do k = 2 to hbound(p1);
df = f + x*df;
f = p1(k) + x * f;
end;
dx = f / df;
x = x - dx;
if abs(dx) < 10*epsilon(dx) then leave;
end;
r(1,i) = x;
r(2,i) = 2/((1-x**2)*df**2);
end;
end gaussquad;
end Integration_Gauss; </lang>
1 6.0000000000000000 -1.40E+0001
2 17.4874646410555690 -2.55E+0000
3 19.8536919968055822 -1.82E-0001
4 20.0286883952907009 -7.06E-0003
5 20.0355777183855621 -1.72E-0004
6 20.0357469750923439 -2.88E-0006
7 20.0357498197266008 -3.51E-0008
8 20.0357498544945173 -3.25E-0010
9 20.0357498548174338 -2.37E-0012
10 20.0357498548197897 -1.41E-0014
11 20.0357498548198037 -6.94E-0017
12 20.0357498548198037 -6.25E-0017
13 20.0357498548198037 -1.25E-0016
14 20.0357498548198026 -1.16E-0015
15 20.0357498548198144 1.06E-0014
16 20.0357498548198021 -1.74E-0015
17 20.0357498548198359 3.21E-0014
18 20.0357498548198473 4.35E-0014
19 20.0357498548198848 8.10E-0014
20 20.0357498548200728 2.69E-0013
program gave me an error message: D:\ig.pli(19:2) : IBM1937I S Extents for parameters must be asterisks or restricted expressions with computational type. I tried to correct that. ok?
Python
<lang Python>from numpy import *
- Recursive generation of the Legendre polynomial of order n
def Legendre(n,x): x=array(x) if (n==0): return x*0+1.0 elif (n==1): return x else: return ((2.0*n-1.0)*x*Legendre(n-1,x)-(n-1)*Legendre(n-2,x))/n
- Derivative of the Legendre polynomials
def DLegendre(n,x): x=array(x) if (n==0): return x*0 elif (n==1): return x*0+1.0 else: return (n/(x**2-1.0))*(x*Legendre(n,x)-Legendre(n-1,x))
- Roots of the polynomial obtained using Newton-Raphson method
def LegendreRoots(polyorder,tolerance=1e-20): if polyorder<2: err=1 # bad polyorder no roots can be found else: roots=[] # The polynomials are alternately even and odd functions. So we evaluate only half the number of roots. for i in range(1,int(polyorder)/2 +1): x=cos(pi*(i-0.25)/(polyorder+0.5)) error=10*tolerance iters=0 while (error>tolerance) and (iters<1000): dx=-Legendre(polyorder,x)/DLegendre(polyorder,x) x=x+dx iters=iters+1 error=abs(dx) roots.append(x) # Use symmetry to get the other roots roots=array(roots) if polyorder%2==0: roots=concatenate( (-1.0*roots, roots[::-1]) ) else: roots=concatenate( (-1.0*roots, [0.0], roots[::-1]) ) err=0 # successfully determined roots return [roots, err]
- Weight coefficients
def GaussLegendreWeights(polyorder): W=[] [xis,err]=LegendreRoots(polyorder) if err==0: W=2.0/( (1.0-xis**2)*(DLegendre(polyorder,xis)**2) ) err=0 else: err=1 # could not determine roots - so no weights return [W, xis, err]
- The integral value
- func : the integrand
- a, b : lower and upper limits of the integral
- polyorder : order of the Legendre polynomial to be used
def GaussLegendreQuadrature(func, polyorder, a, b): [Ws,xs, err]= GaussLegendreWeights(polyorder) if err==0: ans=(b-a)*0.5*sum( Ws*func( (b-a)*0.5*xs+ (b+a)*0.5 ) ) else: # (in case of error) err=1 ans=None return [ans,err]
- The integrand - change as required
def func(x): return exp(x)
order=5 [Ws,xs,err]=GaussLegendreWeights(order) if err==0: print "Order : ", order print "Roots : ", xs print "Weights : ", Ws else: print "Roots/Weights evaluation failed"
- Integrating the function
[ans,err]=GaussLegendreQuadrature(func , order, -3,3) if err==0: print "Integral : ", ans else: print "Integral evaluation failed"</lang>
- Output:
Order : 5 Roots : [-0.90617985 -0.53846931 0. 0.53846931 0.90617985] Weights : [ 0.23692689 0.47862867 0.56888889 0.47862867 0.23692689] Integral : 20.0355777184
Racket
Computation of the Legendre polynomials and derivatives:
<lang racket> (define (LegendreP n x)
(let compute ([n n] [Pn-1 x] [Pn-2 1])
(case n
[(0) Pn-2]
[(1) Pn-1]
[else (compute (- n 1)
(/ (- (* (- (* 2 n) 1) x Pn-1)
(* (- n 1) Pn-2)) n)
Pn-1)])))
(define (LegendreP′ n x)
(* (/ n (- (* x x) 1))
(- (* x (LegendreP n x))
(LegendreP (- n 1) x))))
</lang>
Computation of the Legendre polynomial roots:
<lang racket> (define (LegendreP-root n i)
; newton-raphson step
(define (newton-step x)
(- x (/ (LegendreP n x) (LegendreP′ n x))))
; initial guess
(define x0 (cos (* pi (/ (- i 1/4) (+ n 1/2)))))
; computation of a root with relative accuracy 1e-15
(if (< (abs x0) 1e-15)
0
(let next ([x′ (newton-step x0)] [x x0])
(if (< (abs (/ (- x′ x) (+ x′ x))) 1e-15)
x′
(next (newton-step x′) x′)))))
</lang>
Computation of Gauss-Legendre nodes and weights
<lang racket> (define (Gauss-Legendre-quadrature n)
;; positive roots
(define roots
(for/list ([i (in-range (floor (/ n 2)))])
(LegendreP-root n (+ i 1))))
;; weights for positive roots
(define weights
(for/list ([x (in-list roots)])
(/ 2 (- 1 (sqr x)) (sqr (LegendreP′ n x)))))
;; all roots and weights
(values (append (map - roots)
(if (odd? n) (list 0) '())
(reverse roots))
(append weights
(if (odd? n) (list (/ 2 (sqr (LegendreP′ n 0)))) '())
(reverse weights))))
</lang>
Integration using Gauss-Legendre quadratures:
<lang racket> (define (integrate f a b #:nodes (n 5))
(define m (/ (+ a b) 2)) (define d (/ (- b a) 2)) (define-values [x w] (Gauss-Legendre-quadrature n)) (define (g x) (f (+ m (* d x)))) (* d (+ (apply + (map * w (map g x))))))
</lang>
Usage:
<lang racket>< > (Gauss-Legendre-quadrature 5) '(-0.906179845938664 -0.5384693101056831 0 0.5384693101056831 0.906179845938664) '(0.23692688505618875 0.47862867049936625 128/225 0.47862867049936625 0.23692688505618875)
> (integrate exp -3 3) 20.035577718385547
> (- (exp 3) (exp -3) 20.035749854819805 </lang>
Accuracy of the method:
<lang racket> > (require plot) > (parameterize ([plot-x-label "Number of Gaussian nodes"]
[plot-y-label "Integration error"]
[plot-y-transform log-transform]
[plot-y-ticks (log-ticks #:base 10)])
(plot (points (for/list ([n (in-range 2 11)])
(list n (abs (- (integrate exp -3 3 #:nodes n)
(- (exp 3) (exp -3)))))))))
</lang>
REXX
version 1
<lang rexx>/*---------------------------------------------------------------------
- 31.10.2013 Walter Pachl Translation from PL/I
- 01.11.2014 -"- see Version 2 for improvements
- --------------------------------------------------------------------*/
Call time 'R' prec=60 Numeric Digits prec epsilon=1/10**prec pi=3.141592653589793238462643383279502884197169399375105820974944592307 exact = exp(3,prec)-exp(-3,prec) Do n = 1 To 20
a = -3; b = 3
r.=0
call gaussquad
sum=0
Do j=1 To n
sum=sum + r.2.j * exp((a+b)/2+r.1.j*(b-a)/2,prec)
End
z = (b-a)/2 * sum
Say right(n,2) format(z,2,40) format(z-exact,2,4,,0)
End
Say ' ' exact '(exact)'
say '... and took' format(time('E'),,2) "seconds"
Exit
gaussquad:
p0.0=1; p0.1=1
p1.0=2; p1.1=1; p1.2=0
Do k = 2 To n
tmp.0=p1.0+1
Do L = 1 To p1.0
tmp.l = p1.l
End
tmp.l=0
tmp2.0=p0.0+2
tmp2.1=0
tmp2.2=0
Do L = 1 To p0.0
l2=l+2
tmp2.l2=p0.l
End
Do j=1 To tmp.0
tmp.j = ((2*k-1)*tmp.j - (k-1)*tmp2.j)/k
End
p0.0=p1.0
Do j=1 To p0.0
p0.j = p1.j
End
p1.0=tmp.0
Do j=1 To p1.0
p1.j=tmp.j
End
End
Do i = 1 To n
x = cos(pi*(i-0.25)/(n+0.5),prec)
Do iter = 1 To 10
f = p1.1; df = 0
Do k = 2 To p1.0
df = f + x*df
f = p1.k + x * f
End
dx = f / df
x = x - dx
If abs(dx) < epsilon then leave
End
r.1.i = x
r.2.i = 2/((1-x**2)*df**2)
End
Return
cos: Procedure /* REXX ****************************************************************
- Return cos(x) -- with specified precision
- cos(x) = 1-(x**2/2!)+(x**4/4!)-(x**6/6!)+-...
- 920903 Walter Pachl
- /
Parse Arg x,prec If prec= Then prec=9 Numeric Digits (2*prec) Numeric Fuzz 3 o=1 u=1 r=1 Do i=1 By 2 ra=r o=-o*x*x u=u*i*(i+1) r=r+(o/u) If r=ra Then Leave End Numeric Digits prec Return r+0
exp: Procedure /***********************************************************************
- Return exp(x) -- with reasonable precision
- 920903 Walter Pachl
- /
Parse Arg x,prec If prec<9 Then prec=9 Numeric Digits (2*prec) Numeric Fuzz 3 o=1 u=1 r=1 Do i=1 By 1 ra=r o=o*x u=u*i r=r+(o/u) If r=ra Then Leave End Numeric Digits (prec) Return r+0</lang>
Output:
1 6.0000000000000000000000000000000000000000 -1.4036E+1 2 17.4874646410555689643606840462449458421154 -2.5483 3 19.8536919968055821921309108927158495960775 -1.8206E-1 4 20.0286883952907008527738054439857661647073 -7.0615E-3 5 20.0355777183855621539285357252750939315016 -1.7214E-4 6 20.0357469750923438830654575585499253741530 -2.8797E-6 7 20.0357498197266007755718729372891903369401 -3.5093E-8 8 20.0357498544945172882260918041683132616237 -3.2529E-10 9 20.0357498548174338368864419454858704839263 -2.3700E-12 10 20.0357498548197898711175766908543458234008 -1.3927E-14 11 20.0357498548198037305529147159697031241994 -6.7396E-17 12 20.0357498548198037976759531014454017742327 -2.7323E-19 13 20.0357498548198037979482458119092690701863 -9.4143E-22 14 20.0357498548198037979491844483599375945130 -2.7906E-24 15 20.0357498548198037979491872317401917248453 -7.1915E-27 16 20.0357498548198037979491872389153958789316 -1.6260E-29 17 20.0357498548198037979491872389316236038179 -3.2517E-32 18 20.0357498548198037979491872389316560624361 -5.7920E-35 19 20.0357498548198037979491872389316561202637 -9.2480E-38 20 20.0357498548198037979491872389316561203561 -1.3311E-40 20.0357498548198037979491872389316561203562082463657269288113 (exact) ... and took 4.97 seconds
version 2
This REXX version is an optimized version (of version 1) which uses:
- a faster cos subroutine (which uses radian normalization)
- a faster exp subroutine
- some simple variables instead of stemmed arrays
- some static variables instead of repeated expressions
- multiplication [...*.5] instead of division [.../2]
- a generic approach for setting the numeric DIGITS
- an increase (by +1) of the number of iterations
Note that the function values for pi and e should have more precision than the number of digits specified.
The use of vertical bars is one of the very few times to use leading comments, as there isn't that many situations where there
exists nested DO loops with different (grouped) indentations, and practically no space on the right side of the statements.
It presents a good visual indication of what's what, but it's the dickens to pay when updating the code.
<lang rexx>/*REXX pgm does numerical integration using Gauss-Legendre Quadrature. */
parse arg digs .; if digs== then digs=40 /*assume the DIGS default?*/
numeric digits digs*2+10 /*use higher working DIGs.*/
pi=pi(); a=-3; b=3; bma =b-a; bpa =b+a; tiny ='1E-' || (digs*2)
true = exp(b)-exp(a); bmaH=bma/2; bpaH=bpa/2; times= digs%2 + 1
numeric digits digs+10 /*use lower working DIGITs*/
say 'step' center("iterative value",digs+3) ' difference ' /*show hdr*/
sep='────' copies('─' ,digs+3) '────────────'; say sep
do step=1 for times; p0z=1; p0.1=1; step_=step+.5
p1z=2; p1.1=1; p1.2=0; r.=0
/*█*/ do k=2 to step; km=k-1 /*█*/ do L=1 for p1z; T.L=p1.L /*█*/ end /*L*/ /*█*/ T.L=0; TT.=0 /*█*/ do L=1 for p0z; L2=L+2; TT.L2=p0.L /*█*/ end /*L*/ /*█*/ /*█*/ do j=1 for p1z+1;T.j=((k+km)*T.j-km*TT.j)/k; end /*j*/ /*█*/ p0z=p1z; do j=1 for p0z; p0.j=p1.j ; end /*j*/ /*█*/ p1z=p1z+1; do j=1 for p1z; p1.j= T.j ; end /*j*/ /*█*/ end /*k*/
/*▓*/ do !=1 for step
/*▓*/ x=cos(pi*(!-.25)/step_)
/*▓*/ do times%2 until abs(dx)<tiny
/*▓*/ f=p1.1; df=0
/*▓*/ do k=2 to p1z
/*▓*/ df=f+x*df
/*▓*/ f=p1.k+x*f
/*▓*/ end /*k*/
/*▓*/ dx=f/df; x=x-dx
/*▓*/ end /*times%2 until···*/
/*▓*/ r.1.!=x
/*▓*/ r.2.!=2/((1-x*x)*df*df)
/*▓*/ end /*!*/
sum=0
/*▒*/ do m=1 for step; sum=sum+r.2.m*exp(bpaH+r.1.m*bmaH)
/*▒*/ end /*m*/
z=bmaH*sum
say right(step,4) format(z,2,digs) translate(format(z-true,3,4,,0),'e',"E")
end /*step*/
say sep; say left(,4) true " {exact}" exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────COS subroutine──────────────────────*/ cos: procedure; x=r2r(arg(1)); _=1; z=1; y=x*x
do k=2 by 2 until p==z; p=z; _=-_*y/(k*(k-1)); z=z+_; end; return z
/*──────────────────────────────────E subroutine──────────────────────────────────────*/ e: return 2.7182818284590452353602874713526624977572470936999595749669676277240766303535 /*──────────────────────────────────EXP subroutine──────────────────────*/ exp: procedure; parse arg x; ix=x%1; if abs(x-ix)>.5 then ix=ix+sign(x)
x=x-ix; z=1; _=1; do j=1 until p==z; p=z; _=_*x/j; z=z+_; end
if z\==0 then z=z*e()**ix; return z
/*──────────────────────────────────PI subroutine──────────────────────────────────────*/ pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862 /*──────────────────────────────────R2R subroutine──────────────────────*/ r2r: return arg(1)//(2*pi()) /*normalize radians►1 unit circle*/</lang> output
step iterative value difference
──── ─────────────────────────────────────────── ────────────
1 6.0000000000000000000000000000000000000000 -1.4036e+1
2 17.4874646410555689643606840462449458421154 -2.5483
3 19.8536919968055821921309108927158495960775 -1.8206e-1
4 20.0286883952907008527738054439857661647073 -7.0615e-3
5 20.0355777183855621539285357252750939315016 -1.7214e-4
6 20.0357469750923438830654575585499253741530 -2.8797e-6
7 20.0357498197266007755718729372891903369401 -3.5093e-8
8 20.0357498544945172882260918041683132616237 -3.2529e-10
9 20.0357498548174338368864419454858704839263 -2.3700e-12
10 20.0357498548197898711175766908543458234008 -1.3927e-14
11 20.0357498548198037305529147159697031241994 -6.7396e-17
12 20.0357498548198037976759531014454017742327 -2.7323e-19
13 20.0357498548198037979482458119092690701863 -9.4143e-22
14 20.0357498548198037979491844483599375945130 -2.7906e-24
15 20.0357498548198037979491872317401917248453 -7.1915e-27
16 20.0357498548198037979491872389153958789316 -1.6260e-29
17 20.0357498548198037979491872389316236038179 -3.2517e-32
18 20.0357498548198037979491872389316560624361 -5.7920e-35
19 20.0357498548198037979491872389316561202637 -9.2480e-38
20 20.0357498548198037979491872389316561203561 -1.3328e-40
21 20.0357498548198037979491872389316561203562 5.6562e-44
──── ─────────────────────────────────────────── ────────────
20.0357498548198037979491872389316561203562082463657269288113065020927852103607165290016530 {exact}
Tcl
<lang tcl>package require Tcl 8.5 package require math::special package require math::polynomials package require math::constants math::constants::constants pi
- Computes the initial guess for the root i of a n-order Legendre polynomial
proc guess {n i} {
global pi
expr { cos($pi * ($i - 0.25) / ($n + 0.5)) }
}
- Computes and evaluates the n-order Legendre polynomial at the point x
proc legpoly {n x} {
math::polynomials::evalPolyn [math::special::legendre $n] $x
}
- Computes and evaluates the derivative of an n-order Legendre polynomial at point x
proc legdiff {n x} {
expr {$n / ($x**2 - 1) * ($x * [legpoly $n $x] - [legpoly [incr n -1] $x])}
}
- Computes the n nodes for an n-point quadrature rule. (i.e. n roots of a n-order polynomial)
proc nodes n {
set x [lrepeat $n 0.0]
for {set i 0} {$i < $n} {incr i} {
set val [guess $n [expr {$i + 1}]] foreach . {1 2 3 4 5} { set val [expr {$val - [legpoly $n $val] / [legdiff $n $val]}] } lset x $i $val
} return $x
}
- Computes the weight for an n-order polynomial at the point (node) x
proc legwts {n x} {
expr {2.0 / (1 - $x**2) / [legdiff $n $x]**2}
}
- Takes a array of nodes x and computes an array of corresponding weights w
proc weights x {
set n [llength $x]
set w {}
foreach xi $x {
lappend w [legwts $n $xi]
} return $w
}
- Integrates a lambda term f with a n-point Gauss-Legendre quadrature rule over the interval [a,b]
proc gausslegendreintegrate {f n a b} {
set x [nodes $n]
set w [weights $x]
set rangesize2 [expr {($b - $a)/2}]
set rangesum2 [expr {($a + $b)/2}]
set sum 0.0
foreach xi $x wi $w {
set y [expr {$rangesize2*$xi + $rangesum2}] set sum [expr {$sum + $wi*[apply $f $y]}]
}
expr {$sum * $rangesize2}
}</lang> Demonstrating: <lang tcl>puts "nodes(5) = [nodes 5]" puts "weights(5) = [weights [nodes 5]]" set exp {x {expr {exp($x)}}} puts "int(exp,-3,3) = [gausslegendreintegrate $exp 5 -3 3]"</lang>
- Output:
nodes(5) = 0.906179845938664 0.5384693101056831 -1.198509146801203e-94 -0.5384693101056831 -0.906179845938664 weights(5) = 0.2369268850561896 0.4786286704993664 0.5688888888888889 0.4786286704993664 0.2369268850561896 int(exp,-3,3) = 20.03557771838559
Ursala
using arbitrary precision arithmetic <lang Ursala>#import std
- import nat
legendre = # takes n to the pair of functions (P_n,P'_n), where P_n is the Legendre polynomial of order n
~&?\(1E0!,0E0!)! -+
^|/~& //mp..vid^ mp..sub\1E0+ mp..sqr,
~~ "c". ~&\1E0; ~&\"c"; ~&ar^?\0E0! mp..add^/mp..mul@alrPrhPX ^|R/~& ^|\~&t ^/~&l mp..mul,
@iiXNX ~&rZ->r @l ^/^|(~&tt+ sum@NNiCCiX+ successor,~&) both~&g&&~&+ -+
~* mp..zero_p?/~& (&&~&r ~&EZ+ ~~ mp..prec)^/~& ^(~&,..shr\8); mp..equ^|(~&,..gro\8)->l @r ^/~& ..shr\8,
^(~&rl,mp..mul*lrrPD)^/..nat2mp@r -+
^(~&l,mp..sub*+ zipp0E0^|\~& :/0E0)+ ~&rrt->lhthPX ^(
^lrNCC\~&lh mp..vid^*D/..nat2mp@rl -+
mp..sub*+ zipp0E0^|\~& :/0E0,
mp..mul~*brlD^|bbI/~&hthPX @l ..nat2mp~~+ predecessor~~NiCiX+-,
@r ^|/successor predecessor),
^|(mp..grow/1E0; @iNC ^lrNCC\~& :/0E0,~&/2)+-+-+-
nodes = # takes precision and order (p,n) to a list of nodes and weights <(x_1,w_1)..(x_n,w_n)>
-+
^H(
@lrr *+ ^/~&+ mp..div/( ..nat2mp 2)++ mp..mul^/(mp..sqr; //mp..sub ..nat2mp 1)+ mp..sqr+,
mp..shr^*DrlXS/~&ll ^|H\~& *+ @NiX+ ->l^|(~&lZ!|+ not+ //mp..eq,@r+ ^/~&+ mp..sub^/~&+ mp..div^)),
^/^|(~&,legendre) mp..cos*+ mp..mul^*D(
mp..div^|/mp..pi@NiC mp..add/5E-1+ ..nat2mp,
@r mp..bus/*2.5E-1+ ..nat2mp*+ nrange/1)+-
integral = # takes precision and order (p,n) to a function taking a function and interval (f,(a,b))
("p","n"). -+
mp..shrink^/~& difference\"p"+ mp..prec,
mp..mul^|/~& mp..add:-0E0+ * mp..mul^/~&rr ^H/~&ll mp..add^\~&lrr mp..mul@lrPrXl,
^(~&rl,-*nodes("p","n"))^|/~& mp..vid~~G/2E0+ ^/mp..bus mp..add+-</lang>
demonstration program:<lang Ursala>#show+
demo =
~&lNrCT (
^|lNrCT(:/'nodes:',:/'weights:')@lSrSX ..mp2str~~* nodes/160 5, :/'integral:' ~&iNC ..mp2str integral(160,5) (mp..exp,-3E0,3E0))</lang>
- Output:
nodes: 9.0617984593866399279762687829939296512565191076233E-01 5.3846931010568309103631442070020880496728660690555E-01 0.0000000000000000000000000000000000000000000000000E+00 -5.3846931010568309103631442070020880496728660690555E-01 -9.0617984593866399279762687829939296512565191076233E-01 weights: 2.3692688505618908751426404071991736264326000220463E-01 4.7862867049936646804129151483563819291229555334456E-01 5.6888888888888888888888888888888888888888888888896E-01 4.7862867049936646804129151483563819291229555334456E-01 2.3692688505618908751426404071991736264326000220463E-01 integral: 2.0035577718385562153928535725275093931501627207110E+01