How It Works: Trapezoid Integration

How It Works : Integrate by the Trapezoid Rule

In related post I used a numerical integration function from Bahder, Mathematica for Scientists and Engineers, to integrate a sawtooth function. This is a great book. Despite having been written in 1995, it's a wealth of 'under the hood' insights into the Wolfram Language. Here is a second function from Bahder for integration via the trapezoid rule that is more sophisticated than the first one. And due to using built-in Listable functionality in Dot and Plus it's about twice as fast

The syntax of the parameter uses a named Pattern to check the data parameter to see if it is a List of 2-element (x,y) pairs, then uses Condition to make sure that data contains more than one point.

Here is the formula for integration by the trapezoid rule:

And here is Bahder's code. The syntax of the parameter uses a named Pattern to check data to see if it is a List of 2-element (x,y) pairs, then uses Condition to make sure that data contains more than one point.

Clear@integrateTrapezoid;
integrateTrapezoid[data:{{_,_}..}/;Length@data>1]:=Module[{xDifferences,ySums},
xDifferences=Dot[{-1,1},#]&/@Partition[First/@data,2,1];
ySums=Apply[Plus,Partition[Last/@data,2,1],{1}];
Plus@@(xDifferences,ySums)/2
]

The Dot product multiplies each vector (List) term by term and sums them.

Dot[{a,b,c},{x,y,z}]

a x+b y+c z

Let's check to make sure the abbreviated operators for Apply and Divide will work as intended; which they will. Apply is 'stickier' than Divide and so Plus will be evaluated before taking the average (dividing by 2).

Precedence/@{Apply,Divide}

{620.,470.}

We test it on his data and.. it doesn't work. There is a bug.

data1={{1.2,51.3},{2.3,72.8},{3.0,87.1},{3.7,53.2}};

integrateTrapezoid@data1

{62.6,80.3,70.5}

Let's add some Print statements, which suffice for debugging all pure functional programs, that is, short ones.

Clear@integrateTrapezoid;integrateTrapezoid[data:{{_,_}..}/;Length@data>1]:=Module[{xDifferences,ySums},
xDifferences=Dot[{-1,1},#]&/@Partition[First/@data,2,1];
Print@xDifferences;
ySums=Apply[Plus,Partition[Last/@data,2,1],{1}];
Print@ySums;
Plus@@(xDifferences,ySums)/2
]

With the Print statement results it's easy to see that xDifferences is indeed the difference between successive x-points and ySums is the sum of each pair of successive y-values.

integrateTrapezoid@data1

{1.1,0.7,0.7}
{124.1,159.9,140.3}
{62.6,80.3,70.5}

I see the bug. We should be multiplying xDifferences by ySums, not listing them. There's a lesson: Always use a multiplication symbol in code rather than a space, which might be mis-read as I did in Bahder's code. We'll add the multiplication and remove the Print statements.

Clear@integrateTrapezoid;integrateTrapezoid[data:{{_,_}..}/;Length@data>1]:=Module[{xDifferences,ySums},
xDifferences=Dot[{-1,1},#]&/@Partition[First/@data,2,1];
ySums=Apply[Plus,Partition[Last/@data,2,1],{1}];
Print@(xDifferences*ySums);
Plus@@(xDifferences*ySums)/2
]

Now we get the correct answer.

integrateTrapezoid@data1

{136.51,111.93,98.21}
173.325

Let's try it on something else.

data2={{0.0,0.0},{1.0,1.0},{2.0,0.0},{3.0,2.0},{4.0,0.0},{5.0,2.0},{6.0,1.0},{7.0,3.0},{8.0,0.0},{9.0,2.0},{10.0,0.0}};
Graphics[Line@data2,Axes->True]

integrateTrapezoid@data2

{1.,1.,2.,2.,2.,3.,4.,3.,2.,2.}
11.

Trusting we got xDifferences and ySums right, we can copy and paste them into the program final line by hand to double-check its validity.

Plus@@{1.`,1.`,2.`,2.`,2.`,3.`,4.`,3.`,2.`,2.`}/2

11.

Recommended: