Where did I go wrong with my implementation of the trapezoidal rule?

Question

One method for doing quadrature, called the trapezoidal rule, improves accuracy by connecting the points on the curve corresponding to the points of subdivision with line segments, forming trapezoidal approximations of the area in place of rectangular approximations.

I tried to use this method to caculate the integral of $x \ln(x)$ over $[1,2]$. Here is what I have:

area[a0_, b0_, n_, function0_] := 
  Module[{a = a0, b = b0, Δx, xStart, function = function0, result},
    Δx = (b - a)/n;
    xStart[i_] := a + i Δx;
    result = 
      N[(Sum[2 function[xStart[i]], {i, 1, n - 1}] + function @ a + function @ b) Δx/2]]

area[1, 2, 100, x Log @ x]

Unfortunately, it doesn't give the result that I want. Where have I gone wrong?

Answer 1

Your code is far from optimal. Perhaps a tutorial that presents a better implementation would be helpful. Despite my code being very different from yours, it is still a straight-forward implementation of the textbook formula for the trapezoid rule. Because it takes advantage of just a few of Mathematica's built-in functions, it is much more compact than your implementation.

area[f_, a_?NumericQ, b_?NumericQ, n_Integer?Positive] :=
  With[{dx = (b - a)/n}, Plus @@ MovingAverage[f /@ Range[a, b, dx], 2] dx]

area[# Log @ #&, 1., 2., 100]

0.6363

This simple code works because f /@ Range[a, b, dx] computes $f(x)$ at all the mesh points and the moving average computes $(f(x_i)+f(x_{i+1}))/2$ for every successive pair of mesh points, $x_i$ and $x_{i+1}$. The $i$-th such average, when multiplied by the mesh interval, $dx$, gives the area of the $i$-th trapezoid.

I think it might be useful to illustrate what I said in the previous paragraph by making a plot much like the one shown in the question, but showing where the successive averages lie.

plot[f_, a_, b_, n_] :=
  With[{dx = (b - a)/n},
    Module[{meshPts, trapPts, p1, p2},
      meshPts = {#, f[#]}& /@ Range[a, b, dx];
      trapPts =
        Transpose[{
          Range[a + dx/2, b - dx/2, dx], 
          MovingAverage[f /@ Range[a, b, dx], 2]}];
      p1 = ListPlot[trapPts, PlotStyle -> {Red, PointSize[Medium]}];
      p2 = ListPlot[meshPts, Filling -> Axis, AxesOrigin -> {a, 0.}];
      Show[p1, p2, 
        AxesOrigin -> {a, 0.}, 
        PlotRange -> {{a, b}, {0., f[b]}},
        Prolog -> {Line[meshPts]}]]]
 plot[# ^4 &, 1., 2., 4]

Making the illustration is perversely more difficult than writing the function implementing the trapezoid rule :)

Answer 2

Perhaps you want something like this:

f[x_] := x Log[x]
N[area[1, 2, 100, f]]
(*
 -> 0.6363
*)

Answer 3

Here are a few procedures I wrote on a very old version of Mathematica many years ago. I believe the last one answers your question on how to pass a function not in pure form. It basically amount to using a replacement rule inside the procedure.

The following definitions show how easy it is to overload procedures in Mathematica. The first form of trapIntegrate works with a one-dimensional list of ordinates (y-values only). It is assumed the abscissas are uniformly spaced with step h

trapIntegrate[data_List, h_] := h*(Plus @@ # - (First[#] + Last[#])/2) &[data]

The second form accepts a list of coordinates {x,y} with uniform spacing h. While the step could be inferred from the x-values, it is explicitly passed to the procedure for three reasons: 1) it simplify the coding; 2) it allows to differentiate this procedure call from the following one and 3) it makes coding the last form of the procedure easier.

trapIntegrate[data : {{_, _} ..}, h_] := h*(Plus @@ # - (First[#] + 
Last[#])/2) & [Last /@ data]

The third form accepts a list of coordinates with variable step. Now the step has to be inferred from the x-values, one trapezoid at the time. Of course it si also possible to pass a uniformly spaced data (but this procedure will be slower then the one expressly thought for uniform spaced points).

trapIntegrate[data : {{_, _} ..}] :=
  Module[
    {xvals, yvals, xdiffs, fsums},
    xvals = First /@ data; yvals = Last /@ data;
    xdiffs = Drop[xvals, 1] - Drop[xvals, -1];
    fsums = (Drop[yvals, 1] + Drop[yvals, -1])/2;
    xdiffs.fsums]

The fourth and final form accepts a function of the variable x on the interval [a,b]. Here I specified the number of steps n = (b-a)/h. It simply computes the data to pass the uniform trapIntegrate procedure. I compute the step here, and then pass it onto that procedure.

trapIntegrate[f_, {x_, a_, b_}, n_] :=
  Module[
    {data},
    data = Table[f /. x -> xk, {xk, a, b, (b - a)/n}];
    trapIntegrate[data, (b - a)/n]
    ]

Now, you can compute your integral. This will give an exact value expression in the form of a sum of rational numbers multiplied by logarithms.

trapIntegrate[x Log[x], {x, 1, 2}, 100]

To speed thing up you can either specify machine precision bounds and/or number of steps (to force the conversion from exact numbers to machine precision numbers)

 trapIntegrate[x Log[x], {x, 1., 2.}, 100.]

or create 'approximate numerical' version like trapNIntegrate that convert data to an approximate numerical form with N[].

 trapNIntegrate[f_, {x_, a_, b_}, n_] :=
   Module[
     {data},
     data = N[ Table[f /. x -> xk, {xk, a, b, (b - a)/n}] ];
     trapIntegrate[data, N[(b - a)/n]]
     ]

Please note that

trapIntegrate[x Log[x], {x, 1, 2}, 100] // N

will give you the numerical result you are after, but this will still require computing the exact numerical value first. Forcing N on your data can result in considerable speed improvements. From what I remember these procedure - when supplied with machine precision data - were considerably faster (but dumber) than built-in procedures. This is probably no longer true with newer versions of Mathematica.