# 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@# &] Oct 21, 2013 at 14:35
• @belisarius,+1 Dear belisarius,how to rewrite to avoid using the pure fuction.Thanks sincerely!
– xyz
Oct 21, 2013 at 14:40
• If you want to use the function as a parameter ... it should be pure function :) Oct 21, 2013 at 14:46
• f[x_]:=x Log @ x. Then pass f in the function. Oct 22, 2013 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]
– xyz
Oct 22, 2013 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]] Oct 22, 2013 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] Oct 22, 2013 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. Oct 22, 2013 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]
– xyz
Oct 22, 2013 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. Oct 22, 2013 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


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.

• ,Dear Peltio, Is it conflicing to define two fuctions such as trapIntegrate[data : {{_, _} ..}, h_] and trapIntegrate[f_, {x_, a_, b_}, n_].Thanks.
– xyz
Oct 23, 2013 at 4:01
• The first procedure has two arguments, the first one being a list of two-elements lists; the second procedure has three arguments. They could not possibly conflict. According to what you pass to trapIntegrate, the correct procedure is chosen. Granted, I did not implement sophisticated argument checking, but that's the way I code(d back then): I did it for myself and so I knew what did what and why (lean > robust). You might as well be worried about conflict in the first two definitions but then, mma sorts that out by itself: more specific definitions take the precedence over generic ones. Oct 23, 2013 at 9:37
• I forgot to add that after your comment I edited the post to reflect the above considerations. Oct 23, 2013 at 15:19