Where did I go wrong with my implementation of the trapezoidal rule?
Perhaps you want something like this:
f[x_] := x Log[x]
N[area[1, 2, 100, f]]
(*
-> 0.6363
*)
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 :)
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.