This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 1

Consider the following approximate quadrature formula [✶]:

(∫_a^b)f(x)dx≈(∑_(k=0)^o)((Subscript[α, k]^((o)))(((∂^k)f(x))/((∂^k)x))(Subscript[|, x=a])+(Subscript[β, k]^((o)))(((∂^k)f(x))/((∂^k)x))(Subscript[|, x=((a+b)/2)])+(Subscript[γ, k]^((o)))(((∂^k)f(x))/((∂^k)x))(Subscript[|, x=b])).

Calculate the (Subscript[α, k]^((o))), (Subscript[β, k]^((o))), and (Subscript[γ, k]^((o))) for 0≤o≤10 by requiring exactness of the above approximation for polynomials of the largest possible degree. What are the errors of the resulting formulas for the integral (∫_0^(π/2))sin(x)dx?

Before trying to calculate the coefficients, we observe two obvious simplifications. From the symmetry of the integral follows (Subscript[α, k]^((o)))=(((-1))^k)(Subscript[γ, k]^((o))), and from inserting units of length follows that the (Subscript[α, k]^((o))), (Subscript[β, k]^((o))), and (Subscript[γ, k]^((o))) are proportional (((b-a))^(k+1)). The last observation is crucial in avoiding large calculations using nonnumerical quantities and reduces all calculations to solving a linear system of equations with rational coefficients.

makeQuadratureFormula[o_, (* initial degree *) d0_:1, 𝒻_:F] := 
Module[{d = d0, f, lhs, rhs, eqs, sol, solO},
While[(* quadrature formula is exact for polynomials *)
 f[x_] := Sum[c[j] x^j, {j, 0, d}];
 lhs = Integrate[f[x], {x, a, b}];
 (* the right-hand side of the quadrature formula *)
 rhs = Sum[(b - a)^(k + 1) (α[k] Derivative[k][𝒻][a] + 
                            β[k] Derivative[k][𝒻][(a + b)/2] +
                            (-1)^k α[k] Derivative[k][𝒻][b]), {k, 0, o}];
 (* solve linear system *)
 vars = Flatten[Table[{α[j], β[j]}, {j, 0, o}]];
 eqs = Flatten[CoefficientList[#, {a, b}]& /@ 
               DeleteCases[Flatten[CoefficientList[lhs - (rhs /. F -> f), 
                                   Table[c[k], {k, 0, o}]]], 0]];
 sol = Factor //@ Together //@ Solve[eqs == 0, vars]; 
 sol =!= {}, (* increase polynomial degree *) solO = sol; d++]; 
 (* format result *) Factor /@ (rhs /. solO[[1]])]

SetOptions[RowReduce, Method -> OneStepRowReduction]; 
Off[Solve::svars];

Some experimentation with makeQuadratureFormula shows that the maximal polynomial degree is 3 (o+1)-o mod 2. Taking this into account allows avoiding the looping over the polynomial degree. makeQuadratureFormula1 is the resulting function.

makeQuadratureFormula1[o_, 𝒻_:F] := 
makeQuadratureFormula[o, 3 (o + 1) - Mod[o, 2], 𝒻]

Here are the first three resulting quadrature formulas.

makeQuadratureFormula1[0]

makeQuadratureFormula1[1]

makeQuadratureFormula1[2]

Here is the resulting quadrature formula for o=10.

makeQuadratureFormula1[10]

For the integral (∫_0^(π/2))sin(x)dx, the error of the resulting formulas is approximately (10^(-4.5o)).

Table[{o, 1`50 - makeQuadratureFormula1[o] /. 
          {a -> 0, b -> Pi/2, F -> Sin}}, {o, 0, 10}] // N

For generalizations of quadrature rules that contain derivatives, see [✶].