RootsEqns is a program that finds the roots of N (<= 10) nonlinear functions. The program uses a straightforward version of Newton's algorithm: x' = x - F(x) / J(x) Where x = array of guesses x' = array of refined guesses F(x) = array of nonlinear functions J(x) = Jacobian matrix of function derivatives, J(i,j) = dFi(x)/dxj The iterations continue until one of these conditions are met: 1. The maximum number of iterations is reached. 2. The expression ||x' - x|| / N^0.5 < XTol is true 3. The expression ||F(x)|| / N^0.5 < FTol is true USAGE: 1. To start the program execute XEQ "RTEQNS" 2. The program prompts you to enter the functions toleance value. 3. The program prompts you to enter the root-guesses toleance value. 4. The program prompts you to enter the maximum number of iterations. 5. The program displays the array of guesses in edit mode. Enter the values (you can inspect the various values and edit them). When you are done, press the R/S key. 6. The program displays the solution vector in edit mode. If the maximum number of iterations is exceeded you get a message first. Press R/S to view the current values in array of guesses. Note: The program divides the tolerance values by the square-root of the number of equations in order to get an average that works better with the norm of the arrays of guesses and function values.. To program a new set of equations perform the following edits: 1. Set line 2 to be equal to the number of equations. 2. Edit or cretae labels 01 to 10 (as needed) to represent the nonlinear functions. Label 01 calculates the value for function F1(x). Label 02 calculates the value for function F2(x), and so on 3, In labels 01 through 10 use registers 01 through 10 to represent x(1) through x(10). 4. Edit label 36 to make sure that all of the elements of array VX are copied to memory registers 01 and up. Here is an example of label 36 copying data to solve for 2 nonlinear equations: LBL 36 INDEX "VX" RCLEL STO 01 I+ RCLEL STO 02 RTN To copy more elements insert the following command pattern as many times needed: I+ RCLEL STO <register_number> After editing for a specific set of equations you may want to save the program under a distinct name. EXAMPLE 1: NOTE: You can load the file RootsEqns.raw to run the next example. Given the functions: F1(x) = x(1)^2 + x(2)^2 - 1 F2(x) = x(1)^2 - x(2)^2 + 0.5 Find the roots near x(1) = 1 and x(2) = 1. Allow 100 iterations at most and use a root tolerance of 1E-5.and function tolerance of 1E-8. The code for labels 01, 02, and 36 for this example are: LBL 01 XEQ 36 RCL 01 X^2 RCL 02 X^2 + 1 - RTN LBL 02 XEQ 36 RCL 01 X^2 RCL 02 X^2 - 0.5 + RTN LBL 36 INDEX "VX" RCLEL STO 01 I+ RCLEL STO 02 RTN Also make sure that line 2 contains the number 2 which is the number of nonlinear equations. XEQ "RTEQNS" XTOL= 1.0000E-5 1E-5 RUN FTOL= 1.0000E-8 1E-8 RUN MAX= 100.0000 100.0000 RUN X1? 1.0000 RUN X2? 1.0000 RUN VX= [ 2x1 Matrix ] 1:1= 0.5000 2:1= 0.8660 The program found roots at x(1) = 0.5 and x(2) = 0.866. EXAMPLE 2: NOTE: You can load the file RootsEqns3.raw to run the next example. Otherwise you need to perform the edits described in this example. Given the functions: F1(x) = x(1) + x(2) + x(3)^2 - 12 F2(x) = x(1)^2 - x(2) + x(3) - 2 F3(x) = 2 * x(1) - x(2)^2 + x(3) - 1 Make sure that line 2 contains the integer 3 (the number of nonlinear equations to solve), as shown in the next code snippet: 01 LBL "RTEQNS" 02 3 03 STO "N" Also edit labels 01, 02, 03, and 36 to be as follows: LBL 01 XEQ 36 RCL 01 RCL 02 + RCL 03 X^2 + 12 - RTN LBL 02 XEQ 36 RCL 01 X^2 RCL 02 - RCL 03 + 2 - RTN LBL 03 XEQ 36 RCL 01 2 × RCL 02 X^2 - RCL 03 + 1 - RTN LBL 36 INDEX "VX" RCLEL STO 01 I+ RCLEL STO 02 I+ RCLEL STO 03 RTN After performing the above edits, you are ready to use the program. Find the roots near x(1) = 0, x(2) = 0, and x(3) = 0. Allow 100 iterations at most and use a root tolerance of 1E-5.and function tolerance of 1E-8. Make sure that line 2 contains the number 2 which is the number of nonlinear equations. XEQ "RTEQNS" XTOL= 1.0000E-5 1E-5 RUN FTOL= 1.0000E-8 1E-8 RUN MAX= 100.0000 100.0000 RUN X1? 0.0000 RUN X2? 0.0000 RUN X3? 0.0000 RUN VX= [ 3x1 Matrix ] 1:1= -0.2337 2:1= 1.3532 3:1= 3.2986 The program found roots at x(1) = -0.2337, x(2) = 1.3532, and x(3) = 3.2986. Run the program again and provide an initial guesas of 1 for the three roots. Use the same values for the tolerances and maximum iteration limits. XEQ "RTEQNS" XTOL= 1.0000E-5 RUN FTOL= 1.0000E-8 RUN MAX= 100.0000 RUN X1? 1.0000 RUN X2? 1.0000 RUN X3? 1.0000 RUN VX= [ 3x1 Matrix ] 1:1= 1.0000 2:1= 2.0000 3:1= 3.0000 The program found roots at x(1) = 1, x(2) = 2, and x(3) = 3. PSEUDO-CODE: N = number of equations ' hard coded Redimention matrix A and arrays X, DX, and F Input XTol, FTol, MaxIter, and array X Iter = 0 Do Increment Iter If Iter > MaxIter then exit For I = 1 to N F(I) = FX(X, I) FF = F(I) For J = 1 to N A(I,J) = FF Next J Next I For I = 1 to N XX = X(I) H = 0.01 * (1 + ABS(XX)) X(I) = XX + H For J = 1 to N FF = FX(X, J) A(J,I) = (FF - A(J,I)) / H ' Jacobian matrix elements Next J X(I) = XX Next I DX = INVERSE(A) * F ' This is a matrix operation X = X - DX ' This is a matrix operation Loop Until ||DX|| / N^0.5 < XTol Or ||F|| / N^0.5 < FTol Display elements of array X MEMORY USAGE: MA= Matrix of function slopes VB= Array of function values VX= Array of root guesses VDX= Array of root guesses refinements FTOL= Functions tolerance XTOL= Tolerance for root refinements N= Number of functions I= Loop control variable J= Loop control variable FF= Used XX= Used H= Increment in a guess value MAX= Maximum number ot iterations ITR= Number of iterations R00 to R10 are used to store X(1) through X(10) R11= Number of digits displayed R12= Square root of the number of equations FLAGS 00 Display control 01 Display control 02 Control of indexing array MA LISTING: LBL "RTEQNS" 2 # Number of equations. MUST BE CHANGED ACCORDING TO ACTUAL PROBLEM STO "N" # Preset number of equations ENTER DIM "MA" # Create and dimention the matrix and arrays 1 DIM "VB" DIM "VX" DIM "VDX" INPUT "XTOL" INPUT "FTOL" INPUT "MAX" RCL "N" # Calculate square root of N SQRT STO 12 # Store in R12 0 STO "ITR" # Initialize number of iterations XEQ 17 # Store display setting XEQ 20 # Set I = 1 to N LBL 19 # Start loop to prompt for array x ---------------------- 19 "X" FIX 00 CF 29 ARCL "I" |- "?" XEQ 18 # Restore display SF 29 PROMPT XEQ 30 # Store input in X(I) ISG "I" GTO 19 # End of loop -------------------------------------------<19> LBL 00 # Main loop ---------------------------------------------- 00 1 STO+ "ITR" RCL "ITR" RCL "MAX" X<Y? GTO 16 # Solution diverged -------------------------------------<16> XEQ 20 # Set I = 1 to N LBL 11 # Outer loop --------------------------------------------- 11 XEQ IND "I" # Calculate function F<I>(x) STO "FF" XEQ 34 # Set F(I) = FF XEQ 21 # Set J = 1 to N CF 02 LBL 12 # Inner loop --------------------------------------------- 12 RCL "FF" XEQ 32 # Set A(I,J) = FF ISG "J" GTO 12 # End of inner loop --------------------------------------<12> ISG "I" GTO 11 # End of outer loop --------------------------------------<11> XEQ 20 # Set I = 1 to N LBL 13 # Outer loop --------------------------------------------- 13 XEQ 31 # Get X(I) STO "XX" # Set XX = X(I) ABS 1 + 0.01 x STO "H" # Calculate and store H RCL "XX" + XEQ 30 # X(I) = X(I) + H XEQ 21 # Set J = 1 to N SF 02 LBL 14 # Inner loop --------------------------------------------- 14 XEQ IND "J" STO "FF" # Calculate FF = F<J>(x) XEQ 33 # Get A(J,I) RCL "FF" - +/- RCL "H" / XEQ 32 # Store A(J,I) = derivative of function I wrt X(J) ISG "J" GTO 14 # End of inner loop --------------------------------------<14> CF 02 RCL "XX" XEQ 30 # Restore X(I) = XX ISG "I" GTO 13 # End of outer loop --------------------------------------<13> RCL "MA" INVRT RCL "VB" x STO "VDX" STO- "VX" # Set array x = array x - array delta x FNRM # Calculate vector norm RCL 12 / RCL "XTOL" X>Y? GTO 15 # ---------------------------------------------------------<15> RCL "VB" FNRM # Calculate vector norm RCL 12 / RCL "FTOL" X>Y? GTO 15 # ---------------------------------------------------------<15> GTO 00 # End of main loop ----------------------------------------<00> LBL 16 # Handle case when exceeded max iterations ---------------- 16 "ITERS EXCEEDED" PROMPT LBL 15 CLV "H" # Clear local variables CLV "XX" CLV "FF" CLV "I" CLV "J" PRV "VX" # Print the array of guesses RCL "VX" # Recall array of guesses onto the stack EDIT # View array of guesses BEEP RTN LBL 20 # Set I = 1 to N ----------------------------------------- 20 1 RCL "N" 1E3 / + STO "I" RTN LBL 21 # Set J = 1 to N ----------------------------------------- 21 1 RCL "N" 1E3 / + STO "J" RTN LBL 30 # Store X(I) --------------------------------------------- 30 INDEX "VX" RCL "I" 1 STOIJ RCL ST Z STOEL RTN LBL 31 # Recall X(I) -------------------------------------------- 31 INDEX "VX" RCL "I" 1 STOIJ RCLEL RTN LBL 32 # Store A(I,J) ------------------------------------------- 32 INDEX "MA" RCL "I" RCL "J" FS? 02 X<>Y STOIJ RCL ST Z STOEL RTN LBL 33 # Recall A(I,J) ------------------------------------------ 33 INDEX "MA" RCL "I" RCL "J" FS? 02 X<>Y STOIJ RCLEL RTN LBL 34 # Store F(I) --------------------------------------------- 35 INDEX "VB" RCL "I" 1 STOIJ RCL ST Z STOEL RTN LBL 35 # Recall F(I) -------------------------------------------- 35 INDEX "VB" RCL "I" 1 STOIJ RCLEL RTN LBL 17 # Switch display mode to showing integers ---------------- 17 CF 00 # Determine display mode (FIX, ENG, or SCI) CF 01 FS? 40 SF 00 FS? 41 SF 01 0 # Determine the number of digits displayed STO 11 1 FS? 39 STO+ 11 2 FS? 38 STO+ 11 4 FS? 37 STO+ 11 8 FS? 36 STO+ 11 # Store number of digits in R08 RTN LBL 18 # Restore display mode --------------------------------------- 18 SF 29 # Show decimal point SCI IND 11 # Set to SCI mode by default FS? 00 # Was it in FIX mode? FIX IND 11 FS? 01 # Was it in ENG mode? ENG IND 1 RTN LBL 01 # Calculate function F1(x) = x(1)^2 + x(2)^2 - 1 XEQ 36 # Copy array X into registers RCL 01 X^2 RCL 02 X^2 + 1 - RTN LBL 02 # Calculate function F2(x) = x(1)^2 - x(2)^2 + 0.5 XEQ 36 # Copy array X into registers RCL 01 X^2 RCL 02 X^2 - 0.5 + RTN LBL 36 # Copy array X in registers 01 and up ----------------------- 36 INDEX "VX" RCLEL STO 01 I+ RCLEL STO 02 RTN .END. LABEL USAGE: 01 Label for function F1(x) 02 Label for function F2(x) ... 10 Label for function F10(x) 00 Main loop 10 loop 11 loop 12 loop 13 loop 14 loop 15 used 16 used 17 Query display mode and digits 18 Restore display mode 20 Set I = 1 + N / 1000 21 Set J = 1 + N / 1000 30 Store value in X(I) 31 Recall X(I) 32 Store value in A(I,J) if flag 2 is clear or in A(J,I) is flag 2 is set 33 Recall A(I,J) if flag 2 is clear or recall A(J,I) is flag 2 is set 34 Store value in F(I) 35 Recall F(I) 36 Copy values of array X into memory registers 01 and up. Namir Shammas Program Version 1.0.0 Documentation Version 1.1.0 June 1, 2005 nshammas@aol.com