Numerical integration/Gauss-Legendre Quadrature

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.

In a general Gaussian quadrature rule, an definite integral of $f(x)$ is first approximated over the interval $[-1,1]$ by a polynomial approximable function $g(x)$ and a known weighting function $W(x)$ .	$\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}$ multiplied by some weights $w_{i}$ :	$\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$ , so we can approximate an integral of $f(x)$ with:	$\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}-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}-{\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

C

<lang C>#include <stdio.h>

include <math.h>

define N 5

double Pi; double lroots[N]; double weight[N]; double lcoef[N + 1][N + 1] = Template:0;

void lege_coef() { int n, i; lcoef[0][0] = lcoef[1][1] = 1; for (n = 2; n <= N; n++) { lcoef[n][0] = -(n - 1) * lcoef[n - 2][0] / n; for (i = 1; i <= n; i++) lcoef[n][i] = ((2 * n - 1) * lcoef[n - 1][i - 1] - (n - 1) * lcoef[n - 2][i] ) / n; } }

double lege_eval(int n, double x) { int i; double s = lcoef[n][n]; for (i = n; i; i--) s = s * x + lcoef[n][i - 1]; return s; }

double lege_diff(int n, double x) { return n * (x * lege_eval(n, x) - lege_eval(n - 1, x)) / (x * x - 1); }

void lege_roots() { int i; double x, x1; for (i = 1; i <= N; i++) { x = cos(Pi * (i - .25) / (N + .5)); do { x1 = x; x -= lege_eval(N, x) / lege_diff(N, x); } while (x != x1); /* x != x1 is normally a no-no, but this task happens to be * well behaved. */ lroots[i - 1] = x;

x1 = lege_diff(N, x); weight[i - 1] = 2 / ((1 - x * x) * x1 * x1); } }

double lege_inte(double (*f)(double), double a, double b) { double c1 = (b - a) / 2, c2 = (b + a) / 2, sum = 0; int i; for (i = 0; i < N; i++) sum += weight[i] * f(c1 * lroots[i] + c2); return c1 * sum; }

int main() { int i; Pi = atan2(1, 1) * 4;

lege_coef(); lege_roots();

printf("Roots: "); for (i = 0; i < N; i++) printf(" %g", lroots[i]);

printf("\nWeight:"); for (i = 0; i < N; i++) printf(" %g", weight[i]);

printf("\nintegrating Exp(x) over [-3, 3]:\n\t%10.8f,\n" "compred to actual\n\t%10.8f\n", lege_inte(exp, -3, 3), exp(3) - exp(-3)); return 0; }</lang>output:

Roots:  0.90618 0.538469 0 -0.538469 -0.90618
Weight: 0.236927 0.478629 0.568889 0.478629 0.236927
integrating Exp(x) over [-3, 3]:
        20.03557772,
compred to actual
        20.03574985

Common Lisp

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

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

(defun legdiff (n x)

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

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

(defun legwts (n x)

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

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

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

Comparison of the 5-point rule with simpler, but more costly methods from the task Numerical Integration:

<lang lisp>(int #'(lambda (x) (expt x 3)) 5 0 1) 0.24999999999999997d0

(int #'(lambda (x) (/ 1 x)) 5 1 100) 4.059147508941519d0

(int #'(lambda (x) x) 5 0 5000) 1.25d7

(int #'(lambda (x) x) 5 0 6000) 1.8000000000000004d7</lang>

D

Translation of: C

<lang d>import std.stdio, std.math;

final /*const*/ class GaussLegendreQuadrature(size_t N, FP=double,

                                             size_t NBITS=50) {
   /*const*/ public double[N] lroots, weight;
   //FP[N + 1][N + 1] lcoef = 0.0;
   FP[N + 1][N + 1] lcoef;

   this() pure nothrow {
       foreach (ref row; lcoef)
           row[] = 0.0;
       legendreCoef();
       legendreRoots();
   }

   private void legendreCoef() pure nothrow {
       lcoef[0][0] = lcoef[1][1] = 1;
       foreach (int n; 2 .. N + 1) { // n must be signed
           lcoef[n][0] = -(n - 1) * lcoef[n - 2][0] / n;
           foreach (i; 1 .. n + 1)
               lcoef[n][i] = ((2 * n - 1) * lcoef[n - 1][i - 1] -
                              (n - 1) * lcoef[n - 2][i]) / n;
       }
   }

   private FP legendreEval(in int n, in FP x)
   const pure nothrow {
       FP s = lcoef[n][n];
       foreach_reverse (i; 1 .. n+1)
           s = s * x + lcoef[n][i - 1];
       return s;
   }

   private FP legendreDiff(in int n, in FP x)
   const pure nothrow {
       return n * (x * legendreEval(n, x) - legendreEval(n - 1, x))/
              (x ^^ 2 - 1);
   }

   private void legendreRoots() pure nothrow {
       foreach (i; 1 .. N + 1) {
           FP x = cos(PI * (i - 0.25) / (N + 0.5));
           FP x1;
           do {
               x1 = x;
               x -= legendreEval(N, x) / legendreDiff(N, x);
           } while (feqrel(x, x1) < NBITS);
           lroots[i - 1] = x;

           x1 = legendreDiff(N, x);
           weight[i - 1] = 2 / ((1 - x ^^ 2) * (x1 ^^ 2));
       }
   }

   public FP integrate(in FP function(FP) f,
                       in FP a, in FP b) const {
       immutable FP c1 = (b - a) / 2,
                    c2 = (b + a) / 2;
       FP sum = 0.0;
       foreach (i; 0 .. N)
           sum += weight[i] * f(c1 * lroots[i] + c2);
       return c1 * sum;
   }

}

void main() {

   const glq = new GaussLegendreQuadrature!(5, real);
   writeln("Roots:  ", glq.lroots);
   writeln("Weight: ", glq.weight);
   writefln("Integrating exp(x) over [-3, 3]: %10.12f",
            glq.integrate(&exp, -3, 3));
   writefln("Compred to actual:               %10.12f",
            exp(3.0) - exp(-3.0));

}</lang>

Roots:  [0.90618, 0.538469, 0, -0.538469, -0.90618]
Weight: [0.236927, 0.478629, 0.568889, 0.478629, 0.236927]
Integrating exp(x) over [-3, 3]: 20.035577718386
Compred to actual:               20.035749854820

Go

Implementation pretty much by the methods given in the task description. <lang go>package main

import (

   "fmt"
   "math"

)

// cFunc for continuous function. A type definition for convenience. type cFunc func(float64) float64

func main() {

   fmt.Println("integral:", glq(math.Exp, -3, 3, 5))

}

// glq integrates f from a to b by Guass-Legendre quadrature using n nodes. // For the task, it also shows the intermediate values determining the nodes: // the n roots of the order n Legendre polynomal and the corresponding n // weights used for the integration. func glq(f cFunc, a, b float64, n int) float64 {

   x, w := glqNodes(n, f)
   show := func(label string, vs []float64) {
       fmt.Printf("%8s: ", label)
       for _, v := range vs {
           fmt.Printf("%8.5f ", v)
       }
       fmt.Println()
   }
   show("nodes", x)
   show("weights", w)
   var sum float64
   bma2 := (b - a) * .5
   bpa2 := (b + a) * .5
   for i, xi := range x {
       sum += w[i] * f(bma2*xi+bpa2)
   }
   return bma2 * sum

}

// glqNodes computes both nodes and weights for a Gauss-Legendre // Quadrature integration. Parameters are n, the number of nodes // to compute and f, a continuous function to integrate. Return // values have len n. func glqNodes(n int, f cFunc) (node []float64, weight []float64) {

   p := legendrePoly(n)
   pn := p[n]
   n64 := float64(n)
   dn := func(x float64) float64 {
       return (x*pn(x) - p[n-1](x)) * n64 / (x*x - 1)
   }
   node = make([]float64, n)
   for i := range node {
       x0 := math.Cos(math.Pi * (float64(i+1) - .25) / (n64 + .5))
       node[i] = newtonRaphson(pn, dn, x0)
   }
   weight = make([]float64, n)
   for i, x := range node {
       dnx := dn(x)
       weight[i] = 2 / ((1 - x*x) * dnx * dnx)
   }
   return

}

// legendrePoly constructs functions that implement Lengendre polynomials. // This is done by function composition by recurrence relation (Bonnet's.) // For given n, n+1 functions are returned, computing P0 through Pn. func legendrePoly(n int) []cFunc {

   r := make([]cFunc, n+1)
   r[0] = func(float64) float64 { return 1 }
   r[1] = func(x float64) float64 { return x }
   for i := 2; i <= n; i++ {
       i2m1 := float64(i*2 - 1)
       im1 := float64(i - 1)
       rm1 := r[i-1]
       rm2 := r[i-2]
       invi := 1 / float64(i)
       r[i] = func(x float64) float64 {
           return (i2m1*x*rm1(x) - im1*rm2(x)) * invi
       }
   }
   return r

}

// newtonRaphson is general purpose, although totally primitive, simply // panicking after a fixed number of iterations without convergence to // a fixed error. Parameter f must be a continuous function, // df its derivative, x0 an initial guess. func newtonRaphson(f, df cFunc, x0 float64) float64 {

   for i := 0; i < 30; i++ {
       x1 := x0 - f(x0)/df(x0)
       if math.Abs(x1-x0) <= math.Abs(x0*1e-15) {
           return x1
       }
       x0 = x1
   }
   panic("no convergence")

}</lang> Output:

   nodes:  0.90618  0.53847  0.00000 -0.53847 -0.90618 
 weights:  0.23693  0.47863  0.56889  0.47863  0.23693 
integral: 20.035577718385564

Maxima

<lang maxima>gauss_coeff(n) := block([p, q, v, w], p: expand(legendre_p(n, x)), q: expand(n/2*diff(p, x)*legendre_p(n - 1, x)), v: map(rhs, bfallroots(p)), w: map(lambda([z], 1/subst([x = z], q)), v), [map(bfloat, v), map(bfloat, w)])$

gauss_int(f, a, b, n) := block([u, x, w, c, h], u: gauss_coeff(n), x: u[1], w: u[2], c: bfloat((a + b)/2), h: bfloat((b - a)/2), h*sum(w[i]*bfloat(f(c + x[i]*h)), i, 1, n))$

fpprec: 40$

gauss_int(lambda([x], 4/(1 + x^2)), 0, 1, 20); /* 3.141592653589793238462643379852215927697b0 */

% - bfloat(%pi); /* -3.427286956499858315999116083264403489053b-27 */

gauss_int(exp, -3, 3, 5); /* 2.003557771838556215392853572527509393154b1 */

% - bfloat(integrate(exp(x), x, -3, 3)); /* -1.721364342416440206515136565621888185351b-4 */</lang>

OCaml

<lang OCaml>let rec leg n x = match n with (* Evaluate Legendre polynomial *)

  | 0 -> 1.0
  | 1 -> x
  | k -> let u = 1.0 -. 1.0 /. float k in
     (1.0+.u)*.x*.(leg (k-1) x) -. u*.(leg (k-2) x);;

let leg' n x = match n with (* derivative *)

  | 0 -> 0.0
  | 1 -> 1.0
  | _ -> ((leg (n-1) x) -. x*.(leg n x)) *. (float n)/.(1.0-.x*.x);;

let approx_root k n = (* Reversed Francesco Tricomi: 1 <= k <= n *)

  let pi = acos (-1.0) and s = float(2*n)
  and t = 1.0 +. float(1-4*k)/.float(4*n+2) in
  (1.0 -. (float (n-1))/.(s*.s*.s))*.cos(pi*.t);;

let rec refine r n = (* Newton-Raphson *)

  let r1 = r -. (leg n r)/.(leg' n r) in
  if abs_float (r-.r1) < 2e-16 then r1 else refine r1 n;;

let root k n = refine (approx_root k n) n;;

let node k n = (* Abscissa and weight *)

  let r = root k n in
  let deriv = leg' n r in
  let w = 2.0/.((1.0-.r*.r)*.(deriv*.deriv)) in
  (r,w);;

let nodes n =

  let rec aux k = if k > n then [] else node k n :: aux (k+1)
  in aux 1;;

let quadrature n f a b =

  let f1 x = f ((x*.(b-.a) +. a +. b)*.0.5) in
  let eval s (x,w) = s +. w*.(f1 x) in
  0.5*.(b-.a)*.(List.fold_left eval 0.0 (nodes n));;</lang>

which can be used in: <lang OCaml>let calc n =

  Printf.printf
     "Gauss-Legendre %2d-point quadrature for exp over [-3..3] = %.16f\n"
     n (quadrature n exp (-3.0) 3.0);;

calc 5;; calc 10;; calc 15;; calc 20;;</lang> This outputs

Gauss-Legendre  5-point quadrature for exp over [-3..3] = 20.0355777183855608
Gauss-Legendre 10-point quadrature for exp over [-3..3] = 20.0357498548197839
Gauss-Legendre 15-point quadrature for exp over [-3..3] = 20.0357498548198052
Gauss-Legendre 20-point quadrature for exp over [-3..3] = 20.0357498548198052

showing convergence to the correct double-precision value of the integral <lang Ocaml>Printf.printf "%.16f\n" ((exp 3.0) -.(exp (-3.0)));; 20.0357498548198052</lang> although going beyond 20 points starts reducing the accuracy, due to accumulated rounding errors.

J

Solution: <lang j>P =: 3 :0 NB. list of coefficients for yth Legendre polynomial

  if. y<:1 do. 1{.~->:y return. end.
  y%~ (<:(,~+:)y) -/@:* (0,P<:y),:(P y-2)

)

getpoints =: 3 :0 NB. points,:weights for y points

  x=. 1{:: p. p=.P y
  w=. 2% (-.*:x)**:(p..p)p.x
  x,:w

)

GaussLegendre =: 1 :0 NB. npoints function GaussLegendre (a,b)

  'x w'=.getpoints x
  -:(-~/y)* +/w* u -:((+/,-~/)y)p.x

)</lang>

Example use: <lang j> 5 ^ GaussLegendre _3 3 20.0356</lang>

PARI/GP

Works with: PARI/GP version 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. <lang parigp>GLq(f,a,b,n)={

 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

GLq(x->exp(x), -3, 3, 5) \\ As of version 2.4.4, this can be written GLq(exp, -3, 3, 5)</lang>

Results:

time = 0 ms.
%1 = 20.035577718385562153928535725275093932 + 0.E-37*I

Tcl

Translation of: Common Lisp

Library: Tcllib (Package: math::constants)

Library: Tcllib (Package: math::polynomials)

Library: Tcllib (Package: math::special)

<lang tcl>package require Tcl 8.5 package require math::special package require math::polynomials package require math::constants math::constants::constants pi

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

proc guess {n i} {

   global pi
   expr { cos($pi * ($i - 0.25) / ($n + 0.5)) }

}

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

proc legpoly {n x} {

   math::polynomials::evalPolyn [math::special::legendre $n] $x

}

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

proc legdiff {n x} {

   expr {$n / ($x**2 - 1) * ($x * [legpoly $n $x] - [legpoly [incr n -1] $x])}

}

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

proc nodes n {

   set x [lrepeat $n 0.0]
   for {set i 0} {$i < $n} {incr i} {

set val [guess $n [expr {$i + 1}]] foreach . {1 2 3 4 5} { set val [expr {$val - [legpoly $n $val] / [legdiff $n $val]}] } lset x $i $val

   }
   return $x

}

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

proc legwts {n x} {

   expr {2.0 / (1 - $x**2) / [legdiff $n $x]**2}

}

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

proc weights x {

   set n [llength $x]
   set w {}
   foreach xi $x {

lappend w [legwts $n $xi]

   }
   return $w

}

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

proc gausslegendreintegrate {f n a b} {

   set x [nodes $n]
   set w [weights $x]
   set rangesize2 [expr {($b - $a)/2}]
   set rangesum2 [expr {($a + $b)/2}]
   set sum 0.0
   foreach xi $x wi $w {

set y [expr {$rangesize2*$xi + $rangesum2}] set sum [expr {$sum + $wi*[apply $f $y]}]

   }
   expr {$sum * $rangesize2}

}</lang> Demonstrating: <lang tcl>puts "nodes(5) = [nodes 5]" puts "weights(5) = [weights [nodes 5]]" set exp {x {expr {exp($x)}}} puts "int(exp,-3,3) = [gausslegendreintegrate $exp 5 -3 3]"</lang> Output:

nodes(5) = 0.906179845938664 0.5384693101056831 -1.198509146801203e-94 -0.5384693101056831 -0.906179845938664
weights(5) = 0.2369268850561896 0.4786286704993664 0.5688888888888889 0.4786286704993664 0.2369268850561896
int(exp,-3,3) = 20.03557771838559

Ursala

using arbitrary precision arithmetic <lang Ursala>#import std

import nat

legendre = # takes n to the pair of functions (P_n,P'_n), where P_n is the Legendre polynomial of order n

~&?\(1E0!,0E0!)! -+

  ^|/~& //mp..vid^ mp..sub\1E0+ mp..sqr,
  ~~ "c". ~&\1E0; ~&\"c"; ~&ar^?\0E0! mp..add^/mp..mul@alrPrhPX ^|R/~& ^|\~&t ^/~&l mp..mul,
  @iiXNX ~&rZ->r @l ^/^|(~&tt+ sum@NNiCCiX+ successor,~&) both~&g&&~&+ -+
     ~* mp..zero_p?/~& (&&~&r ~&EZ+ ~~ mp..prec)^/~& ^(~&,..shr\8); mp..equ^|(~&,..gro\8)->l @r ^/~& ..shr\8,
     ^(~&rl,mp..mul*lrrPD)^/..nat2mp@r -+
        ^(~&l,mp..sub*+ zipp0E0^|\~& :/0E0)+ ~&rrt->lhthPX ^(
           ^lrNCC\~&lh mp..vid^*D/..nat2mp@rl -+
              mp..sub*+ zipp0E0^|\~& :/0E0,
              mp..mul~*brlD^|bbI/~&hthPX @l ..nat2mp~~+ predecessor~~NiCiX+-,
           @r ^|/successor predecessor),
        ^|(mp..grow/1E0; @iNC ^lrNCC\~& :/0E0,~&/2)+-+-+-

nodes = # takes precision and order (p,n) to a list of nodes and weights <(x_1,w_1)..(x_n,w_n)>

-+

  ^H(
     @lrr *+ ^/~&+ mp..div/( ..nat2mp 2)++ mp..mul^/(mp..sqr; //mp..sub ..nat2mp 1)+ mp..sqr+,
     mp..shr^*DrlXS/~&ll ^|H\~& *+ @NiX+ ->l^|(~&lZ!|+ not+ //mp..eq,@r+ ^/~&+ mp..sub^/~&+ mp..div^)),
  ^/^|(~&,legendre) mp..cos*+ mp..mul^*D(
     mp..div^|/mp..pi@NiC mp..add/5E-1+ ..nat2mp,
     @r mp..bus/*2.5E-1+ ..nat2mp*+ nrange/1)+-

integral = # takes precision and order (p,n) to a function taking a function and interval (f,(a,b))

("p","n"). -+

  mp..shrink^/~& difference\"p"+ mp..prec,
  mp..mul^|/~& mp..add:-0E0+ * mp..mul^/~&rr ^H/~&ll mp..add^\~&lrr mp..mul@lrPrXl,
  ^(~&rl,-*nodes("p","n"))^|/~& mp..vid~~G/2E0+ ^/mp..bus mp..add+-</lang>

demonstration program:<lang Ursala>#show+

demo =

~&lNrCT (

  ^|lNrCT(:/'nodes:',:/'weights:')@lSrSX ..mp2str~~* nodes/160 5,
  :/'integral:' ~&iNC ..mp2str integral(160,5) (mp..exp,-3E0,3E0))</lang>

output:

nodes:
9.0617984593866399279762687829939296512565191076233E-01
5.3846931010568309103631442070020880496728660690555E-01
0.0000000000000000000000000000000000000000000000000E+00
-5.3846931010568309103631442070020880496728660690555E-01
-9.0617984593866399279762687829939296512565191076233E-01

weights:
2.3692688505618908751426404071991736264326000220463E-01
4.7862867049936646804129151483563819291229555334456E-01
5.6888888888888888888888888888888888888888888888896E-01
4.7862867049936646804129151483563819291229555334456E-01
2.3692688505618908751426404071991736264326000220463E-01

integral:
2.0035577718385562153928535725275093931501627207110E+01