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:

How to Find Memory Used in Computations

At some point you will need to — or just want to — push the limits of Mathematica and it may be helpful to measure memory consumed during a computation. Here are methods inspired by Wagon, Mathematica in Action (2nd ed., pp 431 et seq.). (See also Memory Management Tools.)

The first method is simple; store the initial MemoryInUse in a variable, do your calc, then subtract it from the current MemoryInUse.

currentMemory=MemoryInUse[];
{(MemoryInUse[]-currentMemory),Plus@@(1./Range@1000000)}

{1368,14.3927}

Second, suppose you want to monitor memory in several places in your code. Initialize the variable currentMemory, then sprinkle

Sow[-currentMemory+(currentMemory=MemoryInUse[])]

through your code like Print statements, and the enclosing Reap will cause each one to output the net memory consumed in bytes between the previous use and the current one. I say 'net' memory consumed because Mathematica is designed to optimize memory use and so clears memory in use every chance it gets. Following the pattern of Timing, I reverse the output to give the memory results first, then the function's output.

currentMemory=MemoryInUse[];
(*initialize s*); s=0;
Reap[Do[If[i<25,Sow[-currentMemory+(currentMemory=MemoryInUse[])]];s+=1./i,{i,1000}];s]//Reverse

{{{3016,11216,16,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}},7.48547}

Did you notice, as Wagon points out, that Do uses minimal memory during a calculation? It cleans up after every iteration.

Like Timing, this third method of checkMemory requires you to wrap your computation in checkMemory, access the result of your computation with First and if desired suppress the output with a semi-colon. Mimicking the structure of a built-in function is an example of a good programming technique. The designers of Mathematica are master programmers and they have created and debugged templates or paradigms of function structures. Often you can follow the structure of built-in functions to create sound ones of your own.

Clear@checkMemory;
checkMemory@computation_:=
Module[{computation1,preComputationMemory},
preComputationMemory=MemoryInUse[];Print["Memory in use at start is ",preComputationMemory];
computation1=computation;
memoryUsedInComputation=MemoryInUse[]-preComputationMemory;
Print["Memory in use at end is ",MemoryInUse[]];{memoryUsedInComputation,computation1}
]

Once it's done, this computation uses 32 bytes of memory and the result is 14.3927.

checkMemory[Sum[1./i,{i,1000000}]]

Memory in use at start is 31286024
Memory in use at end is 31286056
{32,14.3927}

We find the memory use with the function output for a small Table.

s1=0;checkMemory[Table[s1+=1./i,{i,10}]]

Memory in use at start is 31289472
Memory in use at end is 31289752
{280,{1.,1.5,1.83333,2.08333,2.28333,2.45,2.59286,2.71786,2.82897,2.92897}}

Then suppress output, with a semi-colon as is done in Timing, for a large Table.

checkMemory[s1=Table[1./i,{i,10^7}];]

Memory in use at start is 31340688
Memory in use at end is 111344696
{80003992,Null}

I left off the summation of the Table to show MemoryInUse during the calc. If we add Total, we see the max MemoryInUse reclaimed (-80,000,024 bytes).

