29 December 2018

Numerical solver - Newton-Raphson method

The program downloadable from this page is a numerical solver. It attempts to find the zeroes of a function $f(x)$, i.e. it attempts to find $x$ such that $f(x)=0$

There are several methods for this. One that requires little computing power is Newton's Method, or the Newton-Raphson method, a description of which can be found on Wikipedia.

We have no way of knowing the exact value of $f'(x)$ on the HP-41CX used here because that machine has no symbolic math capabilities, so we calculate an approximation of $f'(x)$ by taking a small value $\epsilon$ and calculating:
$f'(x) \approx \frac {f(x+\epsilon)-f(x)} \epsilon$
This then allows us to calculate the next value of $x$ to use and display so that we see the calculator homing in on the zero that we're looking for:
$x_{n+1} = x_n - \frac {f(x_n)} {f'(x_n)}$
Once we find a value of $x$ such that $|f(x)| \lt \epsilon$, we stop and that value of $x$ is deemed to be a zero of $f(x)$.

The same $\epsilon$ is used here as the small value added to $x$ to calculate $f'(x)$ and as the threshold below which we deem $f(x)$ to be zero. It is calculated according to the display precision of the HP-41CX. The calculator's flags 36-39 are examined to ascertain how many significant figures are displayed and we retrieve the digit after the last used FIX, SCI or ENG command. E.g., if the calculator is in FIX 4 mode then we get a 4 back. Let's call this number $p$ for precision.

We then calculate $\epsilon = 10^{-p}$, so the higher the display precision, the greater precision we seek from the algorithm. Similarly, if less precision is required for the display then less precision is demanded from the algorithm and it can complete after fewer iterations.

It may well be that the algorithm is unable to converge on a zero because of the nature of the function studied. It could be constant or it could send the algorithm off on a wild goose chase. There's some discussion of such cases on the Wikipedia page linked to above.

Finally, the link to download this HP-41 program is here: newton-raphson.zip

It is provided in text listing format, a .raw file and a printable PDF with wand bar codes.

Here is a video of a HP-41CX finding various zeros of the function:
$f(x) = x^3-2x^2-11x+12 = (x+3)(x-1)(x-4)$
The three zeroes are therefore $x=-3, x=1$ and $x=4$.

And here is the equivalent in BASIC for HP-71:
10 DESTROY ALL @ DELAY 0
20 REAL X,Y,Y1,E @ INTEGER P,M
30 P=FLAG(-17)+2*FLAG(-18)+4*FLAG(-19)+8*FLAG(-20) @ E=10^(-P) @ M=2*P+8
40 INPUT "GUESS? ";X
50 Y=FNF(X) @ IF ABS(Y)<E THEN 120
60 "X=" & STR$(X) 70 Y1=(FNF(E+X)-FNF(X))/E 80 IF Y1=0 THEN "CONSTANT?" @ END 90 X=X-Y/Y1 100 IF M>0 THEN M=M-1 @ GOTO 50 110 "NO CONVERGENCE" @ END 120 "ZERO: " & STR$(X) @ END
130 DEF FNF(X) = (X+3)*(X-1)*(X-4)
After running, the variable X contains the zero found and Y is the value of $f(x)$ calculated for that value of $x$ - this gives you the error, the absolute value of which should be lower than our $\epsilon$, which is stored in the variable E.