Tutorial 06
Parametric Plot, Modular Programming, and Symbolic Equation Solving

Copyright © 2003 by Bill Titus, Carleton College,
Department of Physics and Astronomy, Northfield, MN  55057-4025
March 1, 2003

Part 1 - Getting Started

Step 1 - Start up Mathematica and open the notebook  Tutorial06.nb.  Remember, this notebook is located in the workshop folder.

Step 2 - OPTIONAL.  Since you're going to edit this notebook, you need to make a copy of this notebook and work with that copy.  To do so, select  File / SaveAs  from the Mathematica menu and save the notebook in the workshop folder under the name  myTutorial06.nb.

Part 2 - Creating a Parametric Plot

Topics and Skills

1.  The Mathematica command  ParametricPlot[]
2.  The option  AspectRatio -> Automatic
3.  The Mathematica functions   Sqrt[]  and  ArcTan[]

Step 1.  Sometimes it's useful to plot curves in terms of their cartesian coordinates using a parameter.  For example, a circle of radius  r  can be described parametrically by the two equations  x[r, ϕ] = r cos[ϕ]  and  y[r, ϕ]  = r sin[ϕ],  where  x  and  y  are the cartesian coordinates of a point on a circle, while  ϕ   is the angle (in radians) that the point makes with respect to the x-axis.   Here  ϕ  is positive for counterclockwise rotations and negative for clockwise ones.  To plot a complete circle, you will need to pick  ϕ  so that it goes around once, e.g.  0 ≤ ϕ < 2π.

Suppose you know  r  and  ϕ  and want to calculate the corresponding values for  x  and  y .   Create two functions,  x[r, ϕ]  and  y[r, ϕ]  in the cell below that define the appropriate functions and verify that your functions work by picking   r  =  2.0  and  ϕ  =  0.35.  Remember to clear all the variables you use as well as the name of the functions.

Comment 1

1.  As soon as I define a function, I try it out numerically.  When I do this, many times I discover a simple error that can easily be corrected.  It's much better to find errors quickly than later on when you've used the function more extensively.

Step 2.  Mathematica allows you to plot functions that are parametrically defined in 2D by invoking the Mathematica function  ParametricPlot[].  Get extended help on this command in the cell below.

Comment 2

1.  Make sure you understand the usage statement for this function.  For example, to plot a circle of radius  r,  you could type  

			ParametricPlot[{x[r,ϕ],y[r,ϕ]},{ϕ,0,2π}] .
			
2.  Many of the options for  ParametricPlot[]  are the same as those for  Plot[]  and  ListPlot[].

Step 3.  Produce a frame plot of a circle with radius  2.0  in the cell below by using   ParametricPlot[]  and the equations that you defined for  x[r, ϕ]  and  y[r, ϕ]  in Step 1.

Comment 3

1.  Your circle really doesn't look like one because of Mathematica's fixation on picking an aspect ratio (the ratio of height to width) equal to the golden ratio  (1/2) ((Sqrt[5]) + 1).  One way to remedy the situation is to use the option  AspectRatio -> Automatic,  which finds the aspect ratio using the actual coordinate values calculated.  Thus a circle should look like a circle.  Add this option to your code above and see what happens.

Step 4.  As "more" practice, produce a frame plot of a circle of radius  2.0,  whose center is located along the x-axis at  x = 3.0.  If you think a bit, you can use your equations for  x[r, ϕ]  and  y[r, ϕ]  defined in Step 1 with a slight change to the first argument in  ParametricPlot[].  Choose a  PlotRange  option so that the origin of your coordinate system is included in the plot.  Remember, copying, pasting, and editing can save time.

Step 5.  Review the following information on the fundamental strain equations in 2D given in the cell below.  For more detailed information, see Structural Geology of Rocks and Regions, 2nd edition, by George H. Davis and Stephen J. Reynolds (John Wiley, New York, 1996, pages 68-69).

