RootsEqns: multiple nonlinear equation solver
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