Example 09
Iterative Equations and 
Modular Programming

Copyright „ 1993 by Bill Titus, Carleton College, 
Department of Physics and Astronomy, Northfield, MN  55057-4025
September 6, 1993

Topics and Skills

1.  The  For  loop and its structure.
2.  Print[]. 
3.  ListPlot[], the option  PlotJoined, and the graphics directive   
     PointSize[].
4.  Comments using  (* ... *).
5.  Writing a function using  Module[]. 
6.  Aborting a calculation.

Problem

Most force functions in classical mechanics, because of their complicated time dependence or their explicit dependence on velocity and/or position, can't be integrated directly to obtain closed-form expressions.  However, such problems can always be handled numerically using an integration algorithm such as Euler's method, often used in introductory physics classes, or the more sophisticated Runge-Kutta method.

Recall that in Euler's method, one breaks up the time domain into  n  equally spaced intervals of size  dt:

to = t0 < t1 < t2 <  ...  < tn-1 < tn = tf
dt = (tf - to)/n

For notational purposes, let

x[ti] = xi, v[ti] = vi, a[x[ti],v[ti],ti] = ai.

Euler's iterative equations then become

ti  =  ti-1 + dt
	   vi  =  vi-1 + ai-1 dt
	   xi  =  xi-1 + vi-1 dt

These equations allow you to "bootstrap" up to  xi  and  vi  from the initial conditions:   x0 = x[t0]  and   v0 = v[t0].

Solution - Using a Sequence of Iterative Equations

Step 1 - Briefly describe how you would implement Euler's algorithm if you were doing linear programming using arrays.

Comment 1

1.  You might do the following:
     a.  Input information:  t0,  tf,  x0,  v0,  and  n.
     b.  Initialize the arrays.
     c.  Iterate the equations using a loop structure like  For,  Do,  or   While.
     d.  Output the data in tabular and/or graphical form.

Step 2 - You'll develop the algorithm symbolically to begin so you can make sure it's working correctly.  In particular, look over the following code and then activate the cell.

Clear[t, a, v, x]
	Clear[to, tf, dt, n, vo, xo, i]

	n = 4;
	t[0] = to;
	v[0] = vo;
	x[0] = xo;

	For[i = 1, i <= n, i++,
	    t[i] = t[i-1] + dt;
	    v[i] = v[i-1] + a[i-1] dt;
	    x[i] = x[i-1] + v[i-1] dt;
	    Print[t[i]];
	    Print["  ", v[i]];
	    Print["  ", x[i]];
	    Print["  "];
	   ];

Comment 2

1.  Does the output look correct?
2.  Some people tend to clear every possible symbol.  It's overkill, but safe.
3.  The code doesn't determine  dt  nor does it have an explicit expression for the acceleration.  You'll do that later.
4.  Be sure to use square brackets for the arrays and not round parentheses.  What would happen if you used round parentheses?

Step 3 -  Get help on  For  and summarize what you find.

Comment 3

1.  Note that the four arguments of  For  are separated by commas.  Each one can be composed of composite statements, with the composite statements separated by semicolons.

Step 4 - Get help on  Print[].  What do you find?

Comment 4

1.  Print[]  is useful for debugging a sequence of code.

Step 5 - Copy the code from Step 2 to a cell below and complete the algorithm by adding in the following lines at appropriate places;  keep  n equal to  4  for a while.  Run your program to see if your numerical output looks reasonable.

to = 0.0;
	    tf = 10.0;
	    vo = 0.0;
	    xo = 1.0;
	    dt = (tf - to) / n;
	    a[i-1] = -x[i-1];
	    a[n] = -x[n];

Comment 5

1.  What physical problem are you modeling?
2.  Why use reals for initial conditions?  What happens if you didn't?  Try it.
3.  Why include the line  a[n] = -x[n]?

Step 6 - Use  Table[]  to create a table of  {t, x}  values.  Give your table the name   tempList   so you can refer to it later.

Step 7 - It would be nice to plot the  {t, x}  data.  What MMA command might plot a list of data?   Find the command and describe what it does.

Comment 7

1.  What are the two versions of  ListPlot[]?  Which one do you need?
2.  Do the options look familiar?
3.  What do you think the option  PlotJoined  does?  Find out.
4.  If you set  PlotJoined -> True, the individual points are not plotted.  How could you get both the points and the lines connecting the points to appear on the same graph?
5.  A useful graphics directive for  PlotStyle  is  PointSize[].  What does this directive do?

