Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

enter image description here

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?

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

3 Answers 3

up vote 6 down vote accepted

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]

Per

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

share|improve this answer
    
,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] –  ShutaoTang Oct 22 '13 at 4:32
1  
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
1  
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
1  
@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
*)
share|improve this answer
    
,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] –  ShutaoTang Oct 22 '13 at 4:35
2  
@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]]
     ]

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.

share|improve this answer
1  
,Dear Peltio, Is it conflicing to define two fuctions such as trapIntegrate[data : {{_, _} ..}, h_] and trapIntegrate[f_, {x_, a_, b_}, n_].Thanks. –  ShutaoTang Oct 23 '13 at 4:01
1  
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. –  Peltio Oct 23 '13 at 9:37
    
I forgot to add that after your comment I edited the post to reflect the above considerations. –  Peltio Oct 23 '13 at 15:19

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.