1 UPJODE.HLP  Version 1.08  08 October 1987
  The University of Pittsburgh/Johnstown ODE Solver.
 
  The original version of this program was written by
  E Runnion, with modifications by L Foster and J Gentilesco
  at the University of Pittsburgh, Johnstown Campus.
 
  The program was rewritten by John Burkardt.
 
  UPJODE is designed as an interactive teaching aid.
  It allows a student with no programming skills to define a
  system of up to four first order differential equations,
  choose a solution procedure, and display the results in
  a table or plot.
 
2 The Problems UPJODE can solve
 
  The problems UPJODE can solve are known as systems of first
  order initial value problems.  They are characterized by the
  following features:
 
  There are N equations of the form
  D Y1/DT = F1(T, Y1, ... , YN)
  D Y2/DT = F2(T, Y1, ... , YN)
  through
  D YN/DT = FN(T, Y1, ... , YN)
  For example, D Y3/DT = 1 + T + Y2*Y4 + SIN(Y2)
  Sometimes D Y3/DT is written Y3' or called "Y 3 Prime" or
  "Y 3 Dot".
 
  There is an initial value of time given, at which the solution
  is known.  In the program, we may sometimes call these
  values TINIT, Y1-INIT, Y2-INIT, and so on.
 
  It is desired to know the values of Y as T increases (or
  decreases) from its current value to some new, specified value.
 
  Interestingly, problems involving higher order derivatives
  can sometimes be rewritten and solved as first order problems.
  For example,
 
  Y" + Y = 0
  with initial values T=0, Y(0)=0, Y'(0)=0.
 
  If we simply make up new names for Y and Y', calling them
  Y1 and Y2, then right off the bat we know that
 
  Y1'=Y2
 
  and since Y" can be written as Y2', our original equation
  can be written as
 
  Y2' + Y1 = 0, or Y2' = -Y1.
 
  and so our original problem is the same as the problem
 
  Y1' = Y2
  Y2' = -Y1
  with initial values T=0, Y1(0)=0, Y2(0)=0.
 
