How and why?
Mathematica Techniques
By Allan Hayes

For this double issue, I've included a collection of techniques related to the following Mathematica concepts: NestWhile and NestWhileList, from Mathematica  4.0;  masking overspills in graphics; unbalanced and expanding brackets, aligning text; some speedups; using HoldPattern; an Until function; using $Post to get MatrixForm.automatically.Arrows on axes; one-liners; N[x,p] and automatic precision;  ContourPlot, contour heights and region coloring; FindRoot and HoldAll;  Orderless and Flat.

Arrows on Axes

Programming note: the package,  Graphics`Arrow`, is loaded internally by the function defined below, so we have to use  the long name Graphics`Arrow`Arrow .rather than just Arrow. This is because contexts are assigned before evaluation, and Arrow would be given the current context, usually "Global`", and would not access the package definitions (we would get a message like when we enter Arrow before loading the package). 
Options are passed on to Arrow.

ArrowsOnAxes[gr_, opts___?OptionQ]:=
	Module[{ox,oy,xmin, xmax,ymin,ymax, xst,yst},
	Needs["Graphics`Arrow`"]; 
	{{ox,oy},{{xmin, xmax},{ymin,ymax}},{xst,yst}}=
 Last/@AbsoluteOptions[Graphics[gr], {AxesOrigin, PlotRange, AxesStyle}];
	Show[gr,
	Prolog->{
	Append[xst,
Graphics`Arrow`Arrow[{xmin,oy}, {xmax,oy},opts]],
	Append[yst,
Graphics`Arrow`Arrow[{ox,ymin}, {ox,ymax}, opts]]}
	]
];

gr1=Plot[Sin[x],{x,0, 7},
	AxesStyle-> Hue[0]];

ArrowsOnAxes[gr1, HeadCenter->.7];

One-Liners

Find values of  a  and  b  for which (e^(-2 x)) (10 b-x+6 a (5 x+1)-3) is always zero.

expr = E^(-2 x)((-3) + 10 b - x + 6 a(1 + 5 x)) ;

SolveAlways[expr== 0,{x}]

expr/.%//Simplify

In the following delete the events that fall into one of the gaps.

gaps={{1,5.5},{2,4.3},{2.7,3.1},{3.002,4.007},{10.001,10.007},{10.101,11.001},{11.007,12.0}};
events={6.7,8.9,2.3,2.789,10,11.002,10.115,3.02,2.75};

Two solutions:

DeleteCases[events,n_/;
IntervalMemberQ[Interval@@gaps,n]]

