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.

x |

f[x, 2] |

f[f[x, 2], 3] |

f[f[f[x, 2], 3], 4] |

f[f[f[f[x, 2], 3], 4], 1] |

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].

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)