2 Summary of Available ODE Methods
 
  We are asked to integrate from some original value T
  to some stopping value TSTOP.  All the methods we will
  describe do this by taking a certain number of steps
  NSTEPS from T to TSTOP.  The way UPJODE works, the
  number of steps is decided beforehand, and the steps are
  all the same size, TDELT=(TSTOP-T)/NSTEPS.  (Some programs
  decide the size of their steps as they go along.  They even
  change the method they will use.  They are known as adaptive
  methods, and are much better, but harder to understand.)
 
  Let's ignore the possibility that we are integrating more
  than one equation.  Then UPJODE will be producing a string
  of times and values: T(0), Y(0), then T(1), Y(1), then T(2), Y(2),
  and so on through T(NSTEPS), Y(NSTEPS).  At any point in the process,
  say at step K, we have all the values from T(0), Y(0) through
  T(K), Y(K) and we want to produce T(K+1), Y(K+1).  Of course,
  T(K+1) is easy, it's just T(K) + TDELT.  But the hard part is
  coming up with a value Y(K+1).  The methods we use can all
  be described by how they do this.
 
  From this point on, we will change notation somewhat.
  Instead of TDELT, we will write H.  We also assume the right hand
  side of our problem is F(T,Y).
 
  Many of the methods produce a first guess of the solution
  or sometimes two or three, before finally getting the
  estimate of Y(K+1).  These tentative values will show up
  labeled with a Z.
 
  Also, most of the methods can make an estimate of how far
  they are off.  Actually, they can only estimate how much
  more error they have added to what has already been made.
  This is known as the local error, and it is listed with
  the method.  Sometimes the actual formula is complicated,
  and so the 'O' notation is used.  The important part of
  the error is the term involving the smallest power of H.
  O(H**2) means that the smallest power of H is 2 in the
  error estimate.  An error estimate of 98*H**3 + 1.7*H**5
  would be described as O(H**3).
 
  Truncation error measures how well the true solution satisfies
  the difference equation which replaces the differential equation.
  Our approximate methods can be viewed as replacing the equation
  y'=f(t,y) with some difference method.  For Euler's method,
  we would write (y(k+1)-y(k))/h.  Then truncation error is
  defined as the value of y'(t+h)-(y(k+1)-y(k))/h which equals
  f(t+h,y(k+1))-(y(k+1)-y(k))/h which we can estimate using
  Taylor's formula, which gives us to 0.5*h*y''(xsi)
  for some xsi between t(k) and t(k+1), so the truncation error
  is of order h.  This also is measuring the local error in
  y', so we have locally an error of order h in y', and an
  error of order h**2 in y.  This is typical for ODE methods,
  and usually results in a global error for y and y' that is
  one order lower.  Thus for the Euler method, we can assume
  that globally, our answers have an error term of order h in
  y, and of order 1 (doesn't go to zero with h) in y'.
 
  Finally, note that some methods require 2 or more old
  values to make a new value.  Such methods cannot start up
  from scratch, and so, at the initial value Y(0) they have
  to use a method like a Runge Kutta method which can get
  going with only one previous value.  See for instance
  the Adams-Bashforth methods.
 
3 Euler's Method
  Taylor Method of order 1
  Runge Kutta Method of order 1
 
  Y(K+1) = Y(K) + H * F(T(K), Y(K)).
  Local error = 0.5*H**2*Y''()
 
3 Midpoint Method
  (Half-step Euler Method)
 
  Z(K+1) = Y(K) + 0.5 * H * F(T(K), Y(K))
  Y(K+1) = Y(K) + H * F(T(K)+0.5 * H, Z(K+1))
 
3 Runge-Kutta Method of Order 2
  Heun's Method
 
  Z(K+1) = Y(K) + H * F(T(K), Y(K))
  Y(K+1) = Y(K) + H * (F(T(K+1), Z(K+1)) + F(T(K), Y(K))) / 2.0
  Local error = O(H**3)
 
3 Runge Kutta Method of Order 4
  (Using Runge's constants)
 
  Z1     = F(T(K), Y(K))
  Z2     = F(T(K)+0.5*H, Y(K)+0.5*H*Z1)
  Z3     = F(T(K)+0.5*H, Y(K)+0.5*H*Z2)
  Z4     = F(T(K)+    H, Y(K)+    H*Z3)
  Y(K+1) = Y(K) + H*(Z1 + 2*Z2 + 2*Z3 + Z4) / 6.0
  Local error = O(H**5)
 
3 Runge Kutta Fehlberg method of orders 5 and 6
 
  Z1     = F(T(K), Y(K))
  Z2     = F(T(K) + 0.25 * H, Y(K) + 0.25 * H * Z1)
  Z3     = F(T(K) + 3 * H / 8, Y(K) + 3 * H * Z1/32
                                    + 9 * H * Z2/32)
  Z4     = F(T(K) + 12 * H / 13, Y(K) + 1932 * H * Z1 / 2197
                                      - 7200 * H * Z2 / 2197
                                      + 7296 * H * Z3 / 2197)
  Z5     = F(T(K) + H, Y(K) +  439 * H * Z1 / 216
                            -    8 * H * Z2
                            + 3680 * H * Z3 / 513
                            -  845 * H * Z4 / 4104)
  Z6     = F(T(K) + H / 2, Y(K) -    8 * H * Z1 / 27
                                +    2 * H * Z2
                                - 3544 * H * Z3 / 2565
                                + 1859 * H * Z4 / 4104
                                -   11 * H * Z5 / 40)
 
  Y(K+1) = Y(K) + H *    25 * Z1 / 216
                + H *  1408 * Z3 / 2565
                + H *  2197 * Z4 / 4104
                - H *         Z5 / 5)
  W(K+1) = Y(K) + H *    16 * Z1 / 135
                + H *  6656 * Z3 / 12825
                + H * 28561 * Z4 / 56430
                - H *     9 * Z5 / 50
                + H *     2 * Z6 / 55 )
 
  Local error = W(K+1) - Y(K+1)
 