DeleteCases[events,n_/;Or@@((#≤n≤#2)&@@@gaps)]

N[x,p]

An important feature for numerical work is that if expr is an exact numeric expression (built up from integers, rationals, constants like π and e, complex numbers with such real and imaginary parts, and functions having the attribute NumericFunction), then,  if p>$MachinePrecision, N[expr,p] will increase its working precision to try to return an output with precision p (this increase can only be up to $MachinePrecision + $MaxExtraPrecision but the default value of $MaxExtraPrecision is 50 but it can be set by the user). The values obtained are cached; the cache can be cleared with Developer`ClearCache["Numeric"].

Let's try it out on finding  Sin[10^3000].

For convenience define

n = 10^3000;

N[Sin[n]]

This looks suspicious, especially when we get

Precision[%]

Lets try

N[Sin[n],25]

Increase $MaxExtraPrecision , as suggested

$MaxExtraPrecision = 5000;

We get

N[Sin[n],25]//Timing

Precision[%[[2]]]

A repeat of the calculation is very quick because of the caching.

N[Sin[n],17]//Timing

The cached value is used to speed up calculations of Sin[n] to other precision values

N[Sin[n],50]//Timing

N[Sin[n], 17]//Timing

Clear the cache

Developer`ClearCache["Numeric"];

Let's build a function to display the increases in the working precision and the precision of the intermediate outputs.

sn[x_]:= (Print[{Precision[x],Precision[Sin[x ] ]}]; Sin[x ])

N[sn[n], 25]

Problem!: sn[n] evaluated to Sin[10^3000] (and printed) before N[.,25] was used. Try

Clear[sn]
sn[x_Real]:= (Print[{Precision[x],Precision[Sin[x ] ]}]; Sin[x ])

N[sn[n], 25]

Problem!: sn needs to be have the attribute NumericalFunction for the increase in working precision to be activated - the output is Sin[ N[n,25]].

Attributes[sn]= NumericFunction;

N[sn[n],25]

Thus, the working precision that was eventually used was 3585; it gave a result with precision 585, which was then reduced to precision 25:

SetPrecision[Sin[SetPrecision[n, 3585]],25]

The expression really has to be exact for this adjustment to work:

N[Sin[1.0n],25]//Timing

See Stan Wagon: Volume 7, No 2, Programming Tips, Bending Plot to Your Needs, for an application of these features.

ContourPlot

How are the contour heights chosen and how are the colors of the regions defined by the contours chosen?
In Version 4.0, a new option, ColorFunctionScaling  lets us choose whether or not to have the inputs to the color function scaled.

First, some examples (default settings of selected options are given in comments).
The color function is set so that we get a printout of the numbers that are supplied to it.

expr=2/(.1+(x+.5)^2+y^2);

With Contours -> 4.

cgr1 = 
ContourPlot[expr, {x, -3, 3}, {y, 0, 2}, 
   PlotPoints -> {4, 3},	(*{15,15}*)
   Contours -> 4, 			(*10*)
   ColorFunction -> ((Print[#];Hue[2/3#])&), 
   FrameLabel->{"x", "y"}           
];

The color function inputs for the regions between contours came first, in order of increasing height, then the inputs for lowest and the highest regions in that order.

With Contours -> 4  and ColorFunctionScaling->False.

cgr1unsc = Show[cgr1, ColorFunctionScaling->False];

With Contours -> {-4,.6, 1.8, 4}

cgr2=Show[cgr1, Contours ->{-4,.6, 1.8,4}] ;

Only two contours were drawn, actually at heights 0.6 and 1.8.  Height -4 is clearly impossible, but why not 4, when the maximum value of  expr in the plotted range is  20 ? - I'll come back to this later

With Contours -> {-4,.6, 1.8, 4} and ColorFunctionScaling->False.

cgr2unsc = Show[cgr2, ColorFunctionScaling->False];

How are contour heights and colors between contours calculated?

Contour heights
First, find the z-PlotRange, {zmin, zmax}.

{zmin,zmax}=PlotRange[cgr1][[3]]

The contour heights to try, say {u1,u2,..,un}, are: 
- with Contours -> k (default  k =10): 
	 k heights  dividing  {zmin, zmax} into k+1 equal intervals
- with Contours -> {z1,z2,..}: 
	the (Subscript[z, i]) that are between zmin and zmax , if any;  
	otherwise the least one above zmax if there is one; 
	otherwise the greatest below zmin.

Check against the examples.

For cgr1  (Contours -> 4) this gives

U1=zmin + Table[k (zmax-zmin)/(4+1), {k,1,4}]

For cgr2 ( Contours ->{-2,.3, .9,2}) it gives

U2 = {.6,1.8};

Color function inputs

Input to cf  for the region: 
	- between contour heights ui and ui+1 is   (ui + ui+1)/2
	- below the lowest contour is  (zmin +  u1)/2
	- above the highest contour is (un + zmax)/2
With  ColorFunctionScaling -> True (the default) the numbers above are linearly scaled so that the least one becomes 0 and the greatest becomes 1. 
So, with Contours -> k and ColorFunctionScaling ->True the color function inputs are determined solely by k.

Check against the examples (for computational convenience I give the unscaled versions first;)

cgr1unsc

Join[ Plus@@@Partition[U1,2,1],
{zmin+First[U1]},{zmax + Last[U1]}]/2

cgr1

Interpolation[{{%[[-2]],0},{%[[-1]],1}},InterpolationOrder ->1]/@%

cgr2unsc

Join[ Plus@@@Partition[U2,2,1],
{zmin+First[U2]},{zmax + Last[U2]}]/2

cgr2

Interpolation[{{%[[-2]],0},{%[[-1]],1}},InterpolationOrder ->1]/@%

The z-PlotRange

Here are some effects of setting the z-PlotRange

Show[cgr1, PlotRange-> {-10,10}];(*Contours -> 4*)

Notice that five color function inputs were generated in  spite of only two contour regions being visible.

Show[cgr2, PlotRange-> {-4,4}];(*Contours ->{-4,.6,1.8,4}*)

The contour at height 4 was added, there can be no contour at height -4.

Show[cgr2, PlotRange-> {0,1}];(*Contours ->{-4,.6,1.8,4}*)

Only the contour at height 0.6 was drawn, and only the color function inputs for the two contour regions defined by the contour actually drawn are generated.

Specifying the colors of contour regions

Suppose that we wish to draw several contour plots with  Contours -> {z1,z2,..}and have the same colors correspond to the the same heights. The following function will construct an appropriate color function

MakeColorFunction[cntrs_?OrderedQ, clrs_]/;Length[clrs]==Length[cntrs]+1:=
Function@@{Which@@Flatten[{{#≤ cntrs[[1]],clrs[[1]]},Table[{cntrs[[i]]<#≤cntrs[[i+1]],clrs[[i+1]]},{i,1, Length[cntrs]-1}],{cntrs[[-1]]≤#,clrs[[-1]]}}
]}

For example, we get

MakeColorFunction[{z1,z2,z3},{c1,c2,c3,c4}]

Let's make a plot.

<<Graphics`Colors`

cntrs= {-.9,-.6,-.3,0,.3,.6,.9};
clrs = {Red,Blue,LightSeaGreen,YellowBrown,Khaki,Green, Purple,Yellow};

ContourPlot[(Cos[x]+Sin[y])/2,{x,0,2Pi},{y,0,2Pi},
Contours->  cntrs,
ColorFunction -> (MakeColorFunction[cntrs,clrs]),
ColorFunctionScaling ->False];

Internals, z-PlotRange

We already noticed that {zmin,zmax}, {-0.0471186,1.93186}, is not the same as {exprmin,exprmax}, which is {0, 20}).  There are three reasons for this:
1. because of PlotPoints -> {4, 3} the values of expr are sampled on the 4 × 3 grid (see below);
2. Mathematica may (and here does) clips the range of these values;
3. The result of 2. is enlarged by 5% - this allows for frame ticks.

The 4 × 3 grid referred to in 1. above is

samp=Table[{x,y},
{y,0,2,(2-0)/(3-1)},
	{x,-3,3,(3-(-3))/(4-1)}]

Here are these points superimposed on the cgr1(I block the printing of color function inputs).

Block[{Print},Show[cgr1, Graphics[{PointSize[.05],Map[Point, samp,{2}]}]]];

The corresponding values of expr are

Table[expr,{y,0,2, (2-0)/(3-1)},{x,-3,3,(3-(-3))/(4-1)}]

This table is stored in the ContourGraphics object cgr1 as its first entry

cgr1[[1]]

We have

{Min[%],Max[%]}

But

PlotRange[cgr1][[3]]

Notice that the only metric information that the ContourGraphics object carries is this table and the value of the option MeshRange (which is used to recover the original x and y ranges)

Options[cgr1, MeshRange]

The construction of the contour lines is based solely on the table, not the original function.

FindRoot

FindRoot has the attribute HoldAll, this lets it check and preprocess the input (it also, as we shall see, allows the variables being solved for to be localized and shielded against any existing assignments to them). 
The evaluation of FindRoot[body,{t,n}] begins by matching  the unevaluated body to the first possible one of the forms {lhs==rhs},{expr}, lhs == rhs, expr. Whichever one matches it will then proceed on the assumption that  lhs-rhs, respectively expr, will give a numerical value when t is assigned a numerical value (the HelpBrowser only mentions the equation forms so perhaps we should stick with these).

In the following

Clear["`*"]

f[t_] := Cos[a t]

FindRoot[f[t]==0 /. a->0.4,{t,1} ]

the match of
(f[t]==0)/. a->0.4
is with {expr} (it is not an equation), so  FindRoot assumes that,  for example, f[t]==0)/.a->0.4  will give a numerical value when t is given the value 1.0 - and it does not do this -  it gives False.

We could get round this by forcing FindRoot to evaluate body before matching.

FindRoot[Evaluate[f[t]==0 /. a->0.4],{t,1} ]

But this way would give trouble if t already had an assignment - this shows up the localizing role of the holding attribute:

t= s;
FindRoot[Evaluate[f[t]==0 /. a->0.4],{t,1} ]

The variable t in the body was replaced by s but the t in {t, 20} was not replaced,

It is better to use

FindRoot[(f[t]/. a->0.4)==0,{t,1} ]

However s has now been substituted for t , after the solution was found

We run into similar problems with function like Plot and Sum that hold entries.

Solve, which does not hold its its entries, still has the substitution problem.

Solve[t== 4,t]

There seems to be no way round this other than clearing values or by working entirely inside Block[{t},..]

Orderless and Flat

The behavior of the attributes Orderless and Flat (which is often used with the attribute OneIdentity) frequently puzzles users. I'll try to sort things out.

Firstly, Flat and Orderless have two independent effects - to rearrange expressions and to modify pattern matching.

Let's look at Flat:

ClearAll[f]

Attributes[f] = {Flat};

Rearranging (flattening) expressions is simple

f[f[1],2]

Modifying pattern matching is more complicated:

f[1,2,3]/.f[x_,y_]-> F[x,y]

f[1]/.f[x_]-> F[x]

It seems that in the first example, although f[1,2,3] does not match f[x_,y_], the pattern matcher uses f[f[1],f[2,3]], which will flatten to f[1,2,3]. But, why does it not use the simpler f[1, f[2,3]]? And why in the second example does it match f[x_] to f[f[x]] and give F[f[x]]? It seems that only patterns of the form f[f[___],f[___],...] are used.
The simpler behavior is obtained if the attribute OneIdentity is added.

Attributes[f] = {Flat, OneIdentity};

f[1,2,3]/.f[x_,y_]-> F[x,y]

f[1]/.f[x_]-> F[x]

With assignments, another feature of Flat emerges:

ClearAll[f];
Attributes[f]={Flat} ; (*attribute set before assignment*)
f[1, x_] = F[x];
{ f[f[1]], f[1,2,3] }

but

ClearAll[f];
f[1, x_] = F[x];
Attributes[f]={Flat} ;(* attribute set after assignment*)
{ f[f[1]], f[1,2,3]}

In both examples f[f[1]] is flattened to f[1]—  because when the rule is applied the evaluator finds Flat in Attributes[f]. But only in the first example is the stored rule used on f[1,2,3]. The reason for this is that when the rules were stored, f was evaluated and the then value of Attributes[f], {Flat}in the first example and { } in  the second one, was invisibly attached to f;  and it is this attachment, not the current value of Attributes[f], that the pattern matcher consults when it has to decide whether or not the rule is applicable.

We can take this further

ClearAll[f];
Attributes[f]={Flat} ; (*attribute set before assignment*)
f[1, x_] = F[x];
Attributes[f] = {};
{f[f[1]], f[1,2,3] }

f[f[1]] was not flattened, because when it was evaluated, Attributes[f] gave {}. But the stored rule was used on f[1,2,3],  because when the rule was stored Attributes[f] gave {Flat}.

The difference between  f's with and different attachments is recognized when new rules are made

ClearAll[f]
f[1, x_] = G[x];
Attributes[f] = {Flat};
f[1, x_] = F[x];

The second rule, with {Flat} attached to f, has not replaced the first one, with { } attached.

?f

We can replace the second rule

f[1,x_]= F2[x]; 
?f

and the first one

Attributes[f] = {};
f[1,x_]= G2[x];
?f

But Unset ( =.) does not maintain the recognition of the difference;  for although the following deletes the first rule, with { } attached to f:

f[1,x_]=.;
?f

nevertheless, a repeat of f[1,x_]=.removes the remaining rule, which has {Flat} attached to f:

f[1,x_]=.;
?f

Orderless

Orderless behaves similarly, but after

ClearAll[f];
Attributes[f] = {Orderless};
f[x__/;!OrderedQ[{x}]] = F[x];

there is more to the following than meets the eye in

f[2,1]

I'll edit and annotate the information provided by TracePrint.

f[2,1]//TracePrint

And, perhaps surprisingly,

Attributes[f] = {};
f[1,2]

Which comes about in the following way

f[1,2]//TracePrint

Replacement rules

The preceding attachment complications do not arise when replacement rules are used directly as at the beginning  of this section, but they do crop up when there the rules are stored.

ClearAll[f];
Attributes[f] = {Flat};
rl = (f[1,x_]->F[x]);
f[1,2,3]/.rl

Attributes[f] = {};
f[1,2,3]/.rl

However, if the right side of the assignment is not evaluated then nothing is attached to f when the rule is stored, and Attributes[f] is consulted when the replacement is subsequently performed with the rule. 
The following checks that {Flat} is not attached to f in the stored rule.

ClearAll[f];
Attributes[f] = {Flat};
rl := (f[1,x_]->F[x]);
Attributes[f] = {};
f[1,2,3]/.rl

And this checks that Attributes[f] is consulted when the rule is used.

Attributes[f] = {Flat};
f[1,2,3]/.rl