Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

{{trans|Nim}}

F prev1(idx, pn1)

R (2 * idx - 1) * @x * pn1

R result * sum

{{out}}

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

#include <iomanip>

#include <cmath>

std::cout << "Integrating Exp(X) over [-3, 3]: " << gl5.integrate(-3., 3., RosettaExp) << std::endl;

std::cout << "Actual value: " << RosettaExp(3) - RosettaExp(-3) << std::endl;

{{out}}

Derived from the C++ and Java versions here.

using System;

//Works in .NET 6+

{{out}}

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

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

getLegendreCoeffs=: verb define M.

if. y<:1 do. 1 {.~ - y+1 return. end.

'nodes wgts'=. getGaussLegendrePoints x

-: (-~/ y) * wgts +/@:* u -: nodes p.~ (+/ , -~/) y

{{out|Example use}}

20.0356

-~/ ^ _3 3 NB. true value

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

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

=={{header|Mathematica}}/{{header|Wolfram Language}}==

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.

%Print the result.

disp(GLGD_int(@(x) exp(x), -3, 3, 5));

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

{{works with|Free Pascal|3.0.4}}

{{works with|Multics Pascal|8.0.4a}}

const Order = 5;

Writeln('Integrating Exp(x) over [-3, 3]: ',LegInt(-3,3):13:10);

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

{{out}}

<pre>

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

<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>

<span style="color: #004080;">integer</span> <span style="color: #000000;">order</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>

<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)-</span><span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)})</span>

<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" compared to actual = %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">res</span><span style="color: #0000FF;">})</span>

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

===With library routine===

One can also use the already invented wheel in NumPy:

# func is a function that takes a list-like input values

for d in range(3, 10):

{{out}}

<pre>3 19.853691996805587

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.

pi= pi(); digs= length(pi) - length(.); numeric digits digs; reps= digs % 2

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

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>

It is about twice as slow as version 2, &nbsp; due to the doubling of the number of decimal digits &nbsp; (precision).

pi= pi(); digs= length(pi) - length(.); numeric digits digs; reps= digs % 2

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

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

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

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

{{trans|Raku}}

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>

{{trans|C}}

{{libheader|Wren-fmt}}

var N = 5

var actual = 3.exp - (-3).exp

Fmt.print("\nIntegrating exp(x) over [-3, 3]:\n\t$10.8f,\n" +

{{out}}

=={{header|zkl}}==

{{trans|Raku}}

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>