Clever Little Programs

EvaluateAt

The function in the next cell takes an expression and a list of positions,
and evaluates in place the parts at the specified positions. This is taken
verbatim from "further examples" for ReplacePart in the Help Browser.

In this case, {1,2} specifies the second element of the held list, and {1,-1,1} the first part of the last element.

Evaluate Pattern

The function below evaluates all parts of a held expression that match a  certain pattern. This is based on some code Robby Villegas of Wolfram  Research discussed at the 1999 Developer Converence.  See "Working With  Unevaluated Expressions" posted at
http://library.wolfram.com/conferences/devconf99/#programming.
At that conference Micheal Trott and Adam Strzebonski of Wolfram  Research are mentioned as the inventors of this trick : see  "Trott-Strzebonski method for In-Place Evaluation".

The next cell creates a held expression and evaluates all sub expressions  with the Head Plus but nothing else evaluates. In this example  Erf[∞]+5, 1+3, 5+4, evaluate to 6, 4, 9, 2 respectively since they each have the head Plus.

ReplaceAll at subexpressions matching a pattern

Comments about the code above:



(1)  (Reverse@Sort...) ensures the later positions in (posn) still apply even
after folding several times.



(2)  We can't let the replacements take effect until after Fold is done
because the

     position of things may change.  To solve that I use Hold@@{expr}.



(3)  I don't use (rules__?OptionQ) because OptionQ[n_Integer->n+2] returns
False.

     Instead I use  rules:(_Rule|_RuleDelayed).



Now some examples:

Here we use the rule (b→bb) on all subexpressions with the head
Gamma.

Here we use the rule (n_Integer→n+2) on all subexpressions with the
head Gamma.

Here we use both rules (b→bb, n_Integer→n+2) on all
subexpressions with the head Gamma.

Making a periodic function

In the next cell I define a periodic function that's piecewize linear.  

This example would still work if I changed

f[x_?(Im[#]===0&)] := f[Mod[x,5]]   to simply

f[x_] := f[Mod[x,5]]   but by using the more complicated the definition (f)
is only used for real arguments.

However, Mathematica can't do things like Integrals or Fourier series on the function defined  above.  If you wanted to find several terms of the Fourier series of the  above function I recommend defining the function over the period (-2.5  ≤ x ≤ 2.5) using UnitStep functions as in the next  cell.  Of course the function isn't actually periodic this way, but the  function FourierTrigSeries doesn't care about that.

Needs["Calculus`FourierTransform`"];

FourierTrigSeries[f[t], t, 4, FourierParameters→{0,5/2}]

In general you should find that Mathematica's full symbolic capabilities including Integrate, LaplaceTransform, etc.  can be used on piecewise continuous functions defined in terms of UnitsStep  functions. This isn't true for the programming style used earlier.

Take generalized

The #& notation is explained in the discussion of Function.

Portions of the function above are implemented below to illustrate how it
works.

The effects of FoldList and Through are hard to grasp.  The lines below
demonstrates what they do.

Thread generalized

This was written by Carl Woll.

Flatten without evaluating subexpressions

This was written by Dave Wagner.

Preventing TrigExpand from expanding as far as it can

You may not want TrigExpand to go as far as possible in expanding the
trigonometric expression below.  Tom Burton prevents TrigExpand from
expanding the products by changing the Head Times, and then changing it back
to Times after using TrigExpand.

Finding the Domain of a defined function

Suppose we have a function (f[_]) that is defined for certain arguments.

The Cells below include a  program designed to indicate the values where (f)
is defined.

Definitions for (f) are stored in DownValues[f].

A certain part of each DownValue indicates the type of argument for which (f)

is defined.  This operation is an important part of the function below which
gives the same

result without the Pattern matching notation.

  

Now Domain[f] can be use to determine the values where (f) is defined.

Finding the constant term(s) in a sum

The short program below will find the constant term in a sum.  The #&  notation used here is explained in the discussion of Function.

Below we see several examples where the constant term is found.  When the
expression has no constant term 0 is returned.

ConstantTerm isn't defined if it's given an argument that isn't a numeric
function.

An alternate definition for ConstantTerm with the same effect is given in the
next cell.

KroneckerProduct of matrices

InterpolatingFunction form for the solution to a system of equations.

Carl Woll gave the solution below on how to find Interpolating functions  that solve a set of equations.  
Suppose you are trying to find  Interpolating functions for y[t], x[t] that solve the equations
  for (t) between 0 and 1.

soln=NDSolve[eqs,{x,y},{t,0,1}]

x1=x/.soln[[1]];

y1=y/.soln[[1]];



Block[{$DisplayFunction=Identity},

    p1=Plot[Re[x1[t]],{t,0,1},PlotLabel->"Real Part of x"];

    p2=Plot[Im[x1[t]],{t,0,1},PlotLabel->"Imaginary Part of x"];

    p3=Plot[Re[y1[t]],{t,0,1},PlotLabel->"Real Part of y"];

    p4=Plot[Im[y1[t]],{t,0,1},PlotLabel->"Imaginary Part of y"];

    ];

Show[GraphicsArray[{{p1,p2},{p3,p4}},GraphicsSpacing ->
0.2,ImageSize->{550,380}]];

Making a list of integers relatively prime to n

Suppose we want to find a list of integers relatively prime to a Positive
Integer (n). An immediate candidate is the function in the next cell, but it
isn't very fast for large n because it has to examine each integer between 2
and (n-1).

Ranko Bojanic found a faster solution that uses the fact  that  k< n is
relatively prime to n if it is not divisible by any prime divisor of n. The
list of  prime divisors of n is obtained by the following function.

The definition for (pdList) in the next cell is faster, but it doesn't  work in Version 3 or earlier. I am sort of splitting hairs here because an  expression has to be quite large for Part[expr, All, 1] to be significantly  faster, and FactorInteger will seldom return a list with more than factors.

If (d) is a prime divisor of (n), then a list of integers between 2 and (n-1)
which are not divisible by n is given by the following line.

We have to take all these lists, for all prime divisors of n, and their
intersection will be the list of all integers between 2 and (n-1) which are
relatively prime to n.  It is now easy to see that the list of relative prime
integers of n can be found from the following function.

The solution in the next cell is a slight variation of a program Alan  Hayes wrote and it's a little faster than the last program. Keep in mind this  variation doesn't work in Mathematica Version 3 or earlier because of the use of Part[expt, All, 1].

How the above implementation of RelativePrimes works

An Algebraic Transformation

A user in the MathGroup wanted Mathematica to change

((-1+a) x-a y)\^2  into  ((a-1)x+a y)\^2

Allan Hayes gave the solution in the next cell.  This solution has to be
entered on a case by case basis.

FactorRule in the next cell works and the rule doesn't need to know what  variables are involved.  The only disadvantage of FactorRule is that it's  more complicated than the first solution.  The use of (expr :) is discussed  in my section on Pattern.

A New Together

A user wrote to the MathGroup indicating that Together takes a very long  time with expressions that have on the order of leaves. They noticed together expands the numerator of the result as in  the next example, and suspected time could be saved if the numerator wasn't  expanded.

Allan Hayes wrote SimpleTogether below which doesn't expand the numerator.
Allan's code is given in the next cell and it's very fast and very slick.

In the next cell we see that SimpleTogether doesn't expand the numerator of
the result.

Notice Allan's program uses Union which has a SameTest option.  Some  nuances of Union and it's option are discussed at:
http://support.wolfram.com/Kernel/Symbols/System/Union.html  where it says the default SameTest setting used by Union is stronger  than using Equal or SameQ! This is demonstrated in the next line.

Next Union returns only one of the numbers when SameQ is used for SameTest.

Considering the lines above I think Alan's SimpleTogether function should use
Union with the option SameTest→(SameQ[{#1,#2}]&).  In some applications
something else might be needed.  With that in mind I wrote a version of
Allan's program that has a SameTest option. I also give SimpleTogether the
options Modulus and Trig which the built-in version has and I give it a usage
message. The only thing missing is the extension option which the built-in
version of Together does have.

Next we see that SimpleTogether still works.

I don't demonstrate the options but the default settings are shown below.

Making a Tensor into a Matrix

Consider the tensor (t1) in the next cell.

A user in the MathGroup wanted to express this tensor as a matrix with
MatrixForm in the next cell. Visually  this is simple. However, the
complicated structure of a tensor makes it a challenge to write an elegant
program that makes the conversion.

Allan Hayes provided the following two solutions in the MathGroup.

Now suppose you want to do the same thing on the higher rank tensor (t1)
below

Mathematica Version 4 has a NestWhile function that comes in handy here. The
code in the next cell will merge together tensors of any rank.

Distribute - A slick application

Dr. John Erb sent a problem to the MathGroup.
He had several pieces of  plastic of different thickness and wanted elegant Mathematica code that would determine what thicknesses can be made by stacking  together one or more of the pieces of plastic.

Robby Villegas replied to  the MathGroup in   [mg3203] Re: array ordered, of 18 Feb 1996:
This problem  amounts to finding all subsets of the list of thicknesses T, and for each  subset adding up its elements. Any subset of T can be described by giving a  status to each of T's elements: absent or present (as noted in Jorma  Virtamo's solution). In terms of a contribution to the total thickness, the  ith element adds 0 if it is absent, or adds (ti) if it is present. Thus, if  you take the Cartesian product of these ordered pairs: {0, t1} x {0, t2} x .  . . x {0, tn} you get all possible combinations of plates, e.g. {t1, 0, 0,  t4, t5, ...}, and you can add the elements of each combination. 'Distribute'  is perfect for forming Cartesian products. .....
What follows is a slight  variation of Robby's solution.
First we get a list of thicknesses paired  with zero in s2.

Next Distribute returns the cartesian product of all the ordered pairs.

As Robby Villegas pointed out Distribute can take third and fourth arguments
which specify heads that should be used to replace the outer and inner heads
that were distributed. In this case we want the outer head to remain as List
and we want to add the sublists. Hence Distribute in the next line does the
job.

Alternatively we can use the next cell to get lists of the form {{sum,
list},...} where 'sum' is the total of the elements in 'list', and Union has
removed elements with same sum, and sorted the list so the values of sum are
increasing.

Union without sorting

In June 2000 there was a long thread in the MathGroup and you can read  about the discussion at
http://library.wolfram.com/mathgroup/archive/2000/Jun/msg00115. html  on how to write an efficient function that does the same thing as Union  without sorting the elements. Many agree that the best way to do this is with  the following ingeniuos function written by Carl Woll.

DeleteRepititions is demonstrated below.

Actually the version Carl Woll posted used Block where I use Module above.  
The function runs a bit faster when defined using Block, but then it gives
the wrong result in the following example.

Created by Mathematica  (May 17, 2004)