The fundamental strain equations allow you to determine the stretch  S  and shear strain  γ  for any line in a strained body if the angle  (Subscript[θ, d])  that line makes with respect to the maximum stress direction is known as well as the stretches  (Subscript[S, 1]) and  (Subscript[S, 3])  along the  (Subscript[S, 1])  and  (Subscript[S, 3])  direction, respectively.  The sign convention for  (Subscript[θ, d])  is such that  (Subscript[θ, d])  is positive for clockwise rotations from the  (Subscript[S, 1]) direction and negative for counterclockwise rotations; the range for  (Subscript[θ, d])  is between  -π/2  to  π/2.  The shear strain  γ  is related to the angular shear by  tan[γ].  

In particular, if  the quadratic elongations are defined by   λ = (S^2),   (Subscript[λ, 1]) = ((Subscript[S, 1])^2),  (Subscript[λ, 3]) = ((Subscript[S, 3])^2), then  λ  and  γ  satisfy the equations

		(1/λ)= (1/2) ((1/(Subscript[λ, 3])) + 
              (1/(Subscript[λ, 1]))) - (1/2) ((1/(Subscript[λ, 3])) - 
              (1/(Subscript[λ, 1]))) cos[2  (Subscript[θ, d])],  		  (γ/λ)=(1/2) ((1/λ)3 - (1/λ)1) 
            sin[2  (Subscript[θ, d])].

Step 6.  Assume you know  (Subscript[S, 1]),  (Subscript[S, 3]),  and  (Subscript[θ, d]).  Define  α   to be  1/λ   and   β   to be  γ/λ .  Use the fundamental strain equations to create the functions  α[(Subscript[s, 1]), (Subscript[s, 3]), θd]   and  β[(Subscript[s, 1]), (Subscript[s, 3]), θd]  by activating the code below.   Verify that your functions work by picking reasonable numerical values for  (Subscript[S, 1]),  (Subscript[S, 3]),  and  (Subscript[θ, d]).

Clear[α, β, s1, s3, θd]
α[s1_, 
          s3_, θd_] = (1/2) ((1/(s3^2)) + 
                (1/(s1^2))) - (1/2) ((1/(s3^2)) - (1/(s1^2))) 
            Cos[2*θd];
β[s1_, s3_, θd_] = (1/2) ((1/(s3^2)) - (1/(s1^2))) 
          Sin[2*θd];

Comment 6

1.  Note that you've used the Mathematica functions  Sin[]  and  Cos[]  in this function definition.

Step 7.  Davis and Reynolds (pages 69-70) have an example where they find the stretch  S  and the angular shear  ψ  for a line when  (Subscript[S, 1])= 1.55,  (Subscript[S, 3])= 0.65,  and   (Subscript[θ, d])= -26.5°.  In particular, they determine that  S  is  1.1  and  ψ  is  -43°.  You'll shortly use your equations for  α[(Subscript[s, 1]), (Subscript[s, 3]), θd]  and  β[(Subscript[s, 1]), (Subscript[s, 3]), θd]  to verify that the author's have done the calculations almost correctly.

Your calculation will require some thought since to find  S  and  ψ,  you don't just evaluate  α  and  β.   In particular, with a little algebra, you can show that  S = 1/(Sqrt[α])  and  ψ = arctan[β/α].   Also recall that the conversion from radians to degrees requires a multiplication by  180/π  while from degrees to radians requires the reciprocal.   Finally you'll need to use two new  Mathematica functions,  Sqrt[]  and  ArcTan[],  whose purpose should be self evident.  

Look over the code in the cell below, make sure you understand what it's doing, and then active the cell to find the numerical values for  S  and  ψ.

Clear[s1, s3, θd, θdRadians, 
      s, ψRadians, ψ]
s1 = 1.55; s3 = 0.65; θd = -26.5;
θdRadians = (π/
              180)*θd;
s = 
        1/Sqrt[α[s1, 
              s3, θdRadians]];
ψRadians = 
        ArcTan[β[s1, s3, θdRadians]/α[s1, 
              s3, θdRadians]];
ψ = (180/π)*ψRadians;
{s, ψ}

Comment 7