Timing[checkMemory[s1=Table[1./i,{i,10^7}]//Total]]

Memory in use at start is 31359736
Memory in use at end is 31359864
{0.343202,{128,16.6953}}

Compare memory used in a Do loop. Pretty interesting!

Timing[s1=0;checkMemory[Do[s1+=1./i,{i,10^7}];s1]]

Memory in use at start is 31364400
Memory in use at end is 31364432
{10.296066,{32,16.6953}}

First@%/First@%%

30.0000

So there you have it, Do takes 30 times as long as the functional command Table[...]//Total, but uses a minute fraction of the latter's memory.

Examples of Different Programming Styles with Timing

This comparison of different programming styles is adapted from Roman Maeder's in his Programming in Mathematica course. Here is a List of expressions to square, including different types on numbers and an undefined symbol.

alist={a,1.1,2+3I,4,573297329847};

Functional Style

Here is the shortest solution, which works due to the function Attribute, Listable, of Power that maps Power over a List. You can add Listable to your own functions to achieve the same simplicity.

squareListListable@list_List:=list^2

squareListListable@alist
{a^2,1.21,-5+12 I,16,328669828409699917043409}

As Maeder says, it really cannot get much simpler than that. Here are two more functional-style solution. These are also compact. The first uses Table, a very powerful function that leads us to not deconstruct lists (arrays) but to think of applying a function to them as a whole. The simple squaring of each List item could be replaced with a function of arbitrary complexity.

squareFunctional@list_List:=Table[i^2,{i,list}]

Beginners are often unaware of Table's ability to iterate directly over List items, as above, without this unnecessary clutter. However, see the Timing analysis below for a surprise.

squareFunctionalWithIterator@list_List:=Table[list[[i]]^2,{i,Length@list}]

This third functional example uses Map, another powerful function that, along with Listable and Table, replaces the Do loop. A pure Function expresses the squaring operation concisely.

squareListMapped@list_List:=Map[#^2&,list]

Rule-Based Style

Programming with replacement rules is the way Mathematica itself works — Every Change is a Transformation — and rule-based functions are often very concise and easy to understand. I prefer rule-based functions over other styles.

squareRuleBasedShorter@list_List:=list/.{first___,x_,rest___}->{first,x^2,rest}

Recursive Style

Here is a recursive rule-based approach — the function calls itself and keeps marching down a 'rest' of List that keeps getting smaller. Note that flattening deeply nested Lists is a very fast operation (Log@n) whose timing you rarely need to worry about. However this recursive approach is only practical for small examples.

squareRuleRecursive@list_List:=list/.{x_,rest___}:>{x^2,squareRuleRecursive@{rest}}//Flatten

Procedural Style

While I really think which style to use is a personal decision ('de gustibus non est disputandum'), this procedural solution does seem cumbersome and old-fashioned compared to the functional and rule-based ones. But sometimes procedural functions are easier and faster to write than laboring to find more concise functional or rule-based ones.

squareProcedural@list_List:=Module[{array1={},squaresList,
listLength=Length@list},
Do[array1=Prepend[array1,0],{i,listLength}];
Do[array1[[i]]=list[[i]]^2,{i,listLength}];
Return@array1
]

Maeder cheats a bit, using Table to initialize his array, which speeds up his function considerably as you'll see in the Timing analysis.

squareProceduralMaeder@list_List:=Module[
{squaresArray,
listLength=Length@list},
(*initialize the array*)
squaresArray=Table[0,{listLength}];
Do[squaresArray[[i]]=list[[i]]^2,{i,listLength}];
Return@squaresArray
]

Timing Analysis

aLongList=Range[10^6];

The Listable function is the clear winner.

Timing[squareListListable@aLongList;]
{0.,Null}

Let's push it harder to see more of its speed advantage. It appears to be 10 times faster than its competitors below.

Timing[squareListListable@Range[10^7];]
{0.015600,Null}

Timing[squareFunctional@aLongList;]
{0.031200,Null}

Here's the surprise I mentioned. Not sure why, but using the iterator is exactly twice as fast than 'directly' iterating over List items. I'd only worry about this for functions needing optimization though.

Timing[squareFunctionalWithIterator@aLongList;]
{0.015600,Null}

Timing[squareListMapped@aLongList;]
{0.015600,Null}

It surprise me that the rule-based approach is relatively slow.

Timing[squareRuleBased@aLongList;]
{0.062400,Null}

Both the recursive and procedural approaches are very slow.

Timing[squareRuleRecursive@aLongList;]
During evaluation of In[186]:= $RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>
Out[186]= {9.188459,Null}

Timing[squareProcedural@aLongList;]
$Aborted

Timing[squareProceduralMaeder@aLongList;]

{2.012413,Null}

An Advanced Predicate and Analysis of How It Works with Trace and Timing

Wagner gives a different version of allodd@x that is more efficient than the one shown here. This function stops testing as soon as it finds a number that is not odd. Contrary to his predisposition, Wagner uses a procedural loop, "de-constructing," as he calls it, the List to operate on its Parts.

In[269]:= allOddProcedural@aList_List:=Module[{i},
For[i=1,i<=Length@aList,i++,If[!OddQ[aList[[i]]],Break[]]];i>Length@aList]

But I bet OddQ is Listable and therefore its C compilation would be much faster than the For loop. First let's compare the Wagner allOddQ from above with his procedural one.

Wagner's function Breaks after OddQ returns False when testing the third argument, 6 (near the end of the Trace). Note the 100-fold difference in Timing!

In[275]:= list2={3,5,6,7};

In[276]:= allOdd@list2//Trace

Out[276]= {{list2,{3,5,6,7}},allOdd[{3,5,6,7}],Length[Select[{3,5,6,7},OddQ]]==Length[{3,5,6,7}],{{Select[{3,5,6,7},OddQ],{OddQ[3],True},{OddQ[5],True},{OddQ[6],False},{OddQ[7],True},{3,5,7}},Length[{3,5,7}],3},{Length[{3,5,6,7}],4},3==4,False}

In[277]:= allOddProcedural@list2//Trace

Out[277]= {{list2,{3,5,6,7}},allOddProcedural[{3,5,6,7}],Module[{i$},For[i$=1,i$<=Length[{3,5,6,7}],i$++,If[!OddQ[{3,5,6,7}[[i$]]],Break[]]];i$>Length[{3,5,6,7}]],{For[i$41099=1,i$41099<=Length[{3,5,6,7}],i$41099++,If[!OddQ[{3,5,6,7}[[i$41099]]],Break[]]];i$41099>Length[{3,5,6,7}],{For[i$41099=1,i$41099<=Length[{3,5,6,7}],i$41099++,If[!OddQ[{3,5,6,7}[[i$41099]]],Break[]]],{i$41099=1,1},{{i$41099,1},{Length[{3,5,6,7}],4},1<=4,True},{{{{{i$41099,1},{3,5,6,7}[[1]],3},OddQ[3],True},!True,False},If[False,Break[]],Null},{i$41099++,{i$41099,1},{i$41099=2,2},1},{{i$41099,2},{Length[{3,5,6,7}],4},2<=4,True},{{{{{i$41099,2},{3,5,6,7}[[2]],5},OddQ[5],True},!True,False},If[False,Break[]],Null},{i$41099++,{i$41099,2},{i$41099=3,3},2},{{i$41099,3},{Length[{3,5,6,7}],4},3<=4,True},{{{{{i$41099,3},{3,5,6,7}[[3]],6},OddQ[6],False},!False,True},If[True,Break[]],Break[]}},{{i$41099,3},{Length[{3,5,6,7}],4},3>4,False},False},False}

In[262]:= list1=Range[1,10^7,2];

Timing[#@list1]&/@{allOddProcedural,allodd}

{{9.297660,True},{0.109201,<<1>>}}

Why Use Built-in Functions? They're Already Written, De-bugged, and Fast

Now let's try my premise - first, OddQ is Listable as I thought.

In[260]:= Attributes@OddQ

Out[260]= {Listable,Protected}

Let's write the simple function:

In[279]:= allOddFunctional@aList_List:=If[!OddQ@aList,i>Length@aList]

In[280]:= Timing[allOddFunctional@list1]

Out[280]= {0.093601,<<1>>}

allOddFunctional is 10X faster than Wagner's procedural function, even though it appears to check every argument:

In[281]:= allOddFunctional@list2//Trace

Out[281]= {{list2,{3,5,6,7}},allOddFunctional[{3,5,6,7}],If[!OddQ[{3,5,6,7}],i>Length[{3,5,6,7}]],{{OddQ[{3,5,6,7}],{OddQ[3],OddQ[5],OddQ[6],OddQ[7]},{OddQ[3],True},{OddQ[5],True},{OddQ[6],False},{OddQ[7],True},{True,True,False,True}},!{True,True,False,True}},If[!{True,True,False,True},i>Length[{3,5,6,7}]]}

The lesson is simply that built-in functions are almost always faster than ones we make if for no other reason than they are compiled into C.

Here Wagner shows the use of While instead of For to stop when it finds an odd number.

allodd@aList_List:=Module[{i=1},While[i<=Length@aList&&OddQ@aList[[i]],i++];i>Length@aList]

Here is a Dropbox link to Wagner's book. High kudos to "Mr Wizard Todd" for asking McGraw-Hill for the right to distribute Wagner's out-of-print book. Here's his post on the StackExchange Mathematica forum.

https://www.dropbox.com/s/kllwg6y44p8va5g/Wagner%20All%20Parts-RC.pdf

Here is a link to a compressed file: http://www.verbeia.com/mathematica/PowerProgMa.zip