3 Adams-Bashforth method of order 2
 
  Y(0) is the initial condition.
  Y(1) is generated by Runge-Kutta method of order 2.
 
  Y(K+1) = Y(K) + 0.5 * H * (3 * F(T(K),Y(K)) - F(T(K-1),Y(K-1)))
  Local error = 5*Y'''()*H**3/12
 
3 Adams-Bashforth method of order 4
 
  Y(0) is the initial condition.
  Y(1), Y(2) and Y(3) are generated by the Runge Kutta method
  of order 4.
 
  Y(K+1) = Y(K) + H * 55 * F(T(K),Y(K))     / 24
                - H * 59 * F(T(K-1),Y(K-1)) / 24
                + H * 37 * F(T(K-2),Y(K-2)) / 24
                - H *  9 * F(T(K-3),Y(K-3)) / 24
  Local error = 251*Y'''''()*H**5/720
 
3 Adams-Bashforth-Moulton predictor/corrector of order 2
 
  Y(0) is the initial condition.
  Y(1) is generated by Runge Kutta method of order 2.
 
  Z(K+1) = Y(K) + 0.5 * H * (3 * F(T(K),Y(K)) - F(T(K-1),Y(K-1)))
  Y(K+1) = Y(K) + 0.5 * H * (F(T(K+1),Z(K+1)) + F(T(K),Y(K)))
  Local error = (Z(K+1)-Y(K+1))/6.0
 
3 Adams-Bashforth-Moulton Predictor/Corrector of order 4
 
  Y(0) is the initial condition.
  Y(1) through Y(3) are gotten by Runge-Kutta method of order 4.
 
  Z(K+1) = Y(K) + H * 55 * F(T(K),Y(K))     / 24
                - H * 59 * F(T(K-1),Y(K-1)) / 24
                + H * 37 * F(T(K-2),Y(K-2)) / 24
                - H *  9 * F(T(K-3),Y(K-3)) / 24
  Y(K+1) = Y(K) + H *  9 * F(T(K+1),Z(K+1)) / 24.0
                + H * 19 * F(T(K),Y(K))     / 24.0
                - H *  5 * F(T(K-1),Y(K-1)) / 24.0
                + H *      F(T(K-2),Y(K-2)) / 24.0
  Local error = 19 * (Z(K+1)-Y(K+1))/ 270.0
 
3 Milne/Simpson Predictor/Corrector of order 3
 
  Y(0) is the initial condition.
  Y(1), Y(2) and Y(3) are generated by Runge Kutta method of order 4.
 
  This method is "weakly unstable".
 
  Z(K+1) = Y(K-3) + H * 8 * F(T(K),Y(K))     / 3.0
                  - H * 4 * F(T(K-1),Y(K-1)) / 3.0
                  + H * 8 * F(T(K-2),Y(K-2)) / 3.0
 
  Y(K+1) = Y(K-1) + H *     F(T(K+1),Z(K+1)) / 3.0
                  + H * 4 * F(T(K),Y(K))     / 3.0
                  + H *     F(T(K-1),Y(K-1)) / 3.0
  Local error = (Z(K+1)-Y(K+1))/29.0
 
2 UPJODE commands
 
  All action taken by UPJODE must be specified by you through
  commands from a menu.  These commands tell UPJODE you want
  to set up a problem, specify a solution method, integrate
  from the starting point, continue integrating, output a file,
  and so on.

  Notice in particular the following:
 
  If you have solved a problem, and want to solve the same problem
  with just a few parameters changed, you must reset the value of T.
  For example, to resolve a problem with METHOD=3, which started at
  T=0, use the commands:

  S T=0.0
  S METHOD=3
  C
 
  On the other hand, you might want to simply continue the integration.
  You've solved from 0.0 to 3.0, and now you want to continue to 6.0.
  In that case, you must reset TSTOP before trying to continue.
  Use the commands:
 
  S TSTOP=6.0
  C

3 B - Set up new problem
 
  Before you can do anything with UPJODE, a certain amount of
  information must be entered to define a problem.  Thereafter,
  new, interesting problems can be made by changing just one
  of the items.  Usually, all of this information is input through
  the S command.  But that's hard to get used to at first.  If
  you want UPJODE to lead you through all the steps required to
  set up a new problem, this is the command to use.
 
