Unexpected result of summation

Question

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].*)

       weights = GaussianQuadratureWeights[n, -1, 1];
       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]

Answer 1

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},
  weights = GaussianQuadratureWeights[n, -1, 1];
  (b-a)/2 * Sum[weights[[i, 2]] f[(a + b)/2 + (b - a)/2 weights[[i, 1]]], {i, 1,n}]]    

GaussLegendreQuadrature[0, 1, 8, Sin]

(* 0.459698 *)

Answer 2

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[NIntegrate`GaussRuleData[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[NIntegrate`GaussRuleData[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 NIntegrate`GaussRuleData[] 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.