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:

Apply a List of Functions to a Corresponding List of Arguments

Here are two ways to apply a List of functions to a List of arguments. The first method, from Maeder, uses Inner. Note the use of a purely symbolic example to reveal what Inner is doing.

Inner[#@#2&,{f1,f2,f3,f4},{arg1,arg2,arg3,arg4},List]

{f1[arg1],f2[arg2],f3[arg3],f4[arg4]}

In #1@# &, Slot1 (#1 or abbreviated #) pulls an argument from the first List, which are the functions, and Slot2 (#2 where you do need the "2") pulls an argument from the second List, which are the arguments.

If you wanted to apply a second function, g, to the results rather than just listing them, Inner provides for that:

Inner[#@#2&,{f1,f2,f3,f4},{arg1,arg2,arg3,arg4},g]

g[f1[arg1],f2[arg2],f3[arg3],f4[arg4]]

The second method, from my Mathematica teacher Harry Calkins, uses MapThread, which requires wrapping the function List and the argument List together in a List:

(arguments={{a,b,c},{d,e,f},{r,s,t},{u,v,w}});

(functions={3 #&,3+#&,#/2&,Plus@@@Partition[#,2,1,{1,1},#]&});

He gave an example using some built-in functions and structured the functions and arguments in matrices.By enclosing the assignments in parentheses he showed them as matrices without inadvertently making the assignment to the MatrixForm-formatted List.

arguments//MatrixForm
functions//MatrixForm

(results=MapThread[#1@#2&,{functions,arguments}])//MatrixForm

Grid - Partition - Identify Primes with Formatting

Starting with an example from Mangano, Mathematica Cookbook, and the Doc Center, here are examples of using Partition with Grid and formatting to highlight the values of interest.

Grid@Partition[If[PrimeQ@#,Framed[#,FrameStyle->Red],#]&/@Range@200,20]

Grid@Partition[If[PrimeQ@#,Style[#,FontColor->Red],#]&/@Range@200,20]

Grid@Partition[If[PrimeQ@#,Style[#,FontColor->Red,FontWeight->Bold],Style[#,FontColor->Gray]]&/@Range@200,20]

This sets up a clickable square from which you can interactively choose a color.

Clear@highlightColor;ColorSetter@Dynamic@highlightColor

Then the selected color is used here:

Grid@Partition[If[PrimeQ@#,Style[#,FontColor->highlightColor,FontSize->16,FontWeight->Bold],Style[#,FontColor->Gray]]&/@Range@200,20]

Grid@Partition[If[!PrimeQ@#,Style[#,FontColor->White],#]&/@Range@200,20]

Change the Partition from 20 to 10 to accommodate larger numbers.

Grid@Partition[If[!PrimeQ@#,Style[#,FontColor->White],#]&/@Range[500000,500200],10]

How to See the Equivalence of Select and Cases

MatchQ and PatternTest Equate Cases and Select

Select and Cases are essential Mathematica filtering functions. Select filters expressions, usually Lists, with predicates, such as ones that end with Q (evaluate ?*Q to see them) or those that don't, such as Greater, Less, and others that I list here. Cases uses patterns such as testing for a Head (e.g. _Integer) or a more general pattern (e.g. for a 2-element List, {x_,y_}.

Cases and Select are equivalent and a versatile Mathematica programmer knows how to use either one in any circumstance according to the need of the moment. The keys to seeing their equivalence are MatchQ, which translates a pattern into a predicate, and PatternTest (_?test), which translates a predicate into a pattern.

MatchQ[4.1,_Real]

True

MatchQ[{2,2.1},{_Integer,_Real}]

True

Thus you can take any usage of Cases and translate it into Select.

Cases[{a,4,4.2,6,7.5,8},x_Integer/;x>4]

{6,8}

integerGreaterThan4Q@n_:=MatchQ[n,_Integer&&n>4];

Select[{a,4,4.2,6,7.5,8},integerGreaterThan4Q]

{6,8}

And you can take any usage of Select and translate it into Cases.

Select[{1,2,4,7,6,2},EvenQ]

{2,4,6,2}

Cases[{1,2,4,7,6,2},_?EvenQ]

{2,4,6,2}

If we built our own predicate for EvenQ using MatchQ, it works the same way.

Select[{1,2,3,4},MatchQ[#/2,_Integer]&]

{2,4}

However note the need for parentheses around the test, which is often the case with PatternTest - the no parens usage fails.

Cases[{1,2,3,4},_?(MatchQ[#/2,_Integer]&)]

{2,4}

Cases[{1,2,3,4},_?MatchQ[#/2,_Integer]&]

{}

The reason is that Blank and PatternTest both have a higher Precedence ('stickiness') than Function, and so they stick together and PatternTest ends up inside the pure Function instead of using the pure Function as the pattern test. You don't need to worry about that — just as a rule of thumb put parens around any 'homemade' test.

Precedence/@{PatternTest,Blank,Function}

{680.,730.,90.}

The Q Predicates

This is a handy way to use Partition, especially using {} so a final short sublist won't be dropped from a 'ragged' List. The syntax says

Partition[list, sublist length, offset length, start Position, padding]

and here means divide the List into Partitions of length 4, moving the window that same length 4 with each Partition so there is no overlap, start the List at Position 1,  and the empty set {} for padding tells Partition to return the final sublist even if it's shorter than what is specified by the 2nd argument, and don't pad that final short list. Evaluate this code in your Notebook:

Qpredicates=Names@"*Q";

Print@Length@Qpredicates;

Partition[Qpredicates,4,4,1,{}]//TableForm

Here's the same result in a List.

Names@"*Q"

{AcyclicGraphQ,AlgebraicIntegerQ,AlgebraicUnitQ,AntihermitianMatrixQ,AntisymmetricMatrixQ,ArgumentCountQ,ArrayQ,AssociationQ,AtomQ,BinaryImageQ,BipartiteGraphQ,BooleanQ,BoundaryMeshRegionQ,BoundedRegionQ,BusinessDayQ,ColorQ,CompatibleUnitQ,CompleteGraphQ,CompositeQ,ConnectedGraphQ,ConstantRegionQ,ContinuousTimeModelQ,ControllableModelQ,CoprimeQ,DateObjectQ,DaylightQ,DayMatchQ,DeviceOpenQ,DiagonalizableMatrixQ,DigitQ,DirectedGraphQ,DirectoryQ,DiscreteTimeModelQ,DisjointQ,DispatchQ,DistributionParameterQ,DuplicateFreeQ,EdgeCoverQ,EdgeQ,EllipticNomeQ,EmptyGraphQ,EulerianGraphQ,EvenQ,ExactNumberQ,FileExistsQ,firedQ,FreeQ,GeoWithinQ,GraphQ,GroupElementQ,HamiltonianGraphQ,HermitianMatrixQ,HypergeometricPFQ,ImageQ,IndefiniteMatrixQ,IndependentEdgeSetQ,IndependentVertexSetQ,InexactNumberQ,integerGreaterThan4Q,IntegerQ,IntersectingQ,IntervalMemberQ,InverseEllipticNomeQ,IrreduciblePolynomialQ,IsomorphicGraphQ,KEdgeConnectedGraphQ,KeyExistsQ,KeyFreeQ,KeyMemberQ,KnownUnitQ,KVertexConnectedGraphQ,LeapYearQ,LegendreQ,LetterQ,LinkConnectedQ,LinkReadyQ,ListQ,LoopFreeGraphQ,LowerCaseQ,MachineNumberQ,ManagedLibraryExpressionQ,MandelbrotSetMemberQ,MarcumQ,MatchLocalNameQ,MatchQ,MatrixQ,MemberQ,MeshRegionQ,MixedGraphQ,MultigraphQ,NameQ,NegativeDefiniteMatrixQ,NegativeSemidefiniteMatrixQ,NormalMatrixQ,NumberQ,NumericQ,ObservableModelQ,OddQ,OptionQ,OrderedQ,OrthogonalMatrixQ,OutputControllableModelQ,PartitionsQ,PathGraphQ,PermutationCyclesQ,PermutationListQ,PlanarGraphQ,PolynomialQ,PositiveDefiniteMatrixQ,PositiveSemidefiniteMatrixQ,PossibleZeroQ,PrimePowerQ,PrimeQ,ProcessParameterQ,QHypergeometricPFQ,QuadraticIrrationalQ,QuantityQ,RegionQ,RegularlySampledQ,RootOfUnityQ,SameQ,SatisfiableQ,ScheduledTaskActiveQ,SimpleGraphQ,SquareFreeQ,SquareMatrixQ,StringFreeQ,StringMatchQ,StringQ,SubsetQ,SymmetricMatrixQ,SyntaxQ,TautologyQ,TensorQ,TreeGraphQ,TrueQ,UnateQ,UndirectedGraphQ,UnitaryMatrixQ,UnsameQ,UpperCaseQ,URLExistsQ,ValueQ,VectorQ,VertexCoverQ,VertexQ,WeaklyConnectedGraphQ,WeightedGraphQ}