(HP71B) Newton's method

[Albert Chan]

Instead of calculating f(x±h), we can do f(x+h), f(x+h2), with x+h, x+h2 closer to true root.

From points {x, h=fnh(x)}, {x+h, h2=fnh(x+h)}, we solve for secant root
Note: h = fnh(x) is simply scaled f(x), so that x + h is closer to true root.

x - ((x+h)-x) / (fnh(x+h)-fnh(x)) * h = x + (-h/(h2-h)) * h = x + s * h

Old s can be used to improve fnh(x), to get closer to true root (*)

lua> f = fn'x:exp(x)-3*x*x'
lua> f(3), f(4)
-6.914463076812332      6.598150033144236

f root = 3 .. 4. For h, we need to scale down f(x), and flip its direction.

10 DEF FNF(X)=EXP(X)-3*X*X
20 INPUT "X,S = ";X,S @ H=X @ C=0
30 FOR I=1 TO 30 @ H=H/X @ IF 1+H*H=1 THEN EXIT I
40 H=FNF(X)*S @ H2=FNF(X+H)*S @ C=C+2
50 IF H2 THEN S2=H/(H-H2) @ H=H*S2 @ S=S*S2
60 X=X+H @ DISP X,H,C @ NEXT I

Code:

X,S = 3, -0.1
 3.77904928733        .77904928733         2
 3.73056592227       -4.84833650565E-2     4
 3.73307897801        2.51305573645E-3     6
 3.73307902862        5.06081304277E-8     8

Code:

X,S = 4, -0.1
 3.63243460816       -.367565391835        2
 3.7332302642         .100795656045        4
 3.73307902653       -1.5123767276E-4      6
 3.73307902862        2.09164536676E-9     8

Here is previous post 'flattened' FNF. We just need to switch direction.

>10 DEF FNF(x) = 1 - 3*x^2/exp(x)

Code:

X,S = 3, -1
 3.73923735544        .739237355436        2
 3.73307901812       -6.15833732419E-3     4
 3.73307902863        1.05123533817E-8     6

Code:

X,S = 4, -1
 3.7289957962        -.271004203799        2
 3.73307908199        4.08328579414E-3     4
 3.73307902862       -5.33668310581E-8     6

(*) I had made a mistake replacing old s with new s (negated reciprocal slope)

It should have been s = product of previous s's, and should converge to 1 value.
Even if secant root is off, secant line slope should be close, and can be reused.

Comment: This simple algorithm convergence rate is faster than Newton!
For both examples, rate of convergence approaches (1+√2) ≈ 2.41421
Distributed per function call, we have √(1+√2) ≈ 1.55377
This is close to Secant's rate, φ = (1+√5)/2 ≈ 1.61803

Per iteration convergence rate calculations, mpmath with mp.dps = 10000
Numbers are log ratio of |relative errors| between iteration

lua> ok = 3.73307902863
lua> log(abs(1-3.77904928733/ok)) / log(abs(1-3/ok)) -- 1st col, 1st row
2.701295751437436

Code:

f(x) = exp(x) - 3*x*x           f(x) = 1 - 3*x*x/exp(x)
x,s = 3,-0.1  x,s = 4,-0.1      x,s = 3,-1    x,s = 4,-1
2.70130       1.36973           3.93626       2.58454
1.66101       2.79900           3.07258       2.64935
2.48047       2.10556           2.45092       2.51723
2.26705       2.36375           2.45940       2.45177
2.38094       2.37393           2.42749       2.43010
2.39471       2.40055           2.42055       2.42066
2.40703       2.40795           2.41668       2.41689
2.41106       2.41171           2.41526       2.41532
2.41294       2.41316           2.41464       2.41467