Problem Statement

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.

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

f ( x )

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

Let's start with the solution:

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 */