user
|  |  |
| --- | --- |
| In a general Gaussian quadrature rule, an definite integral of       f ( x ) {\displaystyle {\displaystyle f(x)}}  is first approximated over the interval       [ − 1 , 1 ] {\displaystyle {\displaystyle [-1,1]}}  by a polynomial approximable function       g ( x ) {\displaystyle {\displaystyle g(x)}}  and a known weighting function       W ( x ) {\displaystyle {\displaystyle W(x)}} . | ∫  − 1  1 f ( x )  d x =  ∫  − 1  1 W ( x ) g ( x )  d x {\displaystyle {\displaystyle \int \_{-1}^{1}f(x)\,dx=\int \_{-1}^{1}W(x)g(x)\,dx}} |
| Those are then approximated by a sum of function values at specified points        x  i {\displaystyle {\displaystyle x\_{i}}}  multiplied by some weights        w  i {\displaystyle {\displaystyle w\_{i}}} : | ∫  − 1  1 W ( x ) g ( x )  d x ≈  ∑  i = 1  n  w  i g (  x  i ) {\displaystyle {\displaystyle \int \_{-1}^{1}W(x)g(x)\,dx\approx \sum \_{i=1}^{n}w\_{i}g(x\_{i})}} |
| In the case of Gauss-Legendre quadrature, the weighting function       W ( x ) = 1 {\displaystyle {\displaystyle W(x)=1}} , so we can approximate an integral of       f ( x ) {\displaystyle {\displaystyle f(x)}}  with: | ∫  − 1  1 f ( x )  d x ≈  ∑  i = 1  n  w  i f (  x  i ) {\displaystyle {\displaystyle \int \_{-1}^{1}f(x)\,dx\approx \sum \_{i=1}^{n}w\_{i}f(x\_{i})}} |

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 simple numerical integration methods.

|  |  |
| --- | --- |
| The       n {\displaystyle {\displaystyle n}}  evaluation points        x  i {\displaystyle {\displaystyle x\_{i}}}  for a n-point rule, also called "nodes", are roots of n-th order Legendre Polynomials        P  n ( x ) {\displaystyle {\displaystyle P\_{n}(x)}} . Legendre polynomials are defined by the following recursive rule: | P  0 ( x ) = 1 {\displaystyle {\displaystyle P\_{0}(x)=1}}   P  1 ( x ) = x {\displaystyle {\displaystyle P\_{1}(x)=x}}         n  P  n ( x ) = ( 2 n − 1 ) x  P  n − 1 ( x ) − ( n − 1 )  P  n − 2 ( x ) {\displaystyle {\displaystyle nP\_{n}(x)=(2n-1)xP\_{n-1}(x)-(n-1)P\_{n-2}(x)}} |
| There is also a recursive equation for their derivative: | P  n ′ ( x ) =   n   x  2 − 1  (  x  P  n ( x ) −  P  n − 1 ( x ) ) {\displaystyle {\displaystyle P\_{n}'(x)={\frac {n}{x^{2}-1}}\left(xP\_{n}(x)-P\_{n-1}(x)\right)}} |
| The roots of those polynomials are in general not analytically solvable, so they have to be approximated numerically, for example by Newton-Raphson iteration: | x  n + 1 =  x  n −    f (  x  n )   f ′ (  x  n ) {\displaystyle {\displaystyle x\_{n+1}=x\_{n}-{\frac {f(x\_{n})}{f'(x\_{n})}}}} |
| The first guess        x  0 {\displaystyle {\displaystyle x\_{0}}}  for the       i {\displaystyle {\displaystyle i}} -th root of a       n {\displaystyle {\displaystyle n}} -order polynomial        P  n {\displaystyle {\displaystyle P\_{n}}}  can be given by | x  0 = cos ⁡  (  π     i −   1 4  n +   1 2 ) {\displaystyle {\displaystyle x\_{0}=\cos \left(\pi \,{\frac {i-{\frac {1}{4}}}{n+{\frac {1}{2}}}}\right)}} |
| After we get the nodes        x  i {\displaystyle {\displaystyle x\_{i}}} , we compute the appropriate weights by: | w  i =   2   (  1 −  x  i  2 ) [  P  n ′ (  x  i )  ]  2 {\displaystyle {\displaystyle w\_{i}={\frac {2}{\left(1-x\_{i}^{2}\right)[P'\_{n}(x\_{i})]^{2}}}}} |
| After we have the nodes and the weights for a n-point quadrature rule, we can approximate an integral over any interval       [ a , b ] {\displaystyle {\displaystyle [a,b]}}  by | ∫  a  b f ( x )  d x ≈    b − a 2  ∑  i = 1  n  w  i f  (     b − a 2  x  i +    a + b 2 ) {\displaystyle {\displaystyle \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)}} |

**Task description**

Similar to the task Numerical Integration, the task here is to calculate the definite integral of a function 

f
(
x
)
{\displaystyle {\displaystyle f(x)}}
, but by applying an n-point Gauss-Legendre quadrature rule, as described 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.

To demonstrate the calculation, compute the weights and nodes for an 5-point quadrature rule and then use them to compute:

```
         
  
    
      
        
          
            
              ∫
              
                −
                3
              
              
                3
              
            
            exp
            ⁡
            (
            x
            )
            
            d
            x
            ≈
            
              ∑
              
                i
                =
                1
              
              
                5
              
            
            
              w
              
                i
              
            
            
            exp
            ⁡
            (
            
              x
              
                i
              
            
            )
            ≈
            20.036
          
        
      
    
    {\displaystyle {\displaystyle \int _{-3}^{3}\exp(x)\,dx\approx \sum _{i=1}^{5}w_{i}\;\exp(x_{i})\approx 20.036}}
  


```
assistant
(formerly Perl 6)

**Works with**: rakudo version 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 quadrature 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.

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

```

Output:

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