Numerical integration/Gauss-Legendre Quadrature

In a general Gaussian quadrature rule, an definite integral is first approximated over the interval $[-1,1]$ by a polynomial approximable function $g(x)$ and a weighting function $W(x)$ :

\int _{-1}^{1}f(x)\,dx=\int _{-1}^{1}W(x)g(x)\,dx

In the case of Gauss-Legendre quadrature, the weighting function $W(x)=1$ , so we can easily approximate an integral of $f(x)$ by a sum of function values at specified points $x_{i}$ multiplied by some weights $w_{i}$ :

\int _{-1}^{1}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i})

The $n$ evaluation points $x_{i}$ for a n-point rule, also called "nodes", are roots of n-th order Legendre Polynomials $P_{n}(x)$ . Legendre polynomials are defined by the following recursive rule:

P_{0}(x)=1

P_{1}(x)=x

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)={\frac {n}{x^{2}}}\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}-{\frac {f(x_{n})}{f'(x_{n})}}

The first guess $x_{0}$ for the $i$ -th root of a $n$ -order polynomial $P_{n}$ can be given by

x_{0}=\cos \left(\pi \,{\frac {i-{\frac {1}{4}}}{n+{\frac {1}{2}}}}\right)

After we get the nodes $x_{i}$ , we compute the appropriate weights by:

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

\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)$ , 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

\int _{-3}^{3}\exp(x)\,dx\approx \sum _{i=1}^{5}w_{i}\;\exp(x_{i})\approx 20.036

Common Lisp

Computes the initial guess for the root i of a n-order Legendre polynomial.

(defun guess (n i)

 (cos (* pi
         (/ (- i 0.25d0)
            (+ n 0.5d0)))))

Computes and evaluates the n-order Legendre polynomial at the point x.

(defun legpoly (n x)

 (let ((pa 1.0d0)
       (pb x)
       (pn))
   (cond ((= n 0) pa)
         ((= n 1) pb)
         (t (loop for i from 2 to n do
                 (setf pn (- (* (/ (- (* 2 i) 1) i) x pb) 
                             (* (/ (- i 1) i) pa)))
                 (setf pa pb)
                 (setf pb pn)
                 finally (return pn))))))

Computes and evaluates the derivative of an n-order Legendre polynomial at point x.

(defun legdiff (n x)

 (* (/ n (- (* x x) 1))
    (- (* x (legpoly n x))
       (legpoly (- n 1) x))))

Computes the n nodes for an n-point quadrature rule. (i.e. n roots of a n-order polynomial)

(defun nodes (n)

 (let ((x (make-array n :initial-element 0.0d0)))
   (loop for i from 0 to (- n 1) do
        (let ((val (guess n (+ i 1))) ;Nullstellen-Schätzwert.
              (itermax 5))     
          (dotimes (j itermax)
            (setf val (- val
                         (/ (legpoly n val)
                            (legdiff n val)))))
          (setf (aref x i) val)))
   x))

Computes the weight for an n-order polynomial at the point (node) x.

(defun legwts (n x)

 (/ 2 
    (* (- 1 (* x x))
       (expt (legdiff n x) 2))))

Takes a array of nodes x and computes an array of corresponding weights w.

(defun weights (x)

 (let* ((n (car (array-dimensions x)))
        (w (make-array n :initial-element 0.0d0)))
   (loop for i from 0 to (- n 1) do
        (setf (aref w i) (legwts n (aref x i))))
   w))

Integrates a function f with a n-point Gauss-Legendre quadrature rule over the interval [a,b].

(defun int (f n a b)

 (let* ((x (nodes n))
        (w (weights x)))
   (* (/ (- b a) 2.0d0)
      (loop for i from 0 to (- n 1)
         sum (* (aref w i)
                (funcall f (+ (* (/ (- b a) 2.0d0)
                                 (aref x i))
                              (/ (+ a b) 2.0d0))))))))

</lang>

Example:

<lang lisp> (nodes 5)

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

(weights (nodes 5))

(0.23692688505618917d0 0.47862867049936647d0 0.5688888888888889d0 0.47862867049936647d0 0.23692688505618917d0)

(int #'exp 5 -3 3) 20.035577718385568d0 </lang>