(*^
::[ Information =
"This is a Mathematica Notebook file. It contains ASCII text, and can be
transferred by email, ftp, or other text-file transfer utility. It should
be read or edited using a copy of Mathematica or MathReader. If you
received this as email, use your mail application or copy/paste to save
everything from the line containing (*^ down to the line containing ^*)
into a plain text file. On some systems you may have to give the file a
name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
The line below identifies what version of Mathematica created this file,
but it can be opened using any other version as well.";
FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2";
MacintoshStandardFontEncoding;
fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8, 24, "Times";
fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times";
fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6, 14, "Times";
fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 18, "Times";
fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 14, "Times";
fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 12, "Times";
fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "";
fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times";
fontset = input, noPageBreakBelow, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, 12, "Courier";
fontset = output, output, inactive, noPageBreakBelow, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-4, 12, "Courier";
fontset = message, inactive, noPageBreakBelow, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier";
fontset = print, inactive, noPageBreakBelow, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier";
fontset = info, inactive, noPageBreakBelow, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier";
fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakBelow, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier";
fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, 10, "Times";
fontset = header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "";
fontset = leftheader, inactive, L2, 12, "Times";
fontset = footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, 12, "";
fontset = leftfooter, inactive, L2, 12, "Times";
fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times";
fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "";
fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Courier";
fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 14, "Times";
fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 14, "Times";
fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Courier";
fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "";
fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "";
paletteColors = 128; automaticGrouping; currentKernel;
]
:[font = title; inactive; preserveAspect; startGroup]
Example 09
Iterative Equations and
Modular Programming
:[font = subsubtitle; inactive; preserveAspect]
Copyright ã 1993 by Bill Titus, Carleton College,
Department of Physics and Astronomy, Northfield, MN 55057-4025
September 6, 1993
;[s]
3:0,0;10,1;11,0;133,-1;
2:2,16,12,Times,2,14,0,0,0;1,16,12,Symbol,0,14,0,0,0;
:[font = section; inactive; preserveAspect; startGroup]
Topics and Skills
:[font = special1; inactive; preserveAspect; endGroup]
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.
;[s]
15:0,0;9,1;12,0;42,1;49,0;56,1;66,0;80,1;90,0;127,1;138,0;160,1;169,0;201,1;209,0;241,-1;
2:8,16,12,Times,0,14,0,0,0;7,15,11,Courier,0,14,0,0,0;
:[font = section; inactive; preserveAspect; startGroup]
Problem
:[font = special1; inactive; preserveAspect]
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.
:[font = special1; inactive; preserveAspect]
Recall that in Euler's method, one breaks up the time domain into n equally spaced intervals of size dt:
;[s]
5:0,0;67,1;68,0;104,1;106,0;109,-1;
2:3,16,12,Times,0,14,0,0,0;2,15,11,Courier,0,14,0,0,0;
:[font = special1; inactive; preserveAspect; center]
to = t0 < t1 < t2 < ... < tn-1 < tn = tf
dt = (tf - to)/n
;[s]
19:0,0;1,1;2,0;6,1;7,0;11,1;12,0;16,1;17,0;29,1;33,0;36,1;37,0;41,1;43,0;50,1;52,0;55,1;56,0;60,-1;
2:10,15,11,Courier,0,14,0,0,0;9,23,15,Courier,64,14,0,0,0;
:[font = special1; inactive; preserveAspect]
For notational purposes, let
:[font = special1; inactive; preserveAspect]
x[ti] = xi, v[ti] = vi, a[x[ti],v[ti],ti] = ai.
;[s]
17:0,0;7,1;8,0;13,1;14,0;19,1;20,0;25,1;26,0;33,1;34,0;39,1;40,0;43,1;44,0;49,1;50,2;52,-1;
3:8,15,11,Courier,0,14,0,0,0;8,23,15,Courier,64,14,0,0,0;1,13,9,Times,0,12,0,0,0;
:[font = special1; inactive; preserveAspect]
Euler's iterative equations then become
:[font = special1; inactive; preserveAspect]
ti = ti-1 + dt
vi = vi-1 + ai-1 dt
xi = xi-1 + vi-1 dt
;[s]
17:0,0;8,1;9,2;10,1;16,2;19,1;30,2;32,1;37,2;40,1;44,2;47,1;56,2;57,1;63,2;66,1;70,2;77,-1;
3:1,16,12,Times,0,14,0,0,0;8,15,11,Courier,0,14,0,0,0;8,23,15,Courier,64,14,0,0,0;
:[font = special1; inactive; preserveAspect; endGroup]
These equations allow you to "bootstrap" up to xi and vi from the initial conditions: x0 = x[t0] and v0 = v[t0].
;[s]
19:0,0;48,1;49,2;50,0;57,1;58,2;59,0;92,1;93,2;95,1;100,2;101,1;103,0;110,1;111,2;113,1;118,2;119,1;120,0;122,-1;
3:5,16,12,Times,0,14,0,0,0;8,15,11,Courier,0,14,0,0,0;6,23,15,Courier,64,14,0,0,0;
:[font = section; inactive; preserveAspect; startGroup]
Solution - Using a Sequence of Iterative Equations
:[font = subsection; inactive; preserveAspect]
Step 1 - Briefly describe how you would implement Euler's algorithm if you were doing linear programming using arrays.
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 1
:[font = special1; inactive; preserveAspect; endGroup]
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.
;[s]
21:0,0;61,1;62,2;63,0;66,1;67,2;68,0;71,1;72,2;73,0;76,1;77,2;78,0;86,1;87,0;181,1;184,0;187,1;189,0;197,1;202,0;263,-1;
3:9,16,12,Times,0,14,0,0,0;8,15,11,Courier,0,14,0,0,0;4,23,15,Courier,64,14,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
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.
:[font = input; preserveAspect; endGroup]
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[" "];
];
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 2
:[font = special1; inactive; preserveAspect; endGroup]
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?
;[s]
3:0,0;145,1;147,0;361,-1;
2:2,16,12,Times,0,14,0,0,0;1,15,11,Courier,0,14,0,0,0;
:[font = subsection; inactive; preserveAspect]
Step 3 - Get help on For and summarize what you find.
;[s]
3:0,0;23,1;26,0;57,-1;
2:2,16,12,Times,1,14,0,0,0;1,15,11,Courier,1,14,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 3
:[font = special1; inactive; preserveAspect; endGroup]
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.
;[s]
3:0,0;37,1;40,0;173,-1;
2:2,16,12,Times,0,14,0,0,0;1,15,11,Courier,0,14,0,0,0;
:[font = subsection; inactive; preserveAspect]
Step 4 - Get help on Print[]. What do you find?
;[s]
3:0,0;22,1;29,0;50,-1;
2:2,16,12,Times,1,14,0,0,0;1,15,11,Courier,1,14,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 4
:[font = special1; inactive; preserveAspect; endGroup]
1. Print[] is useful for debugging a sequence of code.
;[s]
3:0,0;4,1;11,0;57,-1;
2:2,16,12,Times,0,14,0,0,0;1,15,11,Courier,0,14,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
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.
;[s]
5:0,0;141,1;143,0;153,1;154,0;237,-1;
2:3,16,12,Times,1,14,0,0,0;2,15,11,Courier,1,14,0,0,0;
:[font = special3; inactive; preserveAspect; endGroup]
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];
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 5
:[font = special1; inactive; preserveAspect; endGroup]
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]?
;[s]
3:0,0;150,1;162,0;164,-1;
2:2,16,12,Times,0,14,0,0,0;1,15,11,Courier,0,14,0,0,0;
:[font = subsection; inactive; preserveAspect]
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.
;[s]
7:0,0;14,1;21,0;45,1;51,0;89,1;97,0;130,-1;
2:4,16,12,Times,1,14,0,0,0;3,15,11,Courier,1,14,0,0,0;
:[font = subsection; inactive; preserveAspect]
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.
;[s]
3:0,0;39,1;45,0;144,-1;
2:2,16,12,Times,1,14,0,0,0;1,15,11,Courier,1,14,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 7
:[font = special1; inactive; preserveAspect; endGroup]
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?
;[s]
11:0,0;34,1;44,0;138,1;148,0;183,1;201,0;379,1;388,0;394,1;405,0;437,-1;
2:6,16,12,Times,0,14,0,0,0;5,15,11,Courier,0,14,0,0,0;
:[font = subsection; inactive; preserveAspect]
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.
;[s]
7:0,0;15,1;25,0;41,1;47,0;122,1;177,0;267,-1;
2:4,16,12,Times,1,14,0,0,0;3,15,11,Courier,1,14,0,0,0;
:[font = subsection; inactive; preserveAspect]
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.
;[s]
12:0,0;162,1;170,0;255,1;265,0;315,1;316,0;323,1;324,0;379,1;383,0;393,1;398,-1;
2:6,16,12,Times,1,14,0,0,0;6,15,11,Courier,1,14,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 9
:[font = special1; inactive; preserveAspect; endGroup; endGroup]
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.
;[s]
13:0,0;63,1;70,0;221,1;222,2;224,1;228,2;231,1;235,2;237,1;239,0;361,1;370,0;490,-1;
3:4,16,12,Times,0,14,0,0,0;6,15,11,Courier,0,14,0,0,0;3,23,15,Courier,64,14,0,0,0;
:[font = section; inactive; preserveAspect; startGroup]
Solution - Encapsulating the Euler Algorithm and Plotting Routine to a Single Function Using Module[]
;[s]
2:0,0;94,1;103,-1;
2:1,19,14,Times,1,18,0,0,0;1,18,13,Courier,1,18,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
Step 1 - Get help on the MMA function Module[].
;[s]
3:0,0;39,1;47,0;49,-1;
2:2,16,12,Times,1,14,0,0,0;1,15,11,Courier,1,14,0,0,0;
:[font = input; preserveAspect; startGroup]
?Module
:[font = info; inactive; preserveAspect; endGroup; endGroup]
Module[{x, y, ...}, expr] specifies that occurrences of
the symbols x, y, ... in expr should be treated as
local. Module[{x = x0, ...}, expr] defines initial
values for x, ....
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 1
:[font = special1; inactive; preserveAspect; endGroup]
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).
;[s]
7:0,0;4,1;12,0;132,1;144,0;381,1;393,0;494,-1;
2:4,16,12,Times,0,14,0,0,0;3,15,11,Courier,0,14,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
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.
:[font = input; preserveAspect]
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
];
:[font = input; preserveAspect; endGroup]
logisticPlot[0.4, 0.78, 20]
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 2
:[font = special1; inactive; preserveAspect; endGroup]
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[]?
;[s]
7:0,0;96,1;107,0;168,1;175,0;206,1;216,0;218,-1;
2:4,16,12,Times,0,14,0,0,0;3,15,11,Courier,0,14,0,0,0;
:[font = subsection; inactive; preserveAspect]
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.
;[s]
5:0,0;15,1;23,0;119,1;146,0;243,-1;
2:3,16,12,Times,1,14,0,0,0;2,15,11,Courier,1,14,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
Comment 3
:[font = special1; inactive; preserveAspect; endGroup]
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.
:[font = subsection; inactive; preserveAspect; startGroup]
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:
:[font = special2; inactive; preserveAspect; leftWrapOffset = 22]
On the Next:
1. Select Action \ Interrupt \ Interrupt Calculation, or -,
2. Select Action \ Interrupt \ Abort Calculation, or -.
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.
;[s]
22:0,0;27,1;69,0;74,1;84,2;85,0;99,1;137,0;143,1;154,0;168,1;209,0;252,1;269,0;314,1;326,0;386,3;397,0;410,1;414,0;438,1;440,0;573,-1;
4:11,16,12,Times,1,14,0,0,0;9,15,11,Helvetica,1,14,0,0,0;1,15,11,Courier,1,14,0,0,0;1,15,11,Helvetica,3,14,0,0,0;
:[font = special2; inactive; preserveAspect; leftWrapOffset = 22]
On the Mac:
1. Select Action \ Interrupt Calculation, or -,
2. Select Action \ Abort Calculation, or -.
3. Select Action \ Quit/Disconnect Kernel.
;[s]
11:0,0;26,1;58,0;64,1;75,0;89,1;119,0;122,1;133,0;148,1;181,0;183,-1;
2:6,16,12,Times,1,14,0,0,0;5,15,11,Helvetica,1,14,0,0,0;
:[font = special2; inactive; preserveAspect; leftWrapOffset = 22]
Now, activate the following cell and then abort.
:[font = input; preserveAspect; endGroup; endGroup]
For[i = 1, i < 1000000, i++]
:[font = section; inactive; preserveAspect; startGroup]
Save Your Notebook
:[font = special1; inactive; preserveAspect; endGroup; endGroup]
Remove any output cells from this notebook and then use Save As to store your notebook under the name myEx9.
;[s]
5:0,0;57,1;64,0;105,1;110,0;112,-1;
2:3,16,12,Times,0,14,0,0,0;2,15,11,Helvetica,0,14,0,0,0;
^*)