1.  You should find, to two significance figures, that   S  is  1.1  and   ψ = -44;  these are basically the same values as those found by Davis and Reynolds.
2.  You could have written the above code much more compactly as follows, although it would have been more difficult to understand.  Sometimes clarity is better than brevity, especially if you haven't used your code for awhile, or if someone else is trying to understand what you did.

Clear[s1, s3, θd]
s1 = 1.55; s3 = 0.65; θd = -26.5;
{1/Sqrt[α[s1, s3, (π/180)*θd]], (180/π)*
        ArcTan[β[s1, s3, (π/180)*θd]/α[s1, 
              s3, (π/180)*θd]]}

Step 8.  Suppose  (Subscript[S, 1]) = 1.55  and  (Subscript[S, 3]) = 0.65.  Plot the Mohr circle strain diagram,  α  versus  β,  for this case and name it plot1.  Use the plot range  -1.0  ≤  x  ≤  3.0  and  -2.0  ≤  y  ≤  2.0.

Step 9.  It would be nice to show the location of the center of the Mohr circle.  This point has cartesian coordinates  {α[(Subscript[s, 1]), (Subscript[s, 3]), π/4],0}.   One way to do this is to use   ListPlot[]  to plot the single point corresponding to the center of the circle, superimpose that plot and   plot1  from Step 8, and then show both plots together.  Do that work in the cell below.  Be sure you make the point representing the center large enough so you can see it.  Also remember that the point's list must reside inside a list in order for  ListPlot[]  to operate correctly.

Part 3 - Modular Programming and the Fundamental Strain Equation in 2D

Topics and Skills

1.  The Mathematica command   Module[]  and its structure

Step 1.  In Step 6 above, you found the numerical values for  S  and  ψ   from your functions  α  and  β,  assuming you knew  (Subscript[S, 1]),  (Subscript[S, 3]),  and  (Subscript[θ, d]),  where  (Subscript[θ, d])  was expressed in degrees.  However, your calculation was not a direct one.  You had to convert  (Subscript[θ, d])  to radians,  find   S  and  ψ  from  α  and  β,  and finally convert  ψ  to degrees.   If you're doing a number of calculations with different values for  (Subscript[S, 1]),  (Subscript[S, 3]),  and  (Subscript[θ, d]),  it would be nice to put all that code in a single function call.  Mathematica allows you to do that using the function   Module[]   that wraps itself about a section of code to produce a single function.  Get help on  Module[].

Step 2.  OK, that wasn't helpful.  The code in the following cell defines a function   strain[s1,s3,θd]  that does what you want.   Look over the code quickly, activate the cell, and try out the function on the parameters  (Subscript[S, 1]) = 1.55,  (Subscript[S, 3]) = 0.65,  and  (Subscript[θ, d]) = -26.5.

Clear[strain, s1, s3, θd]
strain[s1_, 
          s3_, θd_] := 
Module[
{θRadians, 
            s, ψRadians, ψ}, 

θdRadians = (π/180)*θd; 
s = 
            1/Sqrt[α[s1, 
                  s3, θdRadians]]; 
ψRadians = 
            ArcTan[β[s1, s3, θdRadians]/α[s1, 
                  s3, θdRadians]]; 
ψ = (180/π)*ψRadians; 

{s, ψ}
];

Comment 2

1.  The structure of   Module[]  is not difficult to understand.  First, in defining a function with  Module[],  you must use  :=   between the name of your function and its arguments and  Module[],  i.e.,  you need to use  SetDelayed[].
2.  The first line inside  Module[]  must consists of a list of all the variables that appear in the code (the internal or local variables), except for those in the function call (the external variables).  So in your case,  s1,  s2,  and  θd  are the external variables, while  θRadians,  s,  ψRadians,  and  ψ  are the local variables.  Local variables don't exist outside of   Module[]. 
3.  The list of local variables must also be present, even if it's empty; a comma always follows this list.
4.  Next comes your code, with each individual line separated by a semicolon.
5.  Finally, the last line of the module, which has no semicolon after it, is what the function returns.  In your case, the function returns the list  {S, ψ}.

