# 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]

-
Welcome to Mathematica.SE! Is it possible to come up with a name that would distinguish yourself somewhat from all the other userxxxxxx's floating around? That would make it feel more like a happy family instead of a prison camp. – Sjoerd C. de Vries Jun 18 '12 at 8:11

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. 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.gatessucks 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.gatessucks 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. Jun 18 '12 at 7:58