9
$\begingroup$

Is it possible to compute trapezoidal rule numerical integration? I know that Mathematica has Interpolation, and that a list of points can be interpolated and then integrated simply using Integrate. However, my functions are highly oscillatory (they are based on simulation data), and I am not convinced that the interpolation is perfect, even when I set WorkingPrecision to a very high value. Also, I know that ListIntegrate is deprecated, and even if I use it, I am not certain if it is using the trapezoidal rule, which I would like to use.

Do you know if any resources where I can find Mathematica or pseudocode for trapezoidal integration of lists of points? Or do you have any suggestions about how I can use Mathematica efficiently to program such an algorithm myself?

Thanks!

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Why not start at the obvious place: en.wikipedia.org/wiki/Trapezoidal_rule $\endgroup$
    – rm -rf ♦
    Commented May 16, 2012 at 15:51
  • 2
    $\begingroup$ Actually NIntegrate[] indeed is able to perform the trapezoidal rule (see docs for details). I suspect that it won't be the best method for your problem; why not elaborate a bit more on these oscillatory functions you speak of? $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 16, 2012 at 15:53
  • 2
    $\begingroup$ Might check Documentation Center > Integrationtutorial/NIntegrateIntegrationStrategies#144042466 and tutorial/NIntegrateIntegrationRules#32844337 for some ideas on option setting for NIntegrate of finite region oscillatory functions. $\endgroup$
    – Daniel Lichtblau
    Commented May 16, 2012 at 17:04

3 Answers 3

Reset to default
20
$\begingroup$ 
t = Table[{x, Sin[x]}, {x, 0, Pi, .01}];
1/2 Total[((#[[2, 1]] - #[[1, 1]]) (#[[2, 2]] + #[[1, 2]])) & /@  Partition[t, 2, 1]]
(*
-> 1.99998
*)

Perhaps better

1/2 Total[Differences[t[[All, 1]]] ListCorrelate[{1, 1}, t[[All, 2]]]]

They are just

$$\int_a^b f(x)\,dx\approx\frac12\sum_{k=1}^N (x_{k+1}-x_k)(f(x_{k+1})+f(x_k))$$

Edit

Just for fun, using JM's shorter expression:

Manipulate[
   Column[{
     Show[Plot[Sin[x], {x, 0, Pi}], ListLinePlot[#, Filling -> Axis],
         AspectRatio -> Automatic, ImageSize -> 400], 
     Row[{"Approx Integral = ",N@Differences[#1].MovingAverage[#2, 2]& @@ Transpose[#]}]}]&@
     Table[{x, Sin[x]}, {x, 0, Pi, Pi/a}],
 {{a, 2, Dynamic[a]}, 2, 10, 1}]

enter image description here

Improve this answer
$\endgroup$
4
  • 3
    $\begingroup$ Another possibility: Differences[#1].MovingAverage[#2, 2] & @@ Transpose[t]. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 16, 2012 at 16:16
  • $\begingroup$ In your example using ListCorrelate I'd use Differences[...].ListCorrelate[...]/2 instead of Total[Differences[...] ListCorrelate[...]]/2 as the intent is clearer. $\endgroup$
    – rcollyer
    Commented May 17, 2012 at 13:15
  • $\begingroup$ @rcollyer I was just showing off my knowledge of commutativity :). Anyway, I think J. M.'s MovingAverage[] is better, so I used it in the Manipulate example. And yes, Dot[] is better than Total[] here. $\endgroup$
    – Dr. belisarius
    Commented May 17, 2012 at 13:23
  • $\begingroup$ MovingAverage is definitely superior. $\endgroup$
    – rcollyer
    Commented May 17, 2012 at 13:30
3
$\begingroup$

The following is quite fast on packed arrays: 

(Rest@Last@# + Most@Last@#).(Rest@First@# - Most@First@#)/2 &@
 Transpose[t]

Small example:

t = Developer`ToPackedArray@
   Table[{x, Sin[x]}, {x, Subdivide[0., Pi, 100]}];

Differences[#1].MovingAverage[#2, 2] & @@ Transpose[t] //
 RepeatedTiming
1/2 Total[Differences[t[[All, 1]]] ListCorrelate[{1, 1}, t[[All, 2]]]] //
 RepeatedTiming
(Rest@Last@# + Most@Last@#).(Rest@First@# - Most@First@#)/2 &@ Transpose[t] //
 RepeatedTiming
(*
  {0.000036, 1.99984}
  {0.000014, 1.99984}
  {8.5*10^-6, 1.99984}
*)

Larger example (autocompiled by Table, which produces a packed array):

t = Table[{x, Sin[x]}, {x, Subdivide[0., Pi, 10000]}];

Differences[#1].MovingAverage[#2, 2] & @@ Transpose[t] //
 RepeatedTiming
1/2 Total[Differences[t[[All, 1]]] ListCorrelate[{1, 1}, t[[All, 2]]]] //
 RepeatedTiming
(Rest@Last@# + Most@Last@#).(Rest@First@# - Most@First@#)/2 &@ Transpose[t] //
 RepeatedTiming
(*
  {0.0004, 2.}
  {0.00050, 2.}
  {0.000093, 2.}
*)
Improve this answer
$\endgroup$
2
  • $\begingroup$ Super. Your Method "Rest@Last@ ..." is very fast. Would you kindly show use how to use this for an vector of functions like{Sin[x],Cos[x],Tan[x]} . $\endgroup$
    – Gummala Navneeth
    Commented Dec 28, 2020 at 15:38
  • 1
    $\begingroup$ @GummalaNavneeth Perhaps xx = Subdivide[0., Pi/4., 100]; vv = {Sin[x], Cos[x], Tan[x]} /. x -> xx // Developer`ToPackedArray; (vv[[All, 2 ;;]] + vv[[All, ;; -2]]) . (Rest@xx - Most@xx)/2, if you get to generate the data. Advantage of converting to a packed array diminishes as the length of xx approaches 10^6, at which point it becomes a disadvantage. $\endgroup$
    – Michael E2
    Commented Dec 28, 2020 at 16:09
1
$\begingroup$

The Wolfram Demonstrations Project has this demo. It might provide an idea or two to help.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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