Step 3.   Suppose you know that  (Subscript[S, 1]) = 1.55  and  (Subscript[S, 3]) = 0.65,  and you want to find a table of values of  {S, ψ}  corresponding to values of   (Subscript[θ, d])  going from  -90°  to  +90°  in steps of   10°.   Use   Table[]  and   strain[]  to generate this list.   Call the list  strainList,   display the list in table form, and plot the list so that  S  is on the horizontal axis and  ψ  is on the vertical axis.

Part 4 - Symbolic Equation Solving and the Fundamental Strain Equation in 2D

Topics and Skills

1.  The Mathematica commands  Solve[],  NSolve[],  FindRoot[]
2.  Equations in Mathematica and the use of double equal signs, i.e.,   ==
3.  Replacement rules as solutions to equations
4.  Double bracket notation  [[]]  and the equivalent function  Part[]
5.  The Mathematica function  Chop[]

Step 1.   Davis and Reynolds (pages 74 - 76) have any interesting example of the use of the fundamental strain equations to determine stretches  (Subscript[S, 1])  and  (Subscript[S, 3])  of a finite strain ellipse from measurements of the stretch  S,  the angular shear  ψ,  and the displacement angle  (Subscript[θ, d]) of the hinge lines on two brachiopods.   Suppose Cam, the geologist, has his structural geology students measure  S  and  ψ  corresponding to four hinge lines  (Subscript[H, E]),  (Subscript[B, E]),  (Subscript[H, M]),  and  (Subscript[B, M]) as indicated in Figure 2.49 in the Davis and Reynolds text.  In particular, the students generate the following data:

Line		S (cm)		ψ (°)		(Subscript[θ, d])(°)

		(Subscript[H, E])		0.72		-32		-50
		(Subscript[B, E])		1.58		+32		+10
		(Subscript[H, M])		1.10		+41		+30
		(Subscript[B, M])		1.27		-45		-15
			
Create a list called  strainResults  that corresponds to the data  {S, ψ, (Subscript[θ, d])}  given above.

Step 2.   For each element in  strainResults,  you should, in principle, be able to find the corresponding values for  (Subscript[S, 1])  and  (Subscript[S, 3]),  that is, given  S,  ψ,  and  (Subscript[θ, d]),  you would solve the mathematical equations  {s,ψ} = strain[(Subscript[s, 1]), (Subscript[s, 3]) , (Subscript[θ, d])]  for  (Subscript[S, 1])  and  (Subscript[S, 3]) (remember,  strain[(Subscript[s, 1]), (Subscript[s, 3]) , (Subscript[θ, d])]  produces the list  {S, ψ} ).   You're going to have Mathematica find the solution symbolically by using the function  Solve[].  First get help on this function.

Comment 2

1.  For your information, the Mathematica functions that allow equations to be solved numerically are  NSolve[]  and  FindRoot[].   These are valuable functions when you have to solve equations numerically.

Step 3.   In Mathematica, equations are expressed using double equal signs and evaluate to something that is either true or false.  You can see this in action by evaluating the cell below:

Clear[a, b]
{5 == 3,  Sin[π/2] == 1,  a == b}

Comment 3

1.  If Mathematica can't evaluate an input, it will usually echo the input back to the screen.  In the case above, since  a  and  b  have no assigned values (they are just symbols at this stage), Mathematica can't say whether  a = b  or  a ≠ b.

Step 4.  As practice, use  Solve[]  to symbolically find the values of  x  that satisfy the quadratic equation  a (x^2) + b x + c = 0  where  a,   b,  and  c  are arbitrary constants.  Remember, start your code by clearing all the variables you're going to use.  You don't have to do this, but it's good practice to do so.

Comment 4

1.  A common error in working with equations is to use a single equal sign instead of a double one.  Go back to Step 4 and see what happens if you use  =  instead of  ==.   When you're done, return the equation to its original form. 
2.  Mathematica has outputted a list consisting of two lists.  Inside each of those lists is a solution expressed by a replacement rule.  In reading replacement rule, it's helpful to think of the symbol  ->  as "goes to".

