How and Why?

Mathematica techniques

by Allan Hayes

In this column I shall try to demonstrate Mathematica techniques, from one-liners upwards and to explain how Mathematica works. I shall aim to give fairly full explanations, but, inevitably, some things will be have to be left, so I would encourage you to make the most of the excellent Help Browser.  If you have topics or questions that you would like to see dealt with please email them to me at hay@haystack.demon.co.uk.

Whether to use =  or :=

What is the difference between Set, =, and SetDelayed, :=, and when should one be use rather than the other?

The most important difference between the evaluation of l=r and of l:= r is that, with l = r , the right side, r, is evaluated before the rule is stored but with  l :=r it is not (its  evaluation is delayed until the stored rule is used - hence the name: DetDelayed).
This can alter thebehaviour quite significantly. 
Here are two examples, one where := is needed, the other where = is preferable

:= needed

f[x_]  =  Expand[x];

f[(a + b)^2]

We probably wanted (a^2) + 2  a b + (b^2). We didn't get it because  Expand[x] was evaluated to x before the rule was stored.
Here's what was stored.

?f

However with

f[x_]  :=  Expand[x];

the stored rule is

?f

and we get

f[(a + b)^2]

= preferable

With

g[a_]  :=  Integrate[x/Sqrt[1 + x^2],  {x, 0, a}]

the integration will not be done before the rule is stored and so time will be wasted in doing it every time that the rule is used.

?g

Do[g[n],  {n,  1, 5}] // Timing

However, with

g[a_]  =  Integrate[x/Sqrt[1 + x^2],  {x, 0, a}];

the integration is done before the rule is stored

?g

and the evaluation is much faster

Do[g[n],  {n,  1, 5}] // Timing

Another matter to bear in mind is the possible effect of existing definitions (in fact the problem in the first example was caused by the effect of the existing definiton of Expand).
With

h[x_]  :=  a x

we store

?h

So h[3] will give  (SuperDagger[a])3, where  (SuperDagger[a]) is the value of a when this rule is used.
Whereas, with

a  =  b;

h[x_]  =  a x;

the right side  a x will evalute to b x before storing;  so we store

?h

and,  h[3] will give  (SuperDagger[b])3 where  (SuperDagger[b]) is the value of b when the rule is used (by this time the rule a=b may have been removed)

Evaluation steps

In the box below,a  is an atom, (x^†)is the value of x.
Some features are clearly necessary:to be able to redefine a symbol a it is clear that a must not be evaluated; 
          and to be able to redefine a function, (h^†)[(e^†)..]  must not be further \(evaluated.
Apart\ fom\ \ the\ steps\ marked\ NB, \ 
        where\ the\ right\ sides\ evalute\ before\ the\ left, \ 
        the\ evaluation\) is entirely regular({{a := c , store,        , h[e.. ],  }, {Null, output,  , (h^†)[e..]  :=  c,  }, { ,  ,      , (h^†)[(e^†)..]  :=  c, store}, { ,  ,  , Null, output}, { ,  ,  ,  , \(\ 
	\)}, {a = c,  ,  , h[e.. ] = c,  }, {a=(c^†)(NB), store,  , h[e.. ]=(c^†) (NB),  }, {(c^†)†, output,  , (h^†)[e.. ] = (c^†),  }, { ,  ,  , (h^†)[(e^†).. ] = (c^†), store}, { ,  ,  , (c^†)†, ouput}})(1) Storing takes place only if s or the left –most atom of his an unprotected symbol  	(the stored assignment is tagged by this symbol). 
(2) Evaluation of (h^†)[e..] is subject to any attributes of (h^†).
(3) The ouput is generated even if nothing is stored ( (c^†)† may not be the same  (c^†)).

Usually, (3) will not be significant, but it does show up in the following.

a =  {b, b = 2}

?a

> (Rule) and :> (RuleDelayed)

