home *** CD-ROM | disk | FTP | other *** search
- ╔══════════════════════════════╗ (C) Copyright 1986-1990
- ║ WELCOME to the ODE Solver ! ║ Zvi Shippony
- ╚══════════════════════════════╝ (818) 990-0134
-
- This Program solves Ordinary Differential Equations. The main use is
- solving Initial Condition (I.C) type problems but for second order equation
- the user have the option to solve a Two-Point-Boundary Value problem (B.V).
-
- Option #1:
- ----------
- Solves a set of N First order ODE whose Right-Hand-Side (RHS) is any
- legal expression of the form: F(X,Y0,Y1,Y2,...,Y9) (Thus N <= 10 ).
-
- Option #2:
- ----------
- Solves an N-th order ODE whose RHS is any legal expression of the form:
- F(X,Y,Y',Y'',Y''',..,Y''''''''') or F(X,Y0,Y1,Y2,Y3,..,Y9) where Y0 stands
- for Y, Y1 stands for Y', Y2 stands for Y'' and so on. YOU CAN USE BOTH
- NOTATIONS.
- $$$
-
- 'Expression' is any legal combination of: +, -, *, /, **, (, ), AND any of
- the following functions:
-
- ABS, INT, EXP, SIN, COS, TAN, COT, LOG, LN, FACT or ! (Factorial)
- SQRT, SINH, COSH, TANH, ARCSIN, ARCCOS, ARCTAN, ARCSINH, ARCCOSH, ARCTANH
-
- And the "Special Functions" :
-
- Z(x) { Riemann's Zeta function }
- G(x) { Gamma function, (IF x is an integer then x! = G(x+1)) }
- BJ(n,x) { Bessel Function of the first kind, J(n,x) }
- BY(n,x) { Bessel Function of the second kind, Y(n,x) }
- BI(n,x) { Modified Bessel Function of the first kind, I(n,x) }
- BK(n,x) { Modified Bessel Function of the second kind, K(n,x) }
- SBJ(n,x) { Spherical Bessel Function of the first kind, j(n,x) }
- SBY(n,x) { Spherical Bessel Function of the second kind, y(n,x) }
-
- PI is a reserved name and will be interpeted as Pi = 3.14159265358....
-
- $$$
- Interval of solution (A,B) is needed, where A and/or B are either given
- as numbers or any expression containing numbers. A comma is required between
- A and B .
-
- For example: if A = 0.0 and B = Pi/4.0, the input could be given as:
- 0.0 , pi/4 or: 0.0 , arctan(1.0) .
-
- Initial Conditions (or Boundary Values) could be given as numbers or any
- legal expression containing numbers .
-
- NOTE ON TWO-POINT-BOUNDARY-VALUE PROBLEMS: (Applies for 2-nd order eqn. only)
- -----------------------------------------------------------------------------
- Only explicit boundary values are accepted. That is, boundary value can
- not be a function of Y itself.
- ( Right example: Y(1.0)=5*exp(1) Wrong example: Y(1.0)=5*Y(0.0) )
-
- The method for solving boundary value problem is a combination of FINITE-
- DIFFERENCE Method and the SHOOTING Method. The user can supply the program
- with the partial derivatives of F(X,Y,Y') w.r.t. Y and Y', this way the
- convergence will be faster (using Newton's Method) but if these are not
- supplied, the iterations will be done with the False-Position Method, which
- gives slower convergence (and the FINITE-DIFFERENCE step will NOT be used).
- $$$
- Accuracy of solution:
- ---------------------
- The program will try to achieve the inputted accuracy, but the higher
- the requested accurate is, the longer the execution time will be. Some time
- roundoff errors will stop the execution without achieving the desired
- accuracy. The smallest possible accuracy is: 5.0E-14 (default).
-
- The method used for all I.C. problems is the Bulirsch_Stoer integration
- procedure with 4-th order Runge-Kutta step taken at points of non-convergence.
-
- When the program writes: '** Finished' , the solution is stored in an
- array of up to 100 entries equally spaced.
-
- The user has the option to ask for any specific X-value (as long as this
- value is in the interval [A,B] ) and the program will interpolate for that
- value.
-
- Another way is to store these entries in a file for later use (Plotting,
- Printing, etc.)
- $$$
-
- Option #3: Loading a system from a file
- ----------
- Any ASCII file written in the format bellow could be used to load a pre-
- determined system into the program for solving. The file MUST have that
- format in order to work properly. It could be written in any editing type
- S/W (like: Edlin, WordStar, Word-Perfect, Word etc..). Each line has a colon
- (':') separating between the name of the input and its value.
-
- Note: You can insert comments anywhere in the line as { ... }.
-
-
- There are 2 slightly different file formats, according to the option needed:
- $$$
-
- File format for option 1 (System of N first order equations):
- -------------------------------------------------------------
-
- Option : 1
- Order : {1 to 10}
- Type : IC
- Domain : {X-Range in the format: x1 , x2} (Don't forget the comma !!)
- Accuracy: {Integration accuracy to use}
- StepSize: {Step size to use}
- Eqn. #0 : {The Right-hand-size expression to use for Y0'}
- Eqn. #1 : {The Right-hand-size expression to use for Y1'}
- Eqn. #2 : {The Right-hand-size expression to use for Y2'}
- . .
- . .
- Eqn. #n : {Where n is = order -1 }
- Initial condition Y0(LHS): {Y value at the initial point}
- Initial condition Y1(LHS): {Y' value at the initial point}
- . .
- . .
- Initial condition Yn(LHS): {Where n is = order-1}
-
- $$$
- File format for option 2 (N-th order equation):
- -----------------------------------------------
-
- Option : 2
- Order : {1 to 10}
- Type : {IC or BV}
- Domain : {X-Range in the format: x1 , x2} (Don't forget the comma !!)
- Accuracy: {Integration accuracy to use}
- StepSize: {Step size to use}
- Equation: {The Right-hand-size expression to use}
- dF/dY : {FOR BV PROBLEMS ONLY, the derivative of F(X,Y,Y') w.r.t Y }
- dF/dY' : {FOR BV PROBLEMS ONLY, the derivative of F(X,Y,Y') w.r.t Y'}
- Boundary Value Y0(LHS) : {The Left-Hand_Side Y value }
- Boundary Value Y0(RHS) : {The Right-Hand-Side Y value}
- or:
- Initial condition Y0(LHS): {Y value at the initial point}
- Initial condition Y1(LHS): {Y' value at the initial point}
- . .
- . .
- Initial condition Yn(LHS): {Where n is the order-1}
- $$$
-
- NOTE:
-
- Some times the equation (the 7-th input) is too long to fit on one line.
- In this case you can continue it on the next line by starting it with
- a colon (':') and then type in the rest of the equation, i.e. :
-
- Equation: (1 + sin(y') - cos(y''*y''') + tan(y''''*y''''') -
- : y''*y'''') / sqrt(y''**2 + y'''**2)
-
-
- YOU CAN USE AS MANY CONTINUATIONS YOU WANT,BUT AFTER SQUEEZING ALL THE
- EXPRESSION(s) FROM BLANKS AND CONVERTING ALL THE PRIME NOTATIONS
- (Y',Y'',Y''',...) INTO (Y1,Y2,Y3,..) THE FINAL EQUATION CAN NOT BE LONGER
- THEN 76 CHARACTERS !!
-
- $$$
-
- Option #4: Storing a system into a file
- ----------
- After a system has been solved, you can store it into a file. The program
- will write an ASCII file in the format given in the above example. This file
- can be then edited and used later as loading file (as in option 3 above)
-
- $$$
-
-
-
-
- WE ARE GOING TO WALK THROUGH 2 COMPLETE EXAMPLES NOW ..
-
- $$$
-
- Example #1: After the program is loaded, it starts like:
- -----------
-
- Your options are:
-
- 1 - Solve system of N first order O.D.E.: Yi' = F(X,Yj) (i,j = 0,1,..,N-1)
- 2 - Solve N-th order O.D.E., e.g: Y''' = F(X,Y0,Y1,Y2)
- Where: Y0 ≡ Y, Y1 ≡ Y', Y2 ≡ Y'' (You can use either notation..)
- 3 - Load a system from a file
- 4 - Store last system on file
- M - Memory manipulations
- B - Set Bell ON / OFF for error or info. msg. Bell is now: OFF
- H - Help and information about the program
- - Any other key to quit ..
-
- Enter your Choice (1 - 4, M, B, H or Quit):
-
- ( LET US CHOOSE 2nd-ORDER 2-POINTS BOUNDARY-VALUE PROBLEM: Y''= -(2*Y'+Y)
- Solution: Y(x)=X*EXP(-X) , FOR X=[-1,4], Y(-1)=-EXP(1), Y(4)=4*EXP(-4) )
-
- Enter your Choice (1, 2, 3, B, H or Quit): 2
- Store various inputs in Memory ? (y/n): n (SAY: N)
- Enter order of ODE (1-10): 2
- $$$
- Enter: IC for Initial Condition Problem
- Enter: BV for Two-Point-Boundary Value Problem
- Your choice: BV
-
- (THE PROGRAM WILL NOW PRESENTS YOU WITH THE LINE: Y'' =
- THIS IS THE TIME TO ENTER THE EQUATION(S) ...)
-
- Y'' = -(2*y'+y) ( SO F(X,Y,Y') = -2*Y'-Y )
- (YOU COULD ALSO ENTER: -(2*Y1+Y0) or: -(2*Y'+Y0) etc ..)
-
- ** For faster convergence we need dF/dY and dF/dY'
- Are you supplying those ? (y/n): y (SAY: Y)
-
- dF/dY = -1 (THIS IS DF/DY )
- dF/dY' = -2 (THIS IS DF/DY')
-
- Enter x Domain ( e.g.: -1.0 , 1.0): -1,4
- Enter desired accuracy (Default is: 5.0E-14): 1.e-10
- Compute solution every what Delta-x ? 0.05 (Total of 100 points)
-
- Enter B.V: Y(-1.0) (Default: 0.0): -exp(1)
- Enter B.V: Y(4.0) (Default: 0.0): 4*exp(-4)
- ** Do you wish to stop (pause) after each iteration ? (y/n): n (SAY: N )
- ( YOU HAVE THE OPTION TO CONTROL THE ITERATIONS PROCESS )
- $$$
- ( SINCE WE CHOSE 2-POINT BOUNDARY-VALUE PROBLEM, THE FIRST METHOD IS THE
- FINIT-DIFFERENCE METHOD ..THE PROGRAM WILL SAY:)
-
- ** Iterative Finite-Difference Method **
- # of iterations: 2 Yields: dy/dx(LHS) = 5.43638706110044E+000
- Number of mesh points: 100 Step size: 4.95050E-002
- Solution accuracy at mesh points: Order(2.0E-3)
-
- Are you satisfied with the current accuracy ? (y/n) n (SAY: N )
-
- Start Shooting Method now with the above guess for dy/dx(LHS)
-
- ( Working ... )
- $$$
-
- ** Iterative Shooting Method **
- Given B.V: 7.32625555549367E-002 Achieved B.V: 7.32625555469740E-002
- Rel. Error: 7.96270E-012 Using dy/dx(LHS) = 5.43656365665446E+000
- Integration accuracy: 1.00000E-010 Max. step size used: 5.00000E-002
- Next guess for dy/dx(LHS) is: 5.43656365668701E+000
-
- ( Working ... )
-
- *** Finished ... (HOORAY ! IT CONVERGED ...)
-
- To Plot or do Least-Square-Fit of the solution we need to write the results
- into a file. Write a file ? (y/n): (SAY: Y)
- Enter file name: ODETST
-
- First title for the file is: SOLUTION OF: Y'' = -(2*Y1+Y0)
- Enter (optional) 2-nd title (Initial or Boundary conditions, remarks, etc ..)
- Title:Testing. BV: y(-1)=-exp(1), y(4)=4*exp(-4) eps=1.0e-10
- (THE PROGRAM WILL RESPOND BY:) ** File: ODETST is ready .
- $$$
- Example #2:
- -----------
- Your options are:
-
- 1 - Solve system of N first order O.D.E.: Yi' = F(X,Yj) (i,j = 0,1,..,N-1)
- 2 - Solve N-th order O.D.E., e.g: Y''' = F(X,Y0,Y1,Y2)
- Where: Y0 ≡ Y, Y1 ≡ Y', Y2 ≡ Y'' (You can use either notation..)
- 3 - Load a system from a file
- 4 - Store last system on file
- M - Memory manipulations
- B - Set Bell ON / OFF for error or info. msg. Bell is now: OFF
- H - Help and information about the program
- - Any other key to quit ..
-
- Enter your Choice (1 - 4, M, B, H or Quit):
- ( LET US CHOOSE SYSTEM OF TWO 1-st ORDER EQUATIONS)
-
- Enter your Choice (1, 2, 3, B, H or Quit): 1
- Store various inputs in Memory ? (y/n): n (SAY: N)
- Enter number of equations (1-10): 2
-
- Y0' = Y1
- Y1' = -Y0 ( The solution is: Y0(X)=SIN(X), Y1(X)=COS(X) )
- $$$
-
- Enter x Domain ( e.g.: -1.0 , 1.0): -PI,PI
- Enter desired accuracy (Default is: 5.0E-14): 1.E-8
- Compute solution every what Delta-x ? 0.1
-
- Enter I.C: Y0(-3.14159265358979) (Default: 0.0): 0.0
- Enter I.C: Y1(-3.14159265358979) (Default: 0.0): -1.0
-
- ( Working ... )
-
- *** Finished ...
-
- To Plot or do Least-Square-Fit of the solution we need to write the results
- into a file. Write a file ? (y/n): (SAY: N)
- ( WE ARE SKIPPING THE FILE WRITING OPTION .. )
- $$$
- ( THE PROGRAM WILL SAY:)
-
- For what x you want y(x) ? (Hit Enter to stop) 1.1 (ENTER: 1.1 )
-
- ( THE PROGRAM WILL INTERPOLATE THE SOLUTION FOR X = 1.1 )
-
- x: 1.10000000000000E+000 y0: 8.91207360381939E-001
- y1: 4.53596121602802E-001
-
- ( CHECK IT OUT, THE TRUE SOLUTION IS:
-
- SIN(1.1) = 8.91207360061435E-001 SO THE ERROR IS: 3.21e-11
- COS(1.1) = 4.53596121425577E-001 SO THE ERROR IS: 1.77e-11 )
-
- ( THE PROGRAM WILL CONTINUE TO ASK FOR WHAT X TO INTERPOLATE, HIT ENTER TO
- QUIT ..)
-
- ( THE PROGRAM WILL GIVE YOU ONE MORE CHANCE TO WRITE THE SOLUTION INTO A FILE
- AS DESCRIBED IN EXAMPLE #1 ABOVE ..)
-
- That's all folks ...