How and why?
Mathematica Techniques
By Allan Hayes

For this issue I look at timing computations, using patterns, proving equality, getting the axes to go through the origin (in 2D) and precise evaluation.

Timing

I shall be timing computations a lot in this issue; so to reduce the clutter I define a utility.

Tim  =  Function[x, First[Timing[x]],  HoldAll];

which gives me

Tim[Integrate[Sum[Cos[n x], {n, 0, 10}], x]]

The attribute HoldAll is essential: without it we get

Tim2  =  Function[x, First[Timing[x]]];

Tim2[Integrate[Sum[Cos[n x], {n, 0, 10}], x]]

We cannot set attributes for a pure function in the "usual" way:

\(SetAttributes[Tim3, HoldAll]; 

    Tim3\  = \ Function[x, First[Timing[x]]]; \)

Tim3[Integrate[Sum[Cos[n x], {n, 0, 10}], x]]

And the attribute for Tim is not displayed by Attributes:

Attributes[Tim]

But we can use

\(SetAttributes[Tim4, HoldAll]; 
Tim4[x_]\  := \ First[Timing[x]]; \)

Tim4[Integrate[Sum[Cos[n x], {n, 0, 10}], x]]

Patterns

In the previous issue I dealt with how positions are used to target activity. A powerful feature of Mathematica is that we can also target by using the shape and properties of parts. In this section I will look at the use of names and patterns such as in the following examples

{3,  3.0,  4.0} /. x_Real /; x < 4 -> c[x]

{a[1, 1, 1], a[1, 1,  2]} /. a[1 ..] -> 3

f[x_, y__, z : p[5], x_] /; x < 3   := Times[y, z]f[2, 3, 4, p[5],  2]

Besides extending the language, names and patterns can improve clarity and efficiency: compare the following.

g[x_] := x[[1]] - x[[2]] + x[[3]]; 
gp[{a_,b_,c_}] := a - b + c;  

Do[g[{2,3,4}],{10000}]//Tim
Do[gp[{2,3,4}],{10000}]//Tim

hw[x_] := 
	Which[ 
		Head[x]=== Integer, "I", 
		Head[x]=== Real, "R",
		True,"?"
	]
	
hs[x_] := 
	Switch[ x,
		_Integer, "I", 
		_Real, "R",
		_, "?"
	]
	
hp[_Integer] := "I"   
hp[_Real] := "R"
hp[_] := "?"  

Do[hw[2],{10000}]//Tim
Do[hs[2],{10000}]//Tim
Do[hp[2],{10000}]//Tim

Names

Any expression can be given a name, which must be a symbol:

({{form, FullForm}, {name STARTCOMMENT a symbolENDCOMMENT : expr, Pattern[name,  expr]}})

Special forms like n_  for n:_ , and n__  d n:__ are allowed.

Names have a double role, like in mathematics :
	• for requiring equality
	• for passing values

For example, in evaluating g[a,b,c,d] after having entered

g[x_, y__, z : p[5], x_] /; x < 3   := Times[y, z]

If  g[a,b,c,d,e] is to match g[x_, y__, z:p[5], x_] then not only must d be p[5] but also,  a  must be the same as e, since they must match parts with the name x. 
Assuming that these conditions are met, then the name x in x_ is used to pass a to x in x < 3. 
Then, assuming that a < 3  evaluates to True, the names y in y__  and z in z:d are used to pass the sequence b,c to y  and the element d to z in Times[y,z], making Times[3,4, p[5]].

Blanks   _, __, ___

({{form, matches, FullForm}, {_, one expr, Blank[]}, {__, one or more exprs, BlankSequence[]}, {___, zero or more exprs, BlankNullSequence[]}, {_h, one expr with head h, Blank[h]}, {__h, one or more exprs with head h, BlankSequence[h]}, {___h, zero or more exprs with head h, BlankNullSequence[h]}})

As noted above, special forms are allowed when naming these: for example  n_  for n:_ and n__h for n:__h.

Repeats	p..  and  p...

({{form, matches, FullForm}, {p.., one or more matches for p, Repeated[p]}, {p..., one or more mathches for p, RepeatedNull[p]}})

Necessary spaces and parenthesis:

{MatchQ[{1, 1}, {1...}],       
 MatchQ[{1, 1}, {1 ...}],
 MatchQ[{1, 1}, {(1)...}]}

{MatchQ[{1, 2}, {_...}],  (* _. ⟶ Optional[Blank[]])*)
MatchQ[{1, 2}, {_ ...}],
MatchQ[{1, 2}, {(_)...}]}

Effect of naming

{MatchQ[{1, 2}, {(a : _ )...}],       
MatchQ[{1, 2}, {a : _ ...}]}