3 C - begin or continue integration
 
  This command is used once you have set up the problem by
  defining the number of equations, NEQN, the intial conditions,
  T, Y1, perhaps Y2, Y3, Y4, the number of steps NSTEPS,
  the stopping point TSTOP, the method METHOD,
  and the right hand sides defining Y1', and perhaps Y2', Y3', Y4'.
  The program will carry out the integration from the current
  time T to TSTOP, printing out every IPRINT-th step.
 
3 D - (filename) open/close transcript file
 
  If you want a permanent disk file record of your work,
  using this command will open a file in which an exact
  copy of what appears on the screen will be kept.  A second
  D will close the file. If you don't specify a name, the
  default name UPJODE.LPT will be used.
 
3 F - enter right hand sides and perhaps exact solution
 
  Before issuing this command, you must have set the number
  of equations, NEQN, with a command like 'S NEQN=1'.
  This command is required before you can integrate with the C
  command.  For each equation, you must type in the right hand side,
  using a format for the formulas that is explained elsewhere.  For
  example, a system Y1'=T, Y2'=Y1+2*SIN(Y2) would be entered
  as T on one line, then Y1+2*SIN(Y2) on the next line.  For each
  equation, you can also give the exact solution, or hit the
  return key.
 
3 G - (x variable, y variable) produce graph
 
  This command produces a graph of the data, where you
  can plot one of T Y1 Y1' Y2 Y2' etc against any other variable.
  You can plot several quantities on the same graph.  When you
  have done with your list of pairs of things to plot, just
  hit the return key so the plot will be displayed.  Hit return
  again when you are done looking at the plot.

  For an IBM you should have CGA graphics available.  
  For a Macintosh, graphics is always available.
  For a mainframe machine, you should probably have a Tektronix Plot-10
  capable terminal, such as a VT240.

  If your terminal or PC does not have the proper graphics capability,
  then the G command will stand for 'Garbage'.
 
3 H - Print the help list
 
  This command prints a list of one line summaries of commands.
 
3 I - set initial conditions
 
  This command should be issued before using the B command.
  You must give the initial time, and for each equation,
  the initial value of Y.
 
3 J - set up sample problem.
 
  The sample problem Y1'=COS(T) is set up and solved.
 
3 M - (method or 0) choose method or get list
 
  This command allows you to select the integration method to
  use, or to see a list of what is available.
 
3 N - Make a note
 
  Particularly if you are using a transcript file, you might
  want to use this command.  Having typed an N you may then
  type several lines of comments, which will appear in your
  transcript file too.  Finish up by typing a period in column 1.

3 O - (filename) read saved data from file

  This command can be used to read data from a file created with
  the P command.  While the P command is intended to dump the
  UPJODE data to a file for plotting by another program, this is
  also a handy way to save data and have UPJODE look at it later.
  You can plot the data with UPJODE's internal plotter but, because
  the right hand sides are not saved, this data is NOT sufficient
  to allow you to continue an integration.  If you try to do so,
  you will get a weird mish-mash calculation involving the current
  data plus old data that was lying around in the program.

3 P - (filename) save data to file

  This command will write the values of all the variables in a
  file which might then be passed to a plotting program.  Note
  that there are a few header lines, then the word 'TABLE',
  then the values of T, Y1, Y1', and so on.  This option will
  not work well if there are more than 5 quantities to be saved,
  since only 5 numbers per line are written.
 
3 Q - Quit
 
  The program will ask you nicely if you really want to quit.
  Type Y if you meant it.
 
3 S - (parameter=value) Assign the value of a parameter
 
  This command allows you to assign values to the parameters.
  The parameters you can set include

  IPRINT - the frequency of output during integration.
  METHOD - the integration method.
  NEQN   - the number of equations.
  NSTEPS - the number of steps to take from the current time to TSTOP.
  NTAB   - the total number of steps taken.
  ROVER  - the overflow threshhold.  Computation is stopped if
           any component of Y or Y' exceeds ROVER in absolute value.
  T      - the initial time (automatically restarts problem).
           You must reset T to the original value in order to
           restart a problem at the beginning.
  TSTOP  - the stopping time.  You must reset TSTOP if you want
           to continue an integration from the point where you've
           just stopped.
  Y1     - the initial value of Y1 (automatic restart)
  Y2     - the initial value of Y2, and so on.

  Sample commands include 'S NEQN=1', 'S ROVER=100000.0'.

  Note that to reset the formula for Y1', you must use the F command
  instead of the S command.
 
