## 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 fastThe 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:

## No comments:

## Post a Comment