(HP71B) Newton's method

[Albert Chan]

Inspired by Csaba Tizedes videos about Newton's method

(09-23-2023 05:19 PM)Csaba Tizedes Wrote:  https://www.youtube.com/shorts/1R3IMq227dM

(09-23-2023 05:21 PM)Csaba Tizedes Wrote:  Same, without d/dx()

https://www.youtube.com/shorts/-8pGxTXSk40

f'(x) ≈ (f(x+ε) - f(x-ε)) / (2ε)      → f/f' needed 3 function calls per iteration.

We may reduce this to about 2 calls per iteration.
Redo same example in HP71B, for full convergence.

Quote:1 DESTROY ALL @ C=0
2 DEF FNF(X)=X^4-4*X^3
3 INPUT "Guess, Accuracy = ";X,A
10 FOR I=1 TO 30 @ DISP X
20 E=.00001
30 F1=FNF(X-E) @ F2=FNF(X+E)
40 R=2*E/(F2-F1) ! reciprocal slope
50 D=(F1+F2)/2*R ! secant root = X-D
60 C=C+2 @ IF ABS(D)<E THEN D=D+FNF(X-D)*R @ C=C+1
70 X=X-D @ IF ABS(D)<A THEN EXIT I
80 NEXT I
90 DISP X @ DISP "FNF calls =";C

>RUN
Guess, Accuracy = 10, 1e-6
 10
 7.85714477037
 6.29727032561
 5.20041359757
 4.4911553721
 4.12133235435
 4.00984649554
 4.00007201239
 4.00000000386
 4
FNF calls = 19

For comparison, this is using video's formula for Newton's correction.
There is no gain in rate of convergence ... just more function calls.

>50 D = FNF(X)*R
>60 C = C+3
>RUN
Guess, Accuracy = 10, 1e-6
 10
 7.8571447704
 6.29727032563
 5.20041359761
 4.49115522316
 4.12133220071
 4.00984642245
 4.0000719974
 4.00000000388
 4
FNF calls = 27

[Albert Chan]

Previous HP71B code had a hard-coded absolute epsilon. (line 20)
It may be safer to use relative epsilon instead, in case x±ε round back to x

>20 E=(ABS(X)+(X==0))*1E-3

Let's try example from thread, Namir, Byte and REXX, which also does 3 calls per iteration.

(09-10-2018 01:26 PM)Massimo Gnerucci Wrote:  D:\ReginaREXX391w64>regina.exe root.rex
Enter expression (of variable X)
2 * X**2 + 3 * X - 12
Enter guess and accuracy
5 0.001
guess = 5
guess = 2.695652174
guess = 1.925113151
guess = 1.814140206
guess = 1.811738817
Root   = 1.811737692
Number of iterations = 5

>2 DEF FNF(X)=2*X^2+3*X-12
>RUN
Guess, Accuracy = 5, 1e-6
 5
 2.69565
 1.92511161178
 1.81413944914
 1.81173817456
 1.81173769148
FNF calls = 11

True root = (sqrt(105)-3)/4 = 1.8117376914898995958...

[Albert Chan]

(09-23-2023 09:46 AM)carey Wrote:  The simple looking transcendental equation, e^-x = x, is unsolvable in closed form (try it!)

Newton form of x = f(x):

0 = (f - x)
x = x - (f - x) / (f' - 1)  =  (x*f' - f) / (f' - 1)  =  (x - f/f') / (1 - 1/f')

W(b=1) example, newton iteration formula:  x = (x + 1) / (1 + e^x/b)

This setup is nice, with denominator > 0, if b > -1/e = -0.36787944...
Even crazy guess work. For W(b=1): ∞ --> 0 --> 0.5 --> 0.5663 --> ...

> 2 B=1 @ DEF FNF(X)=B*EXP(-X)-X ! root = W(B)
> RUN
Guess, Accuracy = B, 1E-6
 1
 .537882997929
 .566987045872
 .567143290413
 .56714329041
FNF calls = 10

Interestingly, crazy guess of 1e499 cost only 1 extra FNF call.
Note: this work only with previous post relative epsilon patch.

[Albert Chan]

(10-01-2023 09:47 PM)Albert Chan Wrote:  W(b=1) example, newton iteration formula:  x = (x + 1) / (1 + e^x/b)

Because of e^x term in denominator, this setup is not good for W(huge b)
Convegence is slow and inaccurate, even with good guess.

>2 B=1E10 @ DEF FNF(X)=B*EXP(-X)-X ! root = W(B)
>RUN
Guess, Accuracy = LN(B), RES*1E-6
 23.0258509299
 12.0135445721
 13.0133779957
 14.0128061672
 15.0110431755
 16.005828047
 16.9907276945
 17.9478731089
 18.8306249224
 19.5363990109
 19.9283981583
 20.0242147693
 20.0286862103
 20.0286854427
FNF calls = 28

>X*EXP(X) ! expected to be close to B
 10000000308.6

---

It is better is to setup Newton formula with less curvy f (e^-x is as bad as e^x)

x * e^x = b
ln(x) + x = ln(b)
x = f(x)      → f = ln(b) - ln(x)      → f' = -1/x

Newton: x = (x - f/f') / (1 - 1/f') = (1 + ln(b/x)) / (1 + 1/x)

>2 B=1E10 @ DEF FNF(X)=LN(X/B)+X ! root = W(B)
>RUN
Guess, Accuracy = LN(B), RES*1E-6
 23.0258509299
 20.0197856819
 20.0286854133
 20.0286854133
FNF calls = 8

>X*EXP(X) ! very close to B. no other 12-digits X can get closer!
 9999999999.95