Example of a Package in a Notebook: Trapezoid Integration of Data Points

Integrate by the Trapezoid Rule

This example of writing your own numerical integration function is from Bahder's excellent book, Mathematica for Scientists and Engineers. I use it to show how to write a Package in a Notebook, as opposed to in a .m file.

While it always pays to see if Mathematica has a built-in function to do what you want, I don't see that NIntegrate can integrate over data points, although it has this nice introduction for the TrapezoidRule: "The trapezoidal rule for integral estimation is one of the simplest and oldest rules (possibly used by the Babylonians and certainly by the ancient Greek mathematicians)."

I needed a function to integrate over data points from an irregular, piecewise, sawtooth-type function. Bahder notes that integrating via the trapezoid rule approximates the integral over a list of (x,y) data points that may be irregularly spaced. Here is the formula. Keyboard shortcuts are faster than a palette; here they are Ctrl+- for subscript and Ctrl+spacebar to exit a subscript. Remember to use a double equal sign for 'equals'.

Here is the numerical integral in a Package. The Package is in one Cell and I set the Cell to Initialization (right-click on the Cell and select Initialization Cell) to automatically load the Package when I evaluate any function in the Notebook.

Note Bahder's use of a named Pattern and Repeated to type-check the data argument. It compares a List of repeated two-element Lists to the pairs that we give it. He then makes sure that there is more than one data pair using Condition. Some say the Wolfram Language is weakly-typed. That is incorrect; it has powerful type-creation and type-checking capabilities.

BeginPackage@"BahderMSE`";

IntegrateTrapezoid::usage="IntegrateTrapezoid[data:{{_,_}..}/;Length@data>1] takes a list of (x,y) pairs, x1 < x2 <x3 ..., that are not necessarily evenly spaced and returns a numerical approximation to their integral using the trapezoid rule.";

Begin@"`Private`";

IntegrateTrapezoid[data:{{_,_}..}/;Length@data>1]:=1/2 Sum[(data[[n+1,1]]-data[[n,1]])(data[[n,2]]+data[[n+1,2]]),{n,Length[data]-1}]

End[];
EndPackage[]

Let's check the usage statement.

In[7]:= ?IntegrateTrapezoid

IntegrateTrapezoid[data:{{_,_}..}/;Length@data>1] takes a list of (x,y) pairs, x1 < x2 <x3 ..., that are not necessarily evenly spaced and returns a numerical approximation to their integral using the trapezoid rule.

Now test it on the same data Bahder used...

In[8]:= data1={{1.2,51.3},{2.3,72.8},{3.0,87.1},{3.7,53.2}};

And we do get the same result.

In[9]:= IntegrateTrapezoid@data1

Out[9]= 173.325

Let's try a different example using Mathematica's SawtoothWave function. This is a typical use of Thread to combine the x and y points in their own Lists.

In[10]:= Thread@{Range[0,3,0.1],SawtoothWave@#&/@Range[0,3,0.1]}

Out[10]= {{0.,0.},{0.1,0.1},{0.2,0.2},{0.3,0.3},{0.4,0.4},{0.5,0.5},{0.6,0.6},{0.7,0.7},{0.8,0.8},{0.9,0.9},{1.,0.},{1.1,0.1},{1.2,0.2},{1.3,0.3},{1.4,0.4},{1.5,0.5},{1.6,0.6},{1.7,0.7},{1.8,0.8},{1.9,0.9},{2.,0.},{2.1,0.1},{2.2,0.2},{2.3,0.3},{2.4,0.4},{2.5,0.5},{2.6,0.6},{2.7,0.7},{2.8,0.8},{2.9,0.9},{3.,0.}}

We see that it is a numerical approximation of the full SawtoothWave where the peaks fall short of 1.

There are three triangles of unit length and height = 0.9, so the total area by a = 1/2 b * h is 1.35. The function works for this example but would need to be validated and its accuracy evaluated for others.

IntegrateTrapezoid@%%

Out[12]= 1.35

Thomas Bahder, Mathematica for Scientists and Engineers

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.