3 T - type out parameters
 
  In case you forgot what you have set up, you can use this
  command to see the values of the problem parameters.
 
3 V - what parameters can I set?
 
  If you want to be reminded of the parameters you can set with the
  S command, this command will list them for you.
 
3 W - (filename) write output to file
 
  This command writes a description of the current problem into
  a file which can be read back in later.
 
3 X - (filename) read input from file
 
  You can use this command to have UPJODE read some canned
  input from a file.  In this way, problems with difficult
  to type right hand sides can be prepared beforehand.
  A sample input file named UPJODE.DAT is available.
  If you type USER as the name of the file, the program expects
  to read input from the terminal.  This can be the last line
  of your input file, switching control back to you.
 
3 ? - Display extensive help from UPJODE.HLP
 
  This command can be used to read this information while
  running the program.
 
2 Formulas
 
  To set up the right hand sides correctly, you must use the
  proper abbreviations for mathematical functions, and the right
  symbols for operations like multiplication.
 
3 SYMBOLS AND CONSTANTS
 
  * is multiplication, / is division, + is addition,
  - subtraction, and ** exponentiation.
 
3 FUNCTIONS
 
 
  ABS(*)            - Absolute value of *.
  ACOS(S)           - The arc cosine of S.          -1 < S < 1
  ALOG(S)           - Natural logarithm of S.            S > 0
  ALOG10(S)         - Logarithm base 10 of S.            S > 0
  ASIN(S)           - Arc sine of S.                -1 < S < 1
  ATAN(S)           - Arc tangent of S.
  ATAN2(S1,S2)      - Arc tangent of (S1/S2).  Correctly computes
                      ATAN2(0,0)=PI/2.
  COS(S)            - Cosine of S.
  COSH(S)           - Hyperbolic cosine of S.
  EPS               - Machine epsilon.  This is the smallest power of two
                      such that 1+EPS.GT.1.0
  EXP(S)            - Exponential of S.
  LN(S)             - Natural logarithm of S.   0 < S
  LOG(S)            - Natural logarithm of S.   0 < S
  LOG10(S)          - Logarithm base 10 of S.   0 < S
  MAX(S1,S2)        - Maximum of S1, S2.
  MIN(S1,S2)        - Minimum of S1, S2.
  NEG(S)            - Changes sign of S.
  PI                - 3.14159265...
  RAN(S)            - Computes a random number between 0 and 1.
  SIN(S)            - Sine of S.
  SINE(S)           - Sine of S.
  SINH(S)           - Hyperbolic sine of S.
  SQRT(S)           - Square root of S.           0 < S
  STEP(S)           - Heavyside step function.
                      STEP=0 if S<0, STEP=1 if S>0.
  TAN(S)            - Tangent of S.
  TANH(S)           - Hyperbolic tangent of S.
  I!                - I factorial.  Valid for I=0 to 25.
                      I!=I*(I-1)*(I-2)*...*2*1.
 
2 Sample problems
 
  The following problems can be interesting to run and plot:
 
3 Predator Prey
 
  If the guppy population increases with rate A, and if the piranha
  population would die out at rate C without any guppies, and if
  the chances of a guppy being eaten by a piranha are proportional
  to the product of the number of guppies and piranha, which reduces
  the number of guppies and increases the number of piranha, then
  the equations would be
 
  Y1' =   A * Y1 - B * Y1 * Y2
  Y2' = - C * Y2 + D * Y1 * Y2
 
  Equilibrium solution is Y1=C/D, Y2=A/B.
 
3 Predator Prey with fishing
 
  In the above model, include the effects of fishing by outside forces
  on both populations.
 
  Y1' =   A * Y1 - B * Y1 * Y2 - F * Y1
  Y2' = - C * Y2 + D * Y1 * Y2 - F * Y2
 
  Equilibrium solution is Y1=(C+F)/D, Y2=(A-D)/D
 
3 Linear pendulum
 
  THETA'' = -(G/L)*SIN(THETA)
 
  with THETA the angle from the vertical, G the gravitational acceleration,
  and L the length.  Replacing SIN(THETA) by THETA, assuming that THETA
  is small, and rewriting as two first order equations:
 
  Y1' = Y2
  Y2' = -(G/L)*Y1    or Y2' = -(G/L)*SIN(Y1) for large ranges of Y1.
 
3 Duffing's Equation

  Describes a nonlinear spring with linear damping and a periodic
  forcing term.
 
  X'' + K*X + L*X*X*X = A*COS(OMEGA*T)

  or, in the form we can use, with Y1=X, Y2=X',
 
  Y1' = Y2
  Y2' = - K*Y1 - L*Y1*Y1*Y1 - A*COS(OMEGA*T)
 
3 Van Der Pol's Equation
 
  X'' + EPS*(X*X-1)*X'+ X = 0
  X(0)=X'(0)=1
  or, in the form we can use, with Y1=X, Y2=X',
 
  Y1' = Y2
  Y2' = -EPS*(Y1*Y1-1)*Y2 - Y1
  Y1(0)=Y2(0)=1

  When EPS is equal to 1 or more, the graph of Y1(T) is interesting,
  as well as a phase plane plot of Y1 versus Y1' (=Y2).  Integrate
  from T=0 to T=15.

3 Blow up

  Y1' = T*T + Y1*Y1
  Y1(0)=1
 
  Integrate this equation from T=0 to T=1.  The exact solution
  blows up at 1.  Approximate methods will compute a finite value
  at T=1, but as the stepsize decreases, the solution computed at
  1 will increase.  For a Runge-Kutta method of order 2 or 4, you
  should get overflow if you take about 30 or 40 steps from 0 to 1.
 
2 Compilation and running
 
  VMS VAX:   Type 'RUN MTH:UPJODE'.  Note that using graphics will
  require that your terminal is capable of emulating a Tektronix
  terminal.  The VT240 terminal is acceptable, although you must
  put the terminal in Tektronix mode before doing graphics.  Also,
  if you are using a PC as a terminal, some terminal programs like
  SmarTerm 240, can be put in Tektronix mode.  Otherwise, attempts
  to draw a graph will deposit gibberish on your screen.
 
  IBM PC:    If your disk is NOT self booting, then first
  insert into drive A a disk with the system file COMMAND.COM
  on it.  Turn on the machine, allowing it to go through its
  system checks, until you see the prompt sign 'A>'.  At that
  point, remove the original disk from drive A and insert your
  UPJODE disk.  Now type 'UPJODE' and the program should run.
  When you are done, just turn the machine off.
 
  To make a self-booting copy:  I assume you have access to
  a dual floppy system.  Boot up the machine.  Using your
  DOS disk (the one that came with the machine, which has
  the program FORMAT.COM on it, in particular) in drive A,
  type the command 'FORMAT B:/S'.  The program will tell you
  to insert a blank disk in drive B and hit return.  Do so.
  The disk will be formatted, and eventually the system will
  ask you if you want to format another.  Say 'N'.  Now
  remove your DOS disk from drive A, and replace it with
  your UPJODE disk.  Type 'COPY A:*.* B:'  This should copy
  all the UPJODE files onto the disk in drive B, which is
  a self booting disk.
 
  If you've made a self-booting UPJODE disk, then in order
  to run the program, you just stick the disk in drive A
  and turn on the machine, and type UPJODE.
 
  MACINTOSH: Insert the disk containing UPJODE,
             double click on the UPJODE disk to open it,
             double click on the file 'UPJODE APL'.
 
  On a Macintosh, because there are large arrays declared in
  a subroutine (COMRPN), special attention must be paid at link
  time.  Before linking, the command 'Z 50' to increase the heap
  size (whatever that is) to 50K must be used. 
 
  Also, the use of Macintosh graphics requires that the link list
  include toolbx.sub and toolbx.par be on the link disk.