The differences are very much the same as between = and := : 
	expr/.l -> r evaluates r to (SuperDagger[r]), say, and then replaces every l in expr with (SuperDagger[r])
	expr/.l :> r replaces every l in expr with r (unevaluated

Clear[a, b, ]

{(a + b)^2} /. x_Power ->  Expand[x]

{(a + b)^2} /. x_Power :>  Expand[x]

Here is a use of :> to solve a recent problem sent to mathgroup.
Given

\(gt = {{3, 5}, {{}, 6}, {7, {}}, {{}, {}}, {{}, 56}, {{}, 3}}; 

    b = {a1, a2, a3, a4, a5, a6}; \)

We can replace successive occurrences of {} in gt with successive entries in b as follows

\(n = 1; 
gt //. {} :> \ b[\([\(n++\)]\)]\)

Evaluation steps

Evaluation is absolutely regular - left to right then replacement, then final evaluation.

({{e /. a :> c ,   , e /. a -> c }, {(e^†) /. a :> c,  , (e^†) /. a -> c}, {(e^†) /. (a^†) :> c,    , (e^†) /. (a^†) -> c}, {(e^†) with  (a^†)replaced by c,  , (e^†) /. (a^†) -> (c^†)}, {(SuperDagger[((e^†) with  (a^†)replaced by c)]),  , (e^†) with  (a^†)replaced by (c^†) }, { ,  , (SuperDagger[((e^†) with  (a^†)replaced by (c^†))])}})

New in 3.0x

Split

?Split

So, the splits are made between consecutive elements which do not satisfy the test.

Split[h[1, 1, 1, 4, 4, 4, 3, 3]]

\(Split[h[1, 1, 1, 4, 4, 4, 3, 3], \ 
		Abs[#1 - #2] < 2&]\)

Application: removing repeats

\(removeRepeats[x_] := 
First/@Split[Sort[x]]\)

removeRepeats[ {3, 1, 2, 2, 1, 3, 2, 1}]

However, we have lost the ordering

Here is a function that gives us the first occurences of members in the original order

FirstOccurences[x_] := 
      Last/@Sort[
          Reverse/@
            First/@
              Split[Sort[MapIndexed[List, x]], #[[1]] == #2[[1]]&]]

FirstOccurences[{3, 1, 2, 2, 1, 3, 2, 1}]

The steps used are:

MapIndexed[List, {3, 1, 2, 2, 1, 3, 2, 1}]

Sort[%]

Split[%,  #[[1]] == #2[[1]]&]

First/@%

Reverse/@%

Sort[%]

Last/@%

ReplaceList

?ReplaceList

ReplaceList[{a, a, a, a}, {x___, a, y___, a, z___} -> {x, b, y, c, z}]

Part

The way in which Part interacts with Set has been enhanced

(M  =  Table[0, {4}, {5}]) // MatrixForm

M[[{2, 3}, {2, 4, 5}]]  = {{1, 2, 3}, {4, 5, 6}}; 
    M // MatrixForm

M[[{2, 3}]]  =   M[[{3, 2}]]; M // MatrixForm

\(\(M[\([2]\)]\  = {1, \(-4\)}.M[\([{2, 3}]\)]; M // MatrixForm
\)\)

We are not restricted to matrices

a  =  h[f[1, 2], g[3, 4]]a[[1, {1, 2}]] =  {p, q}; a

Form puzzles

I'll discuss MatrixForm but the general principles apply to othe "wrappers" like TableForm, FortranForm ... .

$Line  =  0

MatrixForm[3  + a]

MatrixForm only determines the form of the display, not the output.
The output is the same as for 3 +a: it does not involve MatrixForm:

%

The output cell label  Out[.]//MatrixForm=  is meant to indicate this situation.

This stripping of MatrixForm occurs only when it is the Head of the expression to be displayed, if it is deeper in the expression then it is retained and can cause surprises.

3  + MatrixForm[a]

The reason for this is that MatrixForm is retained around 3+a, and prevents the addition of 3 going inside the lists, besides affecting the formating of the  display. Notice that the output cell label is simply Out[. ]=.

Here is another consequence

f  =  MatrixForm[a]

The rule is stored before the display of the output and keeps MatrixForm on the right side

?f

So, the following behaves like 3 + MatrixForm[a].

3  +  f

To assign the plain list structure a to f but display in MatrixForm we can use

MatrixForm[f =  a]

or

(f =  a) // MatrixForm

Thus,

?f

No more aliasing?

If we look up Plot in the HelpBrowser we find:

◼ Plot initially evaluates f at a number of equally spaced sample points specified by PlotPoints. Then it uses an adaptive algorithm to choose additional sample points, attempting to produce a curve in which the bend between 
            successive segments is less 	than MaxBend. It subdivides a given interval by a factor of at most PlotDivision. 
		
◼ You should realize that with the finite number of sample points used, it is possible for Plot to miss features in your function. To check your results, you should increase the setting for PlotPoints.

Now,  for Plot,   
    the default setting of PlotPoints is PlotPoints -> 25;  
    so if we plotted

Plot[Sin[x Pi], {x, 0, 24 }]

then the table of values at the initial points would be

Table[Sin[x Pi], {x, 0, 24}]

Consequently, Mathematica should think that it had been given the constant zero function and would not take any more action. This sort of thing is called aliasing, and it is what happened up to Version 3.0
But now we get

Plot[Sin[x Pi], {x, 0, 24 }]

The trick that has been introduced is a simple one - we don't divide evenly.
Let's find the initial points actually used: prevent the adaptive subdivision, by plotting 0, and accumulate the points.

\(lst\  = \ {}; 
Plot[AppendTo[lst, x]; 0, {x, 0, \ 24}]; \)

lst

Not

Table[x , {x, 0., 24.0}]

The points in lst are far from  evenly spaced.
Here are the difference between successive values.

Drop[lst, 1]  -  Drop[lst, -1]

But they are repeatedable - the following computation (with display prevented by DisplayFunction->Identity)was made several days after the preceding.

\(lst\  = \ {}; 

    Plot[AppendTo[lst, x]; 1, {x, 0, \ 24}, \ 
DisplayFunction -> Identity]; 
    \)

lst

ABOUT THE COLUMNIST