Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

 

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

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

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

}

sub approximate-legendre-root(Int $n, Int $k) {

# Approximation due to Francesco Tricomi

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

}

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

take 0 => weight($n, 0) if $n !%% 2;

for 1 .. $n div 2 {

}

}

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 ≈ ",

 

{{out}}