Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

 

{|border=1 cellspacing=0 cellpadding=3

|<math>\int_{-1}^1 f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)</math>

|}

 

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 [[Numerical Integration|simple numerical integration methods]].

|<math>\int_a^b f(x)\,dx \approx \frac{b-a}{2} \sum_{i=1}^n w_i f\left(\frac{b-a}{2}x_i + \frac{a+b}{2}\right)</math>

|}

 

'''Task description'''

Similar to the task [[Numerical Integration]], the task here is to calculate the definite integral of a function <math>f(x)</math>, but by applying an n-point Gauss-Legendre quadrature rule, as described [[wp:Gaussian Quadrature|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.

 

 

 

=={{header|C}}==

#include <math.h>

 

x1 = x;

x -= lege_eval(N, x) / lege_diff(N, x);

lroots[i - 1] = x;

 

lege_inte(exp, -3, 3), exp(3) - exp(-3));

return 0;

{{out}}

<pre>Roots: 0.90618 0.538469 0 -0.538469 -0.90618

compred to actual

20.03574985</pre>

 

=={{header|Common Lisp}}==

(defun guess (n i)

(cos (* pi

(funcall f (+ (* (/ (- b a) 2.0d0)

(aref x i))

{{out|Example}}

#(0.906179845938664d0 0.5384693101056831d0 2.996272867003007d-95 -0.5384693101056831d0 -0.906179845938664d0)

 

 

(int #'exp 5 -3 3)

Comparison of the 5-point rule with simpler, but more costly methods from the task [[Numerical Integration]]:

0.24999999999999997d0

 

 

(int #'(lambda (x) x) 5 0 6000)

 

=={{header|D}}==

{{trans|C}}

 

immutable struct GaussLegendreQuadrature(size_t N, FP=double,

alias FP[N + 1][N + 1] CoefMat;

 

static FP legendreEval(in ref FP[N + 1][N + 1] lcoef,

in int n, in FP x) pure nothrow {

 

static FP legendreDiff(in ref CoefMat lcoef,

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

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

(x ^^ 2 - 1);

}

 

CoefMat lcoef = 0.0;

legendreCoefInit(/*ref*/ lcoef);

}

 

static private void legendreCoefInit(ref CoefMat lcoef)

lcoef[0][0] = lcoef[1][1] = 1;

foreach (immutable int n; 2 .. N + 1) { // n must be signed.

}

 

immutable FP c1 = (b - a) / 2;

immutable FP c2 = (b + a) / 2;

glq.integrate(&exp, -3, 3));

writefln("Compred to actual: %10.12f",

{{out}}

<pre>Roots: [0.90618, 0.538469, 0, -0.538469, -0.90618]

Integrating exp(x) over [-3, 3]: 20.035577718386

Compred to actual: 20.035749854820</pre>

 

=={{header|Fortran}}==

! -assume realloc_lhs

! when compiled with Intel Fortran.

end function

end program

 

<pre>

20 20.0357498548198037979491872388495 -.82E-28

</pre>

 

=={{header|Go}}==

Implementation pretty much by the methods given in the task description.

 

import (

}

panic("no convergence")

{{out}}

<pre>

integral: 20.035577718385564

</pre>

 

=={{header|J}}==

'''Solution:'''

)

 

 

:

{{out|Example use}}

 

code assumes function to be integrated has attribute Listable which is true of most built in Mathematica functions

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]]

{{out}}<pre>20.0356</pre>

 

=={{header|Maxima}}==

p: expand(legendre_p(n, x)),

q: expand(n/2*diff(p, x)*legendre_p(n - 1, x)),

% - bfloat(integrate(exp(x), x, -3, 3));

 

=={{header|OCaml}}==

| 0 -> 1.0

| 1 -> x

let f1 x = f ((x*.(b-.a) +. a +. b)*.0.5) in

let eval s (x,w) = s +. w*.(f1 x) in

which can be used in:

Printf.printf

"Gauss-Legendre %2d-point quadrature for exp over [-3..3] = %.16f\n"

calc 10;;

calc 15;;

{{out}}

<pre>

</pre>

This shows convergence to the correct double-precision value of the integral

although going beyond 20 points starts reducing the accuracy, due to accumulated rounding errors.

 

=={{header|ooRexx}}==

* 31.10.2013 Walter Pachl Translation from REXX (from PL/I)

* using ooRexx' rxmath package

Return

 

Output:

<pre> 1 6.0000000000000000000000000000000000000000 -1.4036E+1

20.03574985481980606 (exact)</pre>

 

 

 

 

 

}

 

 

($m1 - $x * $m0) * $n / (1 - $x**2);

}

 

}

 

}

$r;

}

 

}

 

 

}

}

 

}

 

{{out}}

<pre>Gauss-Legendre 5-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0355777183856

Gauss-Legendre 10-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0357498548198

Gauss-Legendre 20-point quadrature ∫₋₃⁺³ exp(x) dx ≈ 20.0357498548198</pre>

 

=={{header|PL/I}}==

Translated from Fortran.

Integration_Gauss: procedure options (main);

 

end gaussquad;

end Integration_Gauss;

<pre>

1 6.0000000000000000 -1.40E+0001

 

=={{header|Python}}==

##################################################################

print "Integral : ", ans

else:

{{out}}

<pre>

Integral : 20.0355777184

</pre>

 

=={{header|Racket}}==

Computation of the Legendre polynomials and derivatives:

 

(define (LegendreP n x)

(let compute ([n n] [Pn-1 x] [Pn-2 1])

(- (* x (LegendreP n x))

(LegendreP (- n 1) x))))

 

Computation of the Legendre polynomial roots:

 

(define (LegendreP-root n i)

; newton-raphson step

x′

(next (newton-step x′) x′)))))

 

Computation of Gauss-Legendre nodes and weights

 

(define (Gauss-Legendre-quadrature n)

;; positive roots

(if (odd? n) (list (/ 2 (sqr (LegendreP′ n 0)))) '())

(reverse weights))))

 

Integration using Gauss-Legendre quadratures:

 

(define (integrate f a b #:nodes (n 5))

(define m (/ (+ a b) 2))

(define (g x) (f (+ m (* d x))))

(* d (+ (apply + (map * w (map g x))))))

 

Usage:

 

> (Gauss-Legendre-quadrature 5)

'(-0.906179845938664 -0.5384693101056831 0 0.5384693101056831 0.906179845938664)

> (- (exp 3) (exp -3)

20.035749854819805

 

Accuracy of the method:

 

> (require plot)

> (parameterize ([plot-x-label "Number of Gaussian nodes"]

(list n (abs (- (integrate exp -3 3 #:nodes n)

(- (exp 3) (exp -3)))))))))

[[File:gauss.png]]

 

=={{header|REXX}}==

===version 1===

* 31.10.2013 Walter Pachl Translation from PL/I

* 01.11.2014 -"- see Version 2 for improvements

End

Numeric Digits (prec)

Output:

<pre> 1 6.0000000000000000000000000000000000000000 -1.4036E+1

 

===version 2===

 

 

 

 

</pre>

 

=={{header|Tcl}}==

{{tcllib|math::polynomials}}

{{tcllib|math::special}}

package require math::special

package require math::polynomials

}

expr {$sum * $rangesize2}

Demonstrating:

puts "weights(5) = [weights [nodes 5]]"

set exp {x {expr {exp($x)}}}

{{out}}

<pre>

=={{header|Ursala}}==

using arbitrary precision arithmetic

#import nat

 

mp..shrink^/~& difference\"p"+ mp..prec,

mp..mul^|/~& mp..add:-0E0+ * mp..mul^/~&rr ^H/~&ll mp..add^\~&lrr mp..mul@lrPrXl,

 

demo =

~&lNrCT (

^|lNrCT(:/'nodes:',:/'weights:')@lSrSX ..mp2str~~* nodes/160 5,

{{out}}

<pre>

integral:

2.0035577718385562153928535725275093931501627207110E+01</pre>

 

{{omit from|GUISS}}