The Way of Mathematica

Sunday, April 24, 2016

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.

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.

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.

Saturday, April 23, 2016

Computing the Mandelbrot set is a classic example of the need to Compile a Mathematica function (see Ruskeepää, p. 451). You just need to get the syntax right for a successful Compile. Below, the first argument is the parameter, c, and its type, Complex. The entire Module is the second argument, the expression to be compiled. Definitely try this at home!

The CompiledFunction can be applied to individual numbers. The criteria for inclusion in the Mandelbrot set is, after a limited number of iterations of z = z^2 + c (in this case 50), is its absolute value less than a given number (in this case 2). The absolute value of a complex number is its Euclidean distance from the origin. If the absolute value of a number exceeds 2 before 50 iterations, e.g. 4, it is excluded.

mandelbrot[0.2+I]4mandelbrot[0.2+0.2I]50

Let's do some heavy-duty computing. As a precaution we'll increase the RAM allocated to Mathematica.

Needs@"JLink`";ReinstallJava[JVMArguments->"-Xmx6000m"]LinkObject[[]

Start with a couple of individual plots. I use ImageSize->Full so when you do this yourself a large image will automatically be generated. On my Dell Optiplex 900 these take 2.68 seconds each.

I'm kind of curious about all the built-in ColorFunctions. Let's plot Mandelbrot sets using each of them. Here they are. I reverse the order since I like some of the ones at the end best.colorFunctions=ColorData@"Gradients"{AlpineColors,Aquamarine,ArmyColors,AtlanticColors,AuroraColors,AvocadoColors,BeachColors,BlueGreenYellow,BrassTones,BrightBands,BrownCyanTones,CandyColors,CherryTones,CMYKColors,CoffeeTones,DarkBands,DarkRainbow,DarkTerrain,DeepSeaColors,FallColors,FruitPunchColors,FuchsiaTones,GrayTones,GrayYellowTones,GreenBrownTerrain,GreenPinkTones,IslandColors,LakeColors,LightTemperatureMap,LightTerrain,MintColors,NeonColors,Pastel,PearlColors,PigeonTones,PlumColors,Rainbow,RedBlueTones,RedGreenSplit,RoseColors,RustTones,SandyTerrain,SiennaTones,SolarColors,SouthwestColors,StarryNightColors,SunsetColors,TemperatureMap,ThermometerColors,ValentineTones,WatermelonColors}Length@colorFunctions51I'll make a Grid of them that will correspond to the Grid of the plots themselves so we can easily look up the ColorFunction of one we like by finding its name in the corresponding place in the Grid. Although each row will have 3 plots and 3 divides 51 evenly, by habit I use the syntax for Partition that would not drop an overhang of a ragged List (the empty List {} in the position of the padding argument).Reverse@colorFunctions//Partition[#,3,3,1,{}]&//Grid

The iterator, predefined, in this Table iterates over all the ColorFunctions. Note you do not need to use the cumbersome {i, Length@colorFunctions} requiring the iterator to be a number.Table[DensityPlot[mandelbrot[x+I y],{x,-2,1},{y,-1.5,1.5},ColorFunction->predefined,PlotPoints->200,Mesh->False,FrameTicks->{{-2,-1,0,1},{-1,0,1}},ImageSize->Medium],{predefined,Reverse@colorFunctions}]//Partition[#,3,3,1,{}]&//Grid

Thursday, April 21, 2016

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]:= ?IntegrateTrapezoidIntegrateTrapezoid[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@data1Out[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.

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

Thursday, February 4, 2016

Attributes modify the behavior of entire classes of functions. Some Attributes cause functions to have well-known mathematical properties such as commutativity, associativity, or idempotency, notably for valid pattern-matching. Others control the way the Wolfram Language evaluates or does not evaluate a function's arguments, which if you study that you will see is at the very core of Wolfram Language operation.

Listable is a valuable functional language Attribute since it incorporates Thread into a function so it automatically maps over a List and instead of taking a List apart, we apply a function to it as a gestalt.

NumericFunction is an Attribute of all Wolfram Language mathematical functions and is True when all of their arguments are numerical. I suspect it is integral to the interaction of many functions with the underlying knowledge base of Mathematica and therefore to many functions' valid behavior (e.g. Piecewise – always define piecewise mathematical functions with Piecewise).

Attribute: OneIdentity

OneIdentity is the mathematical property of idempotency, which dictates that applying a function any number of times to its argument (i.e. Nest) returns the argument unchanged. Examples are adding zero to any number, or multiplying any number by one, or raising any number to 1st or 0th power (except 0 - someone explain that inconsistency to me!).0^0During evaluation of In[549]:= Power::indet: Indeterminate expression 0^0 encountered. >>Out[549]= Indeterminate

The Doc Center says that OneIdentity is necessary for pattern-matching to work with the default of Times, Plus, and Power, which is a period. I.e., x_+ y_. => x + 0, x_ * y_. => x * 1, and x_^y_. => x^1.

MatchQ[x, Times[n_. , x_]]True

You can see that OneIdentity removes nested Heads. It says to the evaluator, 'If the given pattern x_ matches any of these patterns, return True.'

Attributes modify the behavior of entire classes of functions. Some Attributes cause functions to have well-known mathematical properties such as commutativity, associativity, or idempotency, notably for valid pattern-matching. Others control the way the Wolfram Language evaluates or does not evaluate a function's arguments, which if you study that you will see is at the very core of Wolfram Language operation.

Listable is a valuable functional language Attribute since it incorporates Thread into a function so it automatically maps over a List and instead of taking a List apart, we apply a function to it as a gestalt.

NumericFunction is an Attribute of all Wolfram Language mathematical functions and is True when all of their arguments are numerical. I suspect it is integral to the interaction of many functions with the underlying knowledge base of Mathematica and therefore to many functions' valid behavior (e.g. Piecewise – always define piecewise mathematical functions with Piecewise).

Attribute: Flat

Flat is the mathematical property of associativity, which dictates that different groupings of numbers to specify the order in which the operation is performed, return the same result. E.g. ((a+b)+c) = (a+(b+c)) and (2x3)x4 = 2x(3x4).

a+(b+c)==(a+b)+cTrue

Thomas Bahder uses a neat trick to show how Attributes work and how some affect pattern-matching. Here is the trick revealing the mechanism of Flat. Bahder uses a Condition that will never be True, which causes Wolfram Language to try all permutations of the function's arguments that are valid under Orderless. In the false predicate he inserts a Print statement to reveal how the Wolfram Language evaluator works. HoldForm is necessary to stop the evaluator from putting the arguments into standard order before printing them.

Flat tells the Wolfram Language evaluator, 'If the given pattern x_ + y_ matches any of the following patterns, return True.'

a+b+c/.x_+y_:>u/;Print@FullForm[Plus[x,y]//HoldForm]HoldForm[Plus[a,Plus[b,c]]]HoldForm[Plus[b,Plus[a,c]]]HoldForm[Plus[c,Plus[a,b]]]HoldForm[Plus[Plus[a,b],c]]HoldForm[Plus[Plus[a,c],b]]HoldForm[Plus[Plus[b,c],a]]HoldForm[Plus[a,b]]HoldForm[Plus[b,a]]HoldForm[Plus[a,c]]HoldForm[Plus[c,a]]HoldForm[Plus[b,c]]HoldForm[Plus[c,b]]a+b+cHere are related posts on Orderless and OneIdentity, and this post shows how to reveal all built-in functions with Attributes categorize by their Attribute.

This is an excellent out-of-print book but there is plenty to learn from it. I'd compare it to David Wagner's Power Programming in Mathematica in its insight to how the Wolfram Language works.

Attributes modify the behavior of entire classes of functions. Some Attributes cause functions to have well-known mathematical properties such as commutativity, associativity, or idempotency, notably for valid pattern-matching. Others control the way the Wolfram Language evaluates or does not evaluate a function's arguments, which if you study that you will see is at the very core of Wolfram Language operation.

Listable is a valuable functional language Attribute since it incorporates Thread into a function so it automatically maps over a List and instead of taking a List apart, we apply a function to it as a gestalt.

NumericFunction is an Attribute of all WL mathematical functions and is True when all of their arguments are numerical. I suspect it is integral to the interaction of many functions with the underlying knowledge base of Mathematica and therefore to many functions' valid behavior (e.g. Piecewise – always define piecewise mathematical functions with Piecewise).

Attribute: Orderless

Thomas Bahder uses a neat trick to show how Attributes work and how some affect pattern-matching. Here is the trick revealing the mechanism of Orderless, which is equivalent to the mathematical commutative property, that is, x + y = y + x.

Bahder uses a Condition that will never be True, which causes WL to try all permutations of the function's arguments that are valid under Orderless. In the false predicate he inserts a Print statement to reveal how the WL evaluator works. HoldForm is necessary to stop the evaluator from putting the arguments into standard order before printing them.

Orderless says to the Wolfram Language evaluator, 'If given pattern matches any of these patterns, return True.'

a+b+c/.x_+y_+z_:>u/;Print@FullForm[Plus[x,y,z]//HoldForm]HoldForm[Plus[a,b,c]]HoldForm[Plus[a,c,b]]HoldForm[Plus[b,a,c]]HoldForm[Plus[b,c,a]]HoldForm[Plus[c,a,b]]HoldForm[Plus[c,b,a]]a+b+cHere are related posts on Flat and OneIdentity. This post lists all Wolfram Language functions that have any Attribute except Protected.

This is an excellent out-of-print book but there is plenty to learn from it. I'd compare it to David Wagner's Power Programming in Mathematica in its insight to how the Wolfram Language works.

Now construct a Table of functions with any of the Attributes. Use Select to apply the predicate hasAttributeQ to each built-in function, and return it if it has an Attribute. I could not get the column label in TableHeadings to work so I used Labeled to add one to each group. The formatting and centering of the Label didn't translate to this blog format but if you execute the code you'll see the result.

Table[Select[hasAttributeQ[#,attribute]&]@allSystemSymbols//TableForm[#,TableHeadings->{Automatic,None}]&//Labeled[#,"\n"<>ToString@attribute,Top,LabelStyle->Directive[Large,Bold,FontFamily->"Times New Roman"]]&,{attribute,allAttributes}]//Column