Step 5.  Activate the code below that solves the set of simultaneous equations   a x + b y  =  u   and  c x + d y  =  w   for  x  and  y,  given constants  a,  b,  c,  d,  u,  and  w.

Clear[a, b, c, d, u, w]
Solve[
{
a*x + b*y == 
          u, 
c*x + d*y == w
}, {x, 
        y}]

Comment 5

1.  If you have more than one equation to solve, you need to put the equations in a single list, with each equation separated by a comma.
2.  The code in  Solve[]  has been formatted with line breaks to make its structure clear.  You don't have to do this.

Step 6.  See if Mathematica can symbolically solve the strain equations  for  (Subscript[S, 1])  and  (Subscript[S, 3]),  given  S,   (Subscript[θ, d]),  and  ψ.  To do so, just activate the cell below.

Clear[soln, s, ψ, s1, s3, θd]
soln = 
      Solve[
{
{s, ψ} == 
            strain[s1, s3, θd]
}, {s1, 
          s3}]

Comment 6

1.  There are really two simultaneous equations in the code above:  one for  s  and one for  ψ  since the strain function outputs a list.
2.  Note that the output has been given the name  soln  so you can refer to it later on.
3.  It looks like there are four possible pairs of solutions.  However, by definition  S,  (Subscript[S, 1]),  and  (Subscript[S, 3])  are positive and these conditions eliminate the first 3 solutions.

Step 7.  Let's see how to create two new functions  ss1[s,ψ,θd]  and  ss3[s,ψ,θd]  that correspond to the fourth solution given in Step 6.  As a first step, enter the cell below and use its output to try to understand the purpose of the double bracket notation  [[ ]],  which is equivalent to the function  Part[].

soln[[4]]
soln[[4, 1]]
soln[[4, 1, 2]]

Comment 7

1.  Every expression is Mathematica is made up of different parts.  For lists, each part corresponds to the elements of the list.  For rules, the first element is the symbol while the second element is the expression.  If you want to delve into this formalism in more detail, check out the functions  Part[]  and  FullForm[]  in the help browser.

Step 8.  Now use what you've learned in Step 7 to define the functions  ss1[s,θd,ψ]  and  ss3[s,θd,ψ].  See if your functions are working by using the parameters   S = 0.72,  ψ = -32,  and  (Subscript[θ, d]) = -50.

Comment 8

1.  For output expressed in the form  {ss1,ss3},  you should find the answer  {1.44474+0.i, 0.591267+0.i}, which look a bit strange.  The  0.i  that appears means that there was a very small imaginary part to your answer, probably due to round-off errors.   Here  i = (Sqrt[-1]).

Step 9.  To eliminate a very small contribution to a numerical answer, you can invoke the Mathematica function  Chop[]  as follows:

Clear[s, ψ, θd]
s = 
      0.72; ψ = -32; θd = -50;
{ss1[s, θd, ψ], ss3[s, θd, ψ]} // 
      Chop

Comment 9

1.  Remember,  the  //  notation is the postfix version of the function that is applied after the argument of the function is given.  In this case, the postfix notation applies  Chop[] to the preceding list and is equivalent to the nested commands   Chop[{ss1[s,θd,ψ], ss3[s,θd,ψ]}].

Step 10.  Recall that your original data list  strainResults  consists of four elements, each of the form {s,θd,ψ}.  You can easily transform   strainResults  into a list with elements  {s,θd,ψ,s1,s3}  by using the transformation technique discussed in  Tutorial04.  In particular, look at the code below, figure out what it's doing, and active the cell to produce the desired result.

Clear[newResults, x, y, z]
strainResults
newResults = 
      strainResults /. {x_, y_, z_} -> {x, y, z, ss1[x, y, z], 
              ss3[x, y, z]} // Chop

Step 11.  Display  newResults  in tabular form.

Comment 11

1.  Note that the values of  (Subscript[S, 1])  and fairly consistent, although that's not the case for  (Subscript[S, 3]).

Part 5 - Save Your Notebook (OPTIONAL)

Step 1.  Remove all output cells from this notebook, and to do one final save of your notebook before you quit Mathematica.