# Unexpected result of summation

I wrote a small module that gives me an incorrect output-set. It should be a single number!

I don't understand what went wrong.

This is the form of summation used: $$\frac{1}{2} (b-a) \sum_{i=1}^n w_i\,f\left(\frac{b-a}{2}x_i+\frac{a+b}{2}\right)$$ where ${x,w}$ are generated by GaussianQuadratureWeights command:

  Needs["NumericalDifferentialEquationAnalysis"]
GaussLegendreQuadrature[a_, b_, n_] := Module[{weights, i},
(* GaussianQuadratureWeights[n, a, b] gives a list of the
n pairs {Subscript[x, i], Subscript[w, i]} of the elementary
n-point Gaussian formula for quadrature on the interval a to b,
where Subscript[w, i] is the weight of the abscissa
Subscript[x, i].*)

Return[(b - a)/2 \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$n$$]$$\*SubscriptBox[\(weights$$, $$i, 2$$]\ f[
\*FractionBox[$$a + b$$, $$2$$] +
\*FractionBox[$$b - a$$, $$2$$]
\*SubscriptBox[$$weights$$, $$i, 1$$]]\)\)]];

f[x_] := Sin[x];
GaussLegendre[0, 1, 8]


I can see two problems : you need to include the function as an argument to GaussLegendreQuadrature and you need to call it with the correct name.

GaussLegendreQuadrature[a_, b_, n_, f_] := Module[{weights, i},
(b-a)/2 * Sum[weights[[i, 2]] f[(a + b)/2 + (b - a)/2 weights[[i, 1]]], {i, 1,n}]]

(* 0.459698 *)

• Compare the result of your routine with N[Integrate[Sin[x], {x, 0, 1}], 20]. – J. M.'s discontentment Jun 18 '12 at 7:31
• @J.M. Just trying to make the code work; didn't mean to check whether it does the right thing. And I missed the factor (b-a)/2. – b.gates.you.know.what Jun 18 '12 at 7:33
• for one thing, b.gatessucks missed the (b - a)/2*. Still, it is a good result – newprint Jun 18 '12 at 7:33
• @user465292 Yes, sorry about that, fixed now. – b.gates.you.know.what Jun 18 '12 at 7:36
• There's no strict need to pass the function as argument. Failing to do so makes a very bad interface, but nevertheless works. Mathematica will just use the global version it finds at evaluation time, which is the one defined directly before the call, i.e. the one @user465292 intended to be called. And as far as I can see he did use the correct name (see at the end of the line containing the first \*SubScriptBox. And yes, I looked in the edit history, and that part doesn't seem to have been changed). – celtschk Jun 18 '12 at 8:35

There's absolutely no need to load a package if all you want to do is simple Gauss-Legendre quadrature:

GaussLegendreQuadrature[f_, {x_, a_, b_}, n_Integer: 10, prec_: MachinePrecision] :=
Module[{nodes, weights},
{nodes, weights} = Most[NIntegrateGaussRuleData[n, prec]];
(b - a) weights.Map[Function[x, f], Rescale[nodes, {0, 1}, {a, b}]]]

GaussLegendreQuadrature[Sin[x], {x, -1, 3}, 10, 20]
1.530294802468585173

N[Integrate[Sin[x], {x, -1, 3}], 20]
1.5302948024685851747


OP does not consider the code to be sufficiently self-evident, so here's a few notes on the implementation:

• Most[NIntegrateGaussRuleData[n, prec]] is intended to give the Gauss-Legendre weights for the integral $\int_0^1 f(x)\,\mathrm dx$; this is equivalent to GaussianQuadratureWeights[n, 0, 1, prec].
• Rescale[x, {0, 1}, {a, b}] == a + (b - a) x; figuring what this substitution has to do with the integral $\int_0^1 f(x)\,\mathrm dx$ is left as an exercise.
• Dot products are an excellent way to do sums of the form $\sum_i w_i f_i$.

Since OP has already accepted the other answer, for the purposes of making this post more useful and convenient than it already is, I will just include advanced methods for generating the Gauss-Legendre nodes and weights had NIntegrateGaussRuleData[] or GaussianQuadratureWeights[] not been available. I'll be linking to the papers explaining these methods, so I won't say more about them.

Here's the Golub-Welsch algorithm:

{nodes, weights} = MapAt[Map[First, #]^2 &, Eigensystem[N[
SparseArray[{{k_, k_} -> 1, {j_, k_} /; Abs[j - k] == 1 :>
With[{i = Min[j, k]}, i/Sqrt[4 i^2 - 1]]}, {n, n}]/2,
prec]], 2]


Here's Laurie's algorithm:

{nodes, weights} = {Diagonal[#2]^2/2, First[#1]^2} & @@
Take[SingularValueDecomposition[
N[SparseArray[{{k_, k_} :> k/Sqrt[k (2 k - 1)],
{j_, k_} /; j == k + 1 :> k/Sqrt[k (2 k + 1)]}, {n, n}],
prec], Tolerance -> 0], 2]


The two methods are intimately related, in the sense that the eigenvalues of a symmetric positive definite matrix are related to the singular values of that matrix's Cholesky triangle. Here, the underlying symmetric positive definite matrix (the Jacobi matrix) is constructed such that its characteristic polynomial is a suitably normalized Legendre polynomial, whose roots are needed for Gaussian quadrature. Again, see the linked papers for more details.

• I have too admit, it look very cool, and probably works way better than my version, but way out of my league. Though, I will browse through docs to find out out what code does. Thanks ! – newprint Jun 18 '12 at 7:56
• I'll write a few explanatory notes, then. – J. M.'s discontentment Jun 18 '12 at 7:58