Alternatives	a|b|c

({{form, matches, FullForm}, {e1 | e2 | , ..., a match for at least one e1, Alternatives[e1, e2, ...]	}})

{a, b, c} /. a | c -> p

Some timings

Do[{a, b, c} /. a | c -> p, {1000}] // Tim

Do[{a, b, c} /. {a -> p, c -> p}, {1000}] // Tim

Do[{a, b, c} /. {x_ /; x == a || x == c  -> p}, {1000}] // Tim

Tests	p?(test)

({{form, matches, FullForm}, {p?(test), e matching p with test[e]⟶True
or
e1, e2, ... matching p with 
every test[ei] ⟶ True, PatternTest[p,  test]	
}})

If test is a symbol then the parentheses in p?(test) may be omitted:, but not otherwise:

MatchQ[{1, 1}, {__?Positive}]

MatchQ[{1, 1}, {_ ..?Positive}]

MatchQ[{1, 1}, {__?# > 0&}]

Several tests can be used in sequence.

Clear[g, x]

g[(x__?(# > 1&))?(# < 4&)] := A[x]

{g[1, 2], g[2, 3], g[3, 4]}

Condition  	p/;cond

({{form, matches, FullForm}, {p /; cond, e matching p with  \(cond\^*\) ⟶ True,
where \(cond\^*\) is the result of passing
values to cond from e., Condition[p,  cond]}})

MatchQ[{1, 1}, {x__ /; Length[{x}] > 1}]

Several conditions can be used in sequence

Clear[g]

g[x_, y_] /; x > 1 /; x + y < 6 := A[ x, y]

{g[1, 2], g[2, 3], g[3, 4]}

Notes

• Tests and conditions can be mixed

Clear[h]

h[(x__ /; Length[{x}] > 2 )?(# < 4&)  /; Plus[x] > 4] =  
      A[x];

{h[3, 1], h[3, 2, 1],  h[4, 3, 3]}

• Tests and conditions are evaluated even at held positions.

Clear[f, g]

f[HoldComplete[x_ /; x > 0]] :=  6f[HoldComplete[3]]

g[HoldComplete[x_?Positive]] :=  6g[HoldComplete[3]]

This can produce unexpected side effects:

Clear[x, y]\({Hold[Module[{x = \(-1\)}, x]] /. z_?Positive\  -> 5, \ 
	
      Hold[Module[{y = \(-2\)}, y]] /. z_ /; Positive[z]\  -> 5}\)

In spite of no change being apparent the assignments to x and y have been evaluated.

{x, y}

• In p?(test),  test[p] is evaluated but values are not passed from p to test.

Clear[h]

h[x_?(x > 0&)] :=  6h[3]

• Use a simple pattern shape without tests or conditions if possible, and prefer a test to a condition.

Do[{2.0, 2} /. _Integer -> A, {1000}] // Tim

Do[{2.0, 2} /. _?(Head[#] == Integer&) -> A, {1000}] // Tim

Do[{2.0, 2} /. x_ /; Head[x] == Integer -> A, {1000}] // Tim

Proving expressions equal

Proving two expressions equal can be difficult and time-consuming. Here are some generally applicable ideas:
Define

\(y\  = \ 1/\((x^7 + 1)\); \ 
z = D[Integrate[y, x], \ x]\ ; \)

Mathematics tells us that y == z, but how do we prove this by symbolic manipulation?
It's usually better to try to show that y-z is 0 and to use  some mathematics to reduces the calculation needed:

(w1 = Together[y - z,  Trig -> True]);  // Tim

(w2  =  Numerator[w1]);  // Tim

(w3  = First/@TrigFactorList [w2]);  // Tim

FullSimplify/@w3 // Timing

So, we have a zero factor in the numerator, whence y - z == 0.
In fact the zero factor comes from

w3[[2]]

Try RootReduce on this:

RootReduce[%]

This suggests that we might introduce another general step:

FullSimplify/@RootReduce/@w3 // Tim

FullSimplify is often slow, but it is needed here: Simplify is not enough.

Simplify[w3[[2]]]

In this example FactorList would have been quicker than TrigFactorList

First/@FactorList[w2] // Tim

If we are sure that trigonometric identies are not needed then we will not need the Trig variants like TrigFactorList or  the option Trig -> True.

Integrating a sum

You may run into the following kind of situation, for example in dealing with Fourier series.

sm  =  Sum[Sin[n x], {n, 0, 25}];

Integrate[sm, x] // Tim

We can speed this up by distributing the integration over the summands

Integrate[#, x]&/@sm // Timing

This may suggest that we  should always expand the integrand int a sum and then sum the integrals of the summands, but this is not a universal way of speeding up computation:

Integrate[(1 + Sin[x])^10, x] // Tim

Integrate[#, x]&/@Expand[(1 + Sin[x])^25] // Tim

AxesOrigin

gr = 
Plot[2 + Sin[x], {x,2,2Pi},
   	AxesOrigin -> {0,0}
];

How can we get the axes to go through the origin? The axes are determined by the PlotRange and this has not take into acount the origin that we set:

PlotRange[gr]

Here is one way of dealing with this

\(AxesAndOrigin[gr_, \ origin_] := 
	
      Show[gr, 
	AxesOrigin\  -> origin, \ 
	
        PlotRange -> \ 
          \((\(Through[{Min, Max}[#]]&\)/@Transpose[{PlotRange[gr], origin}])
            \)
]\)

AxesAndOrigin[gr, {0, 0}];

PreciseComputation

In the the last-but-one issue of Mathematica in Education and Research (Volume 7, No2) , Stan Wagon in his column "Programmig Tips" dealt with the problem of plotting a function that required high precision computation to get a correct value. The difficulty was that Plot inputs machine precision numbers. Here is an alternative appoach that modifies an expression so that its internal computations can be done to high precision, even when given low-precision input (the "WorkingPercision" is set)

PreciseExpression[expr_, x_, p_:$MachinePrecision] := 
Function[x, #1][##2]&[
	expr /. 
    r:(_Real|Complex[_,_Real]|Complex[_Real,_]):> 
    	SetPrecision[r, p], 
   	Sequence@@ 
    (Switch[#1, 
    	_?InexactNumberQ, SetPrecision[#1, p], 
        _?NumericQ, #1
   	 ]&)/@ Flatten[{x}]
]

This sets the precision of evey inexact number in expr to p and if x is a symbol then it replaces every free (term to be explained in a later column) occurence of x in expr with

Switch[x,
		_?InexactNumberQ, SetPrecision[x, p]
		_?NumericQ, x
	]

If x is a list of symbols {x1,x2, ...} then each free occurence of an xi in expr is similarly replaced.
The use of _?Numeric rather than  _   stops evaluation until x has a numerical value.
I use r:(_Real|Complex[_,_Real]|Complex[_Real,_]) rather than r_?InexactQ to avoid unwanted evaluations inside the test ( see the notes to section on Patterns)
A variant that gives a function, and can be faster, is given in the notes.

Examples 1

f1 = Sin[x] (10^x) Log[1 + (10^(-x))];

gives

Plot[f1, {x, 15, 20}] // Tim

Using the default precision ($MachinePrecision)

Plot[Evaluate[PreciseExpression[f1, x]], {x, 15, 20}] // Tim

Interestingly, we can evan lower the precision and get a good graph more quickly (this is presumably because setting the precision cause the internal software controls to be used in the computation)

Plot[Evaluate[PreciseExpression[f1, x, 3]], {x, 15, 20}] // Tim

But setting the precision to 2 causes a breakdown:

Plot[Evaluate[PreciseExpression[f1, x, 2]], {x, 15, 20}] // Tim

Example 2

This is taken from Stan Wagon's column:

f2  =  Normal[Series[Cos[x], {x, 0, 200}]];

Plot[f2, {x, 0, 76}] // Tim

Using default precision with PreciseExpression is not good enough

Plot[Evaluate[PreciseExpression[f2, x]], {x, 0, 76}] // Tim

To get a satisfactory result need the precision to be set to about 30

Plot[Evaluate[PreciseExpression[f2, x,  30]], {x, 0, 76}] // Tim

Example 3

PreciseExpression works with any number of variables.

Clear[y];

\(Plot3D[Evaluate[PreciseExpression[f1\ Sin[y], {x, y}]], 

        	{x, 15, 20}, {y, 0, Pi}] // Tim\)

Function form

PreciseExpression greatly enlarges the expression; PreciseFunction, below,  gives more concise and possibly faster result less flexible result.

PreciseFunction[expr_, x_, p_:$MachinePrecision] := 
  Function[x, Function[x, #1][##2]]&[
  	expr /. 
        r:(_Real|Complex[_,_Real]|Complex[_Real,_]):> 
    		SetPrecision[r, p] , 
   	Sequence @@ 
    (Switch[#1,
    	 _?InexactNumberQ, SetPrecision[#1, p], 
        _?NumericQ, #1] & 
     )/@ Flatten[{x}]
   ]

Compare

PreciseExpression[f[x] + x, x, p]

PreciseFunction[f[x] + x, x, p]

For an expression like f2 with many occurrences of x, PreciseFunction should be faster, but the difference is is often very slight:

Plot[Evaluate[PreciseFunction[f2, x,  30][x]], {x, 0, 76}] // Tim