LESSON 6
Useful Mathematica  Techniques

Author's notice

This tutorial was written by 

Fred Goodman
Department of Mathematics
University of Iowa
Iowa City, IA 52242
goodman@math.uiowa.edu

This material may be freely copied and used for non-commercial purposes.  Please retain this author's notice, and do not distribute any altered version of this material without the author's permission  (unless  the material is altered beyond recognition).

Parametric Plots

Sometimes curves are given naturally in parametric form,  x = x(t), y = y(t).  For example, the unit circle is given as y = sin(t), x = cos(t).  Here is the Mathematica technique for plotting a curve given in parametric form:

ParametricPlot[{Cos[t], Sin[t]}, {t, 0, π}, 
        AspectRatio -> Automatic];

Example (Exercise 2.3.32)

In Exercise 2.3.32,  you find the x and y velocities of a projectile moving under the influence of gravity and air resistance to be the following, where 
u 	is the initial velocity, 
A 	the angle at which the projectile is launched, 
r 	the coefficient of air resistance, and 
g 	gravitational acceleration near the surface of the earth :

Clear[v,  w,  g, r,  u, A, t];

    
v[t_] := u Cos[A] (E^(-r t));
w[
          t_] := -(g/r) + (u Sin[A] + 
                (g/r)) (E^(-r t));

Consequently, the x and y position of the projectile are given by the following, where h is the initial elevation of the projectile (and the initial x - coordinate of the projectile is zero):

Clear[x,  y,  h];

    
x[t_]   =  (∫_0^t) v[
            s] \[DifferentialD]s
y[t_]   =   h  +  (∫_0^t) w[
              s] \[DifferentialD]s

We assign values to all the parameters r,  h,  g,  u,  A, and plot the trajectory of the projectile:

r  =   .2; 

    h  =  3; 

    g  =  32;

    u  =  125;

    A  =  Pi/4;

    ParametricPlot[{x[t],  y[t]},  {t,  0, 7} ]

Clear[r, h, g,  
      u,  A,  x,  y,  v,  w]

Mathematica procedures (Modules)

Let's go back to our example with the projectile, where we assign  values to r, h, and g but leave u and A as variables.

Clear[r,  h,  g,  v,  w,  x,  y,  u,  
        A];

v[t_] := u Cos[A] (E^(-r t));
w[
          t_] := -(g/r) + (u Sin[A] + 
                (g/r)) (E^(-r t));
x[t_]   =  (∫_0^t) v[s] \[DifferentialD]s
\(y[t_]\ \  = \ 
      Expand[\ h\  + \ ∫\_0\%t w[
                s] \[DifferentialD]s]
\)
r  =   .2; 
h  =  3; 
g  =  32;

For a given angle A, I would like to find the initial velocity u  which will make the projectile hit the top of a 10 foot high fence located 350 feet away. To do this, I solve for both t and u (and then just take the value of u. )    

The result is that, with an angle A of .8 radians, the projectile hits the top of the fence at time t = 5.566, if the initial velocity is u = 149.634:

A  =   .8;
soln  =  
      FindRoot[{x[t]  ==  350,  y[t]  ==  10},  {t,  
          10},  {u,  125}]
u /. soln

What I have above is essentially a 3 line Mathematica program, which inputs an angle A and outputs the corresponding initial velocity u  (which makes the projectile launched with angle A hit the top of the fence.)

I want to make a function out of the 3 line program which will allow me to input A and find u.

Here is how to do it.  In the following,  
Module[]  
is a Mathematica procedure.  The variable
soln
is a local variable in the Module, while
A, x, y, t, u
are global variables in the Mathematica session.

Clear[f,  Angle, u,  A];
f[Angle_]  :=  
      Module[{soln}, 

A  =  
          Angle; 
soln  =  
          FindRoot[{x[t]  ==  350,  y[t]  ==  10},  {t,  
              5},  {u,  150}]; 
u /. 
          soln
]

Let's try it out:

f[ .8]

We can now even plot u versus A:

Plot[f[A],  {A,  .6,   .7}]

Interpolating functions:  plotting the output of NDSolve

NDSolve gives as output a transformation rule involing an InterpolatingFunction, which is just a function given by a table of numerical data.  The first argument of  InterpolatingFunction (in the example below,  {-0.999999,  1.} ) is the domain of the numerical function.  

You can  plot the function only on this domain; if you try to plot in on a larger domain, you will get error messages and possible strange results.

The following example is from Exercise 2.2.21.  The command
Off[NDSolve::"ndsz"]
suppresses some warning messages about the singularity of the ODE.  It isn't necessary.

Off[NDSolve::"ndsz"]
Clear[x,  y,  f,  soln]

soln=NDSolve[{y'[x]==((1 + 3 (x^2))/(3 y[x]^2 - 6 y[x])),y[0] == 1},y,{x, -10, 10}]
f[x_] := 
      y[x] /. soln[[1]]

Let's attempt to plot the solution on the appropriate domain, and also on a larger domain, to see what happens:

Plot[f[x],  {x,  - .9999,   .9999}]

The following gives screwy results:

Plot[f[x],  {x,  -2,  2}]

Plotting two interpolating functions together

Let's solve the ODE in the example above, with two different initial conditions:

Off[NDSolve::"ndsz"];
Clear[f, x, y, soln1];
soln1=NDSolve[{y'[x]==((1 + 3 (x^2))/(3 y[x]^2 - 
                          6 y[x])),y[0] == 1},y,{x, -10, 10}]
f[
            x_] := y[
              x]/.soln1[[1]]

Clear[g, x, y,  soln2];

soln2=NDSolve[{y'[x]==((1 + 3 (x^2))/(3 y[x]^2 - 6 y[x])),y[0] ==  1/2},y,{x, -10, 10}]
g[x_] := 
      y[x] /. soln2[[1]]

The two solutions are defined on different domains.  How can we plot both of them together?
If we try to plot them both simultaneously on the union of their domains, we'll get nonsense:

Plot[{f[x],  g[x]},  {x,  -1.27,  
        1}]

Here is what to do instead, using the graphics command 
Show[]  
to show two plots together.

plot1  =  
      Plot[f[x],  {x,  - .999,   .999}]

    plot2  =  Plot[g[x],  {x,  -1.27,   .499}]

    Show[plot1,  plot2]