Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

=={{header|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>

{{out|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>

=={{header|D}}==

{{trans|C}}

<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 @safe @nogc 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 @safe @nogc {

return n * (x * legendreEval(lcoef, n, x) -

legendreEval(lcoef, n - 1, x)) /

(x ^^ 2 - 1);

}

CoefMat lcoef = 0.0;

legendreCoefInit(/*ref*/ lcoef);

// Legendre roots:

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));

}

}

static private void legendreCoefInit(ref CoefMat lcoef)

pure nothrow @safe @nogc {

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(/*in*/ FP x) pure nothrow @safe @nogc f,

in FP a, in FP b)

pure nothrow @safe @nogc {

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 - exp(-3.0));

}</lang>

{{out}}

<pre>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</pre>

=={{header|PARI/GP}}==

{{works with|PARI/GP|2.4.2 and above}}

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>

{{out}}

<pre>time = 0 ms.

%1 = 20.035577718385562153928535725275093932 + 0.E-37*I</pre>

{{works with|PARI/GP|2.9.0 and above}}

Gauss-Legendre quadrature is built-in from 2.9 forward.

<lang parigp>intnumgauss(x=-3, 3, exp(x), intnumgaussinit(5))

intnumgauss(x=-3, 3, exp(x)) \\ determine number of points automatically; all digits shown should be accurate</lang>

{{out}}

<pre>%1 = 20.035746975092343883065457558549925374

%2 = 20.035749854819803797949187238931656120</pre>

{{out}}

<pre>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</pre>

=={{header|Perl 6}}==

{{works with|rakudo|2015-09-24}}

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 <tt>quadrature</tt> 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) {

flat 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 flat 5 .. 10, 20;</lang>