Step 8 -  Use  ListPlot[]  to plot your  {t, x}  points on a graph, but don't connect them by lines.  Use the option

    PlotStyle ->
     {RGBColor[1, 0, 0], PointSize[0.02]}
    
to color your points red with point size of  0.02.  Also be sure to label your axes.

Step 9 - Copy your program from Step 5 to a cell below and insert code that allows you to displays two plots:  x  versus  t   and   v  versus  t.  Eliminate the  Print[] commands from your code.  Use a name for the tables so you can refer to them in the  ListPlot[]  command, but don't print out the tables.  Keep  n   at  4  until your code is working correctly; then switch to  n  equal to  100.

Comment 9

1.  Be sure to add the name of your temporary table(s) to the  Clear[]  statement.
2.  Do you notice something wrong with the output? 
3.  Modify the iterative equations to correspond to Euler's Last Point Method where   xi -> xi-1 + vi dt.  Include a comment to this effect after the line of code where you made the change.  Comments in MMA are surrounded by  (* ... *).  
4.  Finally, decrease the point size to  0.01  and then rerun your code to see if the output looks more reasonable.

Solution - Encapsulating the Euler Algorithm and Plotting Routine to a Single Function Using  Module[]

Step 1 - Get help on the MMA function  Module[].

?Module

Comment 1

1.  Module[]  can be used to write what would be called a procedure in Pascal, a subroutine in Fortran, or a function in C.
2.  The {x, y,  ...} are local variables and only exist within the module.  This is in contrast to global variables which are visible outside the module.  Local variables need not to be cleared.
3.  If you set a function equal to a module, you need to use  SetDelayed[]  and then place semicolons after each line in the module (although the last line doesn't need one).

Step 2 - Look over the following code which creates a function which plots the logistic equation and prints out a table of values.  Identify the local and global variables.  Then activate the cell and try out the function with  xo = 0.4,  r = 0.78, and  n = 20.

Clear[xo, r, n, logisticPlot]
logisticPlot[xo_, r_, n_] :=
  Module[
   {i, x, tempList},
   
    x[0] = xo; 
    
    For[i = 1, i <= n, i++,
        x[i] = 4*r*x[i-1]*(1 - x[i-1])];
       
    tempList = Table[x[i], {i, 0, n}];   
    ListPlot[
          tempList,
          PlotStyle  ->
              {RGBColor[1, 0, 0], PointSize[0.02]},
          AxesLabel  -> {"i", "x[i]"},
          PlotRange  -> {0, 1},
          AxesOrigin -> {0, 0}
        ];
        
    Table[{i, x[i]}, {i, 0, n}] // TableForm
  ];

logisticPlot[0.4, 0.78, 20]

Comment 2

1.  What would happen if a semicolon was placed after the last line in the module, i.e., after  TableForm[]?
2.  Would the table of values had been printed out if the  Table[]  command had occurred before  ListPlot[]?

Step 3 -  Use  Module[]  to convert your sequence of code for the Euler proplem into a single function,

              eulerDEQ[to, tf, n, vo, xo]. 
  
Then call your function with different parameters to verify that it's working correctly.

Comment 3

1.  What do you see are the advantages and disadvantages for using modules?
2.  Indents and spacing are very important in making your code easy to read.

Step 4 -  When doing a loop calculation, you may want to abort the calculation either because it's taking too long or you did something to produce an infinite loop.  If this happens, there are several ways to abort, ranging from gentle to harsh:

On the Next:

1.  Select   Action \ Interrupt \ Interrupt Calculation, or <Command>-,
2.  Select   Action \ Interrupt \ Abort Calculation, or  <Command>-.
3.  Select   Action \ Kernels \ Quit/Disconnect Kernel.
4.  Go to the Workspace menu and select  Tools \ Processes.  When the process panel comes up,  select  Applications  from the pop up menu at the top of the panel.  Highlight  Mathematica,  press the Kill button, and then press OK.  This technique looses all the work you haven't saved but sometimes it's necessary if the ominous wheel appears and won't go away.

On the Mac:

1.  Select   Action  \  Interrupt Calculation, or  <Command>-,
2.  Select   Action  \  Abort Calculation, or  <Command>-.
3.  Select   Action  \  Quit/Disconnect Kernel.

Now, activate the following cell and then abort.

For[i = 1, i < 1000000, i++]

Save Your Notebook

Remove any output cells from this notebook and then use  Save As  to store your notebook under the name  myEx9.