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

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?

-
area[1, 2, 100, # Log@# &] –  belisarius Oct 21 '13 at 14:35
@belisarius,+1 Dear belisarius,how to rewrite to avoid using the pure fuction.Thanks sincerely! –  Shutao Tang Oct 21 '13 at 14:40
If you want to use the function as a parameter ... it should be pure function :) –  belisarius Oct 21 '13 at 14:46
f[x_]:=x Log @ x. Then pass f in the function. –  Sjoerd C. de Vries Oct 22 '13 at 5:39

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 :)

-
,Dear m_goldberg,thanks for your detailed solution.However,area[# Log@# &, 1., 2., 100] gives the result of area[#1 Log[#1] &, 1., 2., 100] –  Shutao Tang Oct 22 '13 at 4:32
plot[f_, a_, b_, n_] := Show[{ListPlot[#, Filling -> Axis], Graphics[{Line@#, Red, Point@MovingAverage[#, 2]}]} &@ Transpose@{#, f /@ #} &@Range[a, b, (b - a)/n]] –  belisarius Oct 22 '13 at 5:01
Still shorter: plot[f_, a_, b_, n_] := ListPlot[#, Filling -> Axis, Epilog -> {Line@#, Red, Point@MovingAverage[#, 2]}] &@ Transpose@{#, f /@ #} &@Range[a, b, (b - a)/n] –  belisarius Oct 22 '13 at 13:35
@Tangshutao. There was a syntax error in the argument pattern that I gave for area. I have edited my answer to correct this error. Thanks for pointing out the problem. –  m_goldberg Oct 22 '13 at 15:16

Perhaps you want something like this:

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

-
,Dear belisarius, if I define the function in the module,it cannot give the right result?Can you tell me why?Thanks sincerely:).area1[a0_, b0_, n_, function0_] := Module[{a = a0, b = b0, \[CapitalDelta]x, xStart, result, function}, \[CapitalDelta]x = (b - a)/n; xStart[i_] := a + i \[CapitalDelta]x; function[x_] := function0; result = N[(Sum[2 function[xStart[i]], {i, 1, n - 1}] + function@a + function@b) \[CapitalDelta]x/2]] area1[1, 2, 100, x Log@x] –  Shutao Tang Oct 22 '13 at 4:35
@Tangshutao it's a scoping problem. The pattern x_ inside the Module isn't referencing the same symbol than the x symbol in your x Log@x. To see it add a Print@Information@function in the Module after the function definition and you'll see that Mathematica appends a \$ to the x_ to isolate contexts. –  belisarius Oct 22 '13 at 4:46

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


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

,Dear Peltio, Is it conflicing to define two fuctions such as trapIntegrate[data : {{_, _} ..}, h_] and trapIntegrate[f_, {x_, a_, b_}, n_].Thanks. –  Shutao Tang Oct 23 '13 at 4:01