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.