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:

Recursive Programming: The General Principle

Recursion, where a function calls itself, is a fundamental programming method. Recursion is a different way of iterating over a List than using a loop. Recurse is from the Latin recurrere, to run back (I picture a dog retrieving a stick over and over).

There are two principles of writing a recursive function, which is piecewise:

1. Recursive case: Define a function f(n+1) recursively in terms of f(n), or equivalently, f(n) in terms of f(n-1)
2. Base case: Define one or more base cases, which are stopping points for the recursion, so that it doesn't become infinite. 

An elementary example is exponentiation defined this way:

xⁿ = x•x^(n-1)

First we define the base case, so that Mathematica will match it to prevent an infinite recursive loop. Let's just consider positive integers and zero. The base case is any integer raised to the 0^th power equals 1.

exponentiation[x_,0]:=1;

The recursive case is:

exponentiation[x_,n_Integer/;x>0]:=x exponentiation[x,n-1]

Start testing a recursive function with the simplest examples:

exponentiation[0,0]
1

exponentiation[1,0]
1

exponentiation[1,1]
1

Let's make sure the Condition requiring x > 0 works; it does.

exponentiation[1,-1]
exponentiation[1,-1]

Trace shows the exponentiation hitting the base case and stopping.

exponentiation[1,1]//Trace
{exponentiation[1,1],1 exponentiation[1,1-1],{{1-1,0},exponentiation[1,0],1},1 1,1}

To compute 2², the function calls itself to computer 2^2-1 i.e. 2¹, which in turn calls itself to compute 2⁰, which is the base case, 1. Once it has recursed down to the base case, it has a value to substitute in the preceding computations and the recursion 'unwinds.' All recursive functions work this way.

exponentiation[2,2]//Trace
{exponentiation[2,2],2 exponentiation[2,2-1],{{2-1,1},exponentiation[2,1],2 exponentiation[2,1-1],{{1-1,0},exponentiation[2,0],1},2 1,2},2 2,4}

exponentiation[2,10]
1024

exponentiation[2,100]
1267650600228229401496703205376

exponentiation[2,1000]
10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376

Use Block for Large Recursions and Debugging

You may run into the built-in default recursion limit, $RecursionLimit, of 1024 when performing a large recursion.

In[584]:= exponentiation[2,10000]
During evaluation of In[584]:= $RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>
Out[584]= Hold[exponentiation[2,8978-1]]

The safe way to exceed the limit is using Block. Block temporarily changes the value of a Global or System variable. I've suppressed the output here since it is 3,011 digits long, but you can try this safely at home without "//IntegerDigits//Length".

Block[{$RecursionLimit=10050},exponentiation[2,10000]]//IntegerDigits//Length
3011

Mathematica will warn you and stops an infinite recursion with a default recursion limit specified by $RecursionLimit ($IterationLimit is used similarly to stop infinite loops).

{$RecursionLimit,$IterationLimit}
{1024,4096}

Just as you can use Block to temporarily exceed $RecursionLimit, you can use Block to temporarily 'tighten up' $RecursionLimit while writing a recursive function.

In[592]:= Block[{$RecursionLimit=20},factorial@x_:=x*factorial[x-1];factorial@2]
During evaluation of In[592]:= $RecursionLimit::reclim: Recursion depth of 20 exceeded. >>
Out[592]= Hold[factorial[-14-1]]

Always remember to write the base case rule!

In[593]:= Block[{$RecursionLimit=20},factorial@x_:=x*factorial[x-1];factorial@0=1;factorial@2]
Out[593]= 2

Now we may safely compute the factorial of larger numbers:

factorial@200

788657867364790503552363213932185062295135977687173263294742533244359449963403342920304284011984623904177212138919638830257642790242637105061926624952829931113462857270763317237396988943922445621451664240254033291864131227428294853277524242407573903240321257405579568660226031904170324062351700858796178922222789623703897374720000000000000000000000000000000000000000000000000

Note all the zeroes at the end. They must be from each time a multiple of 10 is incorporated into the factorial.