Reproducible mathematics

Symbolic mathematical computer languages, like MACSYMA, Maple, Mathcad, Mathematica, or Theorist, usually comprise a symbol-manipulating language, numerical routines, a general purpose programming language, graphics tools, and a user interface which incorporates text. In the following paragraph, I annotate the Mathematica script smile.ma, which derives the mathematical results outlined in the paper above.

The script starts by defining the ellipses (1) and (4) in eq1 and eq2:

eq1 = x^2 / a1^2 + z^2 / b1^2 == 1 ;
eq2 = (x-ds)^2 / a2^2 + z^2 / b2^2 == 1 ;

Next, I apply the tangency conditions. First, I eliminate z from eq1 and eq2 and yield eq3:

eq3 = Eliminate[{eq1,eq2},z][[3]] ;

I then define the parameter of eq1 and eq2 as independent of x, take the derivative of both equations and finally eliminate the slope p. The resulting equation eq4 corresponds to equation (8) in the article:

a1/: Dt[a1,x] = 0 ;
a2/: Dt[a2,x] = 0 ;
b1/: Dt[b1,x] = 0 ;
b2/: Dt[b2,x] = 0 ;
ds/: Dt[ds,x] = 0 ;
d1x = Dt[eq1,x] /. Dt[z,x] -> p ;
d2x = Dt[eq2,x] /. Dt[z,x] -> p ;
eq4 = Eliminate[{d1x ,d2x},{p,z}][[1,3]] ;

Combination of the two tangential conditions, eq3, eq4, eliminates x from the equation system, and yields equation eq5, which is represented by equation (9) in the article:

eq5 = eq3 /. Solve[eq4, x][[1]] ;

After introducing the definition (6) of b₂ in eq9, the program rearranges the terms by operating on both sides of any equation simultaneously. Finally, equation eq10 represents the desired quadratic polynomial (10). Mathematica then isolates the coefficients of the polynomial, which correspond to expression (11), (12), (13) in the paper above. They are the same coefficients implemented in the Ratfor program shct.rst:

den = (a2^2 b1^2 - a1^2 b2^2) ;
eq6 = Thread[Simplify[Plus[-(ds^2 + 2 a1^2 b2^2 ds^2 / den),eq5]], Equal];
eq7 = Thread[Simplify[Times[- den b2^2 / a2^2 ,eq6]],Equal] ;
eq8 = Thread[Simplify[eq7],Equal] ;
eq9 = Thread[ReplaceAll[eq8,b2 -> Sqrt[a2^2 - h2^2]],Equal] ;
eq10= Thread[Collect[ExpandAll[eq9],a2],Equal] ;
(* The polynomial *)
ex1 = Simplify[Coefficient[eq10[[2]],a2,0]] ;
ex2 = Simplify[Coefficient[eq10[[2]],a2,2]] ;
ex3 = Simplify[Coefficient[eq10[[2]],a2,4]] ;
eq11= 0 ==  ex3 as^2 + ex2 as + ex1 ;

The rest of the mathematica script smile.ma, which I will not discuss here, illustrates the shot continuation operator by computing its impulse responses for a series of cases (see Figure 5).

The script is executed by piping it into the mathematica start command 'math', or by calling it from within Mathematica. The semicolon behind an expressions prevents Mathematica from printing the value of this expression on the screen. Without the semicolon mathematica renders this line's intermediate result during execution.