Solving Simultaneous Nonlinear Equations for the HP-42S

The program 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 

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(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 R11
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 11
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.

BACK

Copyright (c) Namir Shammas. All rights reserved.