Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

<big><big><math>\int_{-3}^{3} \exp(x) \, dx \approx \sum_{i=1}^5 w_i \; \exp(x_i) \approx 20.036</math></big></big>

<br><br>

 

=={{header|Axiom}}==

{{trans|Maxima}}

RECORD ==> Record(x : List Fraction Integer, w : List Fraction Integer)

 

c := (a+b)/2

h := (b-a)/2

gaussIntegrate(4/(1+x^2), x=0..1, 20)

 

% - %pi

 

 

=={{header|C}}==

#include <math.h>

 

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

return 0;

{{out}}

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

 

Does not quite perform the task quite as specified since the node count, N, is set at compile time (to avoid heap allocation) so cannot be passed as a parameter.

namespace Rosetta {

 

out << ' ' << legpoly.root(i);

}

out << "Weights:";

for (int i = 0; i <= eDEGREE; ++i) {

out << ' ' << legpoly.weight(i);

}

}

private:

 

gl5.print_roots_and_weights(std::cout);

 

{{out}}

Actual value: 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,

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

3.0.exp - exp(-3.0));

{{out}}

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

=={{header|Delphi}}==

 

 

{$APPTYPE CONSOLE}

Writeln('Actual value: ',Exp(3)-Exp(-3):13:10);

Readln;

 

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

=={{header|Haskell}}==

Integration formula

where d = (b - a)/2

m = (b + a)/2

w x = 2/(1-x^2)/(legendreP' n x)^2

roots = map (findRoot (legendreP n) (legendreP' n) . x0) [1..n]

 

Calculation of Legendre polynomials

where go 0 p2 _ = p2

go 1 _ p1 = p1

go n p2 p1 = go (n-1) p1 $ ((2*n-1)*x*p1 - (n-1)*p2)/n

 

 

Universal auxilary functions

 

fixedPoint f x | abs (fx - x) < 1e-15 = x

| otherwise = fixedPoint f fx

 

Integration on a given mesh using Gauss-Legendre quadrature:

 

{{out}}

=={{header|J}}==

'''Solution:'''

)

 

 

:

{{out|Example use}}

 

=={{header|Java}}==

{{trans|C}}

{{works with|Java|8}}

import java.util.function.Function;

 

legeInte(x -> exp(x), -3, 3), exp(3) - exp(-3));

}

<pre>Roots: 0,906180 0,538469 0,000000 -0,538469 -0,906180

Weight: 0,236927 0,478629 0,568889 0,478629 0,236927

 

=={{header|JavaScript}}==

const factorial = n => n <= 1 ? 1 : n * factorial(n - 1);

const M = n => (n - (n % 2 !== 0)) / 2;

}

console.log(gaussLegendre(x => Math.exp(x), -3, 3, 5));

{{output}}

<pre>

20.035577718385575

</pre>

 

=={{header|Julia}}==

This function computes the points and weights of an ''N''-point Gauss–Legendre quadrature rule on the interval (''a'',''b''). It uses the O(''N''²) algorithm described in Trefethen & Bau, ''Numerical Linear Algebra'', which finds the points and weights by computing the eigenvalues and eigenvectors of a real-symmetric tridiagonal matrix:

 

function gauss(a, b, N)

λ, Q = eigen(SymTridiagonal(zeros(N), [n / sqrt(4n^2 - 1) for n = 1:N-1]))

@. (λ + 1) * (b - a) / 2 + a, [2Q[1, i]^2 for i = 1:N] * (b - a) / 2

(This code is a simplified version of the <code>Base.gauss</code> subroutine in the Julia standard library.)

{{out}}

=={{header|Kotlin}}==

{{trans|Java}}

 

class Legendre(val N: Int) {

println("compared to actual:")

println("\t%10.8f".format(exp(3.0) - exp(-3.0)))

{{Out}}

<pre>Roots: 0.906180 0.538469 0.000000 -0.538469 -0.906180

 

=={{header|Lua}}==

 

local legendreRoots = {}

do

print(gaussLegendreQuadrature(function(x) return math.exp(x) end, -3, 3, 5))

{{out}}<pre>20.035577718386</pre>

 

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|MATLAB}}==

Translated from the Python solution.

%Integration using Gauss-Legendre quad

%Does almost the same as 'integral' in MATLAB

end

end

{{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|Nim}}==

{{trans|Common Lisp}}

import math, strformat

 

 

main()

{{out}}

<pre>

 

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

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

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

{{out}}

<pre>time = 0 ms.

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

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

{{out}}

<pre>%1 = 20.035746975092343883065457558549925374

 

=={{header|Pascal}}==

 

=={{header|Perl}}==

{{trans|Raku}}

use constant pi => 3.14159265;

 

printf("Gauss-Legendre %2d-point quadrature ∫₋₃⁺³ exp(x) dx ≈ %.13f\n", $_, quadrature($_, -3, +3) )

for 5 .. 10, 20;

{{out}}

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

=={{header|Phix}}==

{{trans|Lua}}

{{out}}

<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}}==

{{libheader|NumPy}}

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

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

 

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.

 

multi legendre-pair(Int $n, $x) {

my ($m1, $m2) = legendre-pair($n - 1, $x);

say "Gauss-Legendre $_.fmt('%2d')-point quadrature ∫₋₃⁺³ exp(x) dx ≈ ",

 

{{out}}

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

<br>where there's practically no space on the right side of the REXX source statements. &nbsp; It presents a good

<br>visual indication of what's what, &nbsp; but it's the dickens to pay when updating the source code.

!.= .; b= 3; a= -b; bma= b - a; bmaH= bma / 2; tiny= '1e-'digs

trueV= exp(b)-exp(a); bpa= b + a; bpaH= bpa / 2

sep='──────' copies("─", digs+3) '─────────────'; say sep /* " sep*/

 

/*█*/ end /*k*/

$= 0

z= bmaH * $ /*calculate target value (Z)*/

Ndif= translate( format(dif, 3, 4, 2, 0), 'e', "E")

if #\==1 then say center(#, 6) z' ' Ndif /*no display if not computed*/

say 'Using ' digs " digit precision, the" ,

'N-point Gauss─Legendre quadrature (GLQ) had an accuracy of ' xdif-2 " digits."

/*───────────────────────────────────────────────────────────────────────────────────────────*/

e: return 2.718281828459045235360287471352662497757247093699959574966967627724076630353547595

/*───────────────────────────────────────────────────────────────────────────────────────────*/

exp: procedure; parse arg x; ix= x % 1; if abs(x-ix)>.5 then ix= ix + sign(x); x= x-ix; z= 1

{{out|output|text=&nbsp; when using the default inputs:}}

<pre>

 

===version 3, more precision===

 

!.= .; b= 3; a= -b; bma= b - a; bmaH= bma / 2; tiny= '1e-'digs

trueV= exp(b)-exp(a); bpa= b + a; bpaH= bpa / 2

sep='──────' copies("─", digs+3) '─────────────'; say sep /* " sep*/

 

/*█*/ end /*k*/

$= 0

z= bmaH * $ /*calculate target value (Z)*/

Ndif= translate( format(dif, 3, 4, 3, 0), 'e', "E")

if #\==1 then say center(#, 6) z' ' Ndif /*no display if not computed*/

say 'Using ' digs " digit precision, the" ,

'N-point Gauss─Legendre quadrature (GLQ) had an accuracy of ' xdif-2 " digits."

/*───────────────────────────────────────────────────────────────────────────────────────────*/

e: return 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759,

/*───────────────────────────────────────────────────────────────────────────────────────────*/

pi: return 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899,

/*───────────────────────────────────────────────────────────────────────────────────────────*/

cos: procedure expose !.; parse arg x; if !.x\==. then return !.x; _= 1; z=1; y= x*x

/*───────────────────────────────────────────────────────────────────────────────────────────*/

exp: procedure; parse arg x; ix= x % 1; if abs(x-ix)>.5 then ix= ix + sign(x); x= x-ix; z= 1

{{out|output|text=&nbsp; when using the default inputs:}}

 

(Shown at about two-thirds size.)

<pre style="font-size:67%">

 

</pre>

 

=={{header|Scala}}==

{{Out}}Best seen in running your browser either by [https://scalafiddle.io/sf/rrvzhH1/0 ScalaFiddle (ES aka JavaScript, non JVM)] or [https://scastie.scala-lang.org/yYqRqizfSZip2DhYbdfZ2w Scastie (remote JVM)].

 

object GaussLegendreQuadrature extends App {

println(f"compared to actual%n\t${exp(3) - exp(-3)}%10.8f")

 

 

=={{header|Sidef}}==

func legendre_pair( n, x) {

var (m1, m2) = legendre_pair(n - 1, x)

printf("Gauss-Legendre %2d-point quadrature ∫₋₃⁺³ exp(x) dx ≈ %.15f\n",

i, quadrature(i, {.exp}, -3, +3))

{{out}}

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

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

 

=={{header|zkl}}==

if(n==1) return(x,1.0);

m1,m2:=legendrePair(n-1,x);

scale:='wrap(x){ (x*(b - a) + a + b) / 2 };

nds.reduce('wrap(p,[(r,w)]){ p + w*f(scale(r)) },0.0) * (b - a)/2

("Gauss-Legendre %2d-point quadrature "

"\U222B;\U208B;\U2083;\U207A;\UB3; exp(x) dx = %.13f")

.fmt(n,quadrature(n, fcn(x){ x.exp() }, -3, 3))

{{out}}

<pre>