Variant of the Secant Method - Printable Version

Variant of the Secant Method - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Variant of the Secant Method (/thread-16025.html)

Variant of the Secant Method - ttw - 12-09-2020 04:33 AM

The order of the secant method for finding zeros of a real function can be improved by using spline-like interpolants. This gives a method which has order just under 2 and uses only 1 function evaluation per step.

https://arxiv.org/pdf/2012.04248.pdf

RE: Variant of the Second Method - Valentin Albillo - 12-09-2020 05:40 AM

(12-09-2020 04:33 AM)ttw Wrote: The order of the secant method for finding zeros of a real function can be improved by using spline-like interpolants. This gives a method which has order just under 2 and uses only 1 function evaluation per step.

https://arxiv.org/pdf/2012.04248.pdf

Interesting. Perhaps you might consider implementing it yourself for some HP calc and telling us how it compares in terms of speed and memory versus plain old Newton's method {two evaluations of f(x) per iteration), which can be implemented in as few as 34 steps even in the simple HP-10C model (fewer in RPN models with subroutine capability like the HP-11C).

Anyway, thanks for sharing and take care.
V.

RE: Variant of the Second Method - Albert Chan - 12-09-2020 04:38 PM

Let's try solving x, for sinc(x) = 1/6.5

lua> function f(x) return sin(x)/x - 1/6.5 end
lua> function secant(x1,y1,x2,y2) return x1-y1/(y2-y1)*(x2-x1) end
lua> x1, x2 = 2.711, 2.712
lua> y1, y2 = f(x1) , f(x2)
lua> x3 = secant(x1,y1, x2,y2)
lua> y3 = f(x3)
lua> x3
2.7113129474671775

-- Inverse Interpolation, assumed f(x) is strictly increasing or descreasing around root
-- Fit x = g(y), then evaluate g(0)

lua> x12 = (x2-x1)/(y2-y1)
lua> x23 = (x3-x2)/(y3-y2)
lua> x123 = (x23-x12)/(y3-y1)
lua> x3 + (0-y3)*(x23 + (0-y2)*x123)
2.711312910356982

-- Aikten interpolation, as 1 liner: inner secant = x4, outer secant = x5
-- y1 x1
-- y2 x2 x3
-- y3 x3 x4 x5

lua> secant(secant(x3,y3, x1,y1),y3, x3,y2)
2.711312910356982

-- improved secant, interpolate for slope

lua> y23 = (y3-y2)/(x3-x2)
lua> y12 = (y2-y1)/(x2-x1)
lua> y123 = (y23-y12)/(x3-x1)
lua> slope = y23 + (x3-x2)*y123
lua> x3 - y3/slope
2.7113129103569826

-- Newton's method, actual slope

lua> x, y = x3, y3
lua> slope = (x*cos(x) - sin(x))/x^2
lua> x - y/slope
2.711312910356983

For reference, I used ASINC, for Free-42 Decimal.
Comparing apples to apples, input is float(1/6.5) = 0x1.3b13b13b13b14p-3

5542891849071380 2 55 X^Y ÷ XEQ "ASINC"

2.71131291035698331469084928874754

All 3 methods gives excellent results, doubling accuracy of x3

RE: Variant of the Secant Method - Albert Chan - 12-09-2020 11:51 PM

I extend previous post, with another round of iteration, using all points.
(First iteration matched previous post, done in double precision)

inverse-interpolation:
2.71131291035698220988289207111279
2.71131291035698331469084928874683

improved secant:
2.71131291035698244759230088334763
2.71131291035698331469084928874750

newton's method:
2.71131291035698307703744887554087
2.71131291035698331469084928874753

RE: Variant of the Secant Method - Namir - 12-10-2020 01:54 PM

I implemented the Generalized Secant algorithm in Excel VBA. I used the same order as in the article's example (divided differences for first and second derivatives). The algorithm did better than Newton's method in most cases. The Generalized Secant algorithm needs two initial guesses, while Newton's method needs one. It made a difference in some test if I gave Newton's method the second guess OR averaged the two guesses. The closer of these two values to the neighboring root gave better results.

The algorithm is very interesting!!!

Thanks!

Namir

RE: Variant of the Secant Method - Namir - 12-10-2020 02:56 PM

I was able to implement in Excel VBA a version of the Generalized Secant algorithm using divided differences for the first, second and third derivatives.

Here is the VBA listing for the version that uses two initial guesses and divided differences for the first and second derivatives.

Code:

Cell Mapping

===========

A1 : "X0"

A2 : value for X0

A3 : "X1"

A4 : value for X1

A5 : "Tolerance"

A6 : value for tolerance

A7 : "F(x)"

A8: String containing expression for f(x)=0, such as exp(x)-3*x^2.

B1 : "X"

C1 : "F(X)"  (for Generalized Secant method)

D1 : "X"

E1 : "F(X)"  (for Newton's method, used for comparison)

VBA Listing

===========

Function Fx(ByVal sFx As String, ByVal X As Double) As Double

  sFx = Replace(sFx, "EXP", "!")

  sFx = Replace(sFx, "X", "(" & CStr(X) & ")")

  sFx = Replace(sFx, "!", "EXP")

  Fx = Evaluate(sFx)

End Function

Sub genSecant2()

  Dim sFx As String, X0 As Double, X1 As Double, X2 As Double, X3 As Double

  Dim Toler As Double, Df1 As Double, Df2 As Double

  Dim Row As Integer, h As Double

  Dim Fx0 As Double, Fx1 As Double, Fx2 As Double, Diff As Double

  Range("B2:X1000").Clear

  X0 = [A2].Value

  X1 = [A4].Value

  Toler = [A6].Value

  sFx = [A8].Value

  sFx = Replace(sFx, " ", "")

  sFx = UCase(sFx)

  Fx0 = Fx(sFx, X0)

  Fx1 = Fx(sFx, X1)

  Df1 = Fx0 / (X0 - X1) + Fx1 / (X1 - X0)

  X2 = X1 - Fx1 / Df1

  Fx2 = Fx(sFx, X2)

  Diff = X2 - X1

  Row = 2

  Cells(Row, 2) = X2

  Cells(Row, 3) = Fx(sFx, X2)

  Row = Row + 1

  Do Until Abs(Diff) < Toler

    Df1 = Fx2 / (X2 - X1) + Fx1 / (X1 - X2)

    Df2 = Fx2 / (X2 - X1) / (X2 - X0) + Fx1 / (X1 - X2) / (X1 - X0) + Fx0 / (X0 - X2) / (X0 - X1)

    Diff = Fx2 / (Df1 + Df2 * (X2 - X1))

    X3 = X2 - Diff

    X0 = X1

    Fx0 = Fx1

    X1 = X2

    Fx1 = Fx2

    X2 = X3

    Fx2 = Fx(sFx, X3)

    Cells(Row, 2) = X3

    Cells(Row, 3) = Fx2

    Row = Row + 1

  Loop

  ' Newton's method

  X1 = [A4].Value

  Fx1 = Fx(sFx, X1)

  Row = 2

  Diff = 2 * Toler

  Do Until Abs(Diff) < Toler

    h = 0.01 * (1 + Abs(X1))

    Diff = h * Fx1 / (Fx(sFx, X1 + h) - Fx1)

    X1 = X1 - Diff

    Fx1 = Fx(sFx, X1)

    Cells(Row, 4) = X1

    Cells(Row, 5) = Fx1

    Row = Row + 1

  Loop

End Sub

Output

======

Columns B and C show the refined guesses and their function values for the Generalized Secant Algorithm.

Columns D and E show the refined guesses and their function values for Newton's method.

Here is the version of the VBA code that uses three initial guesses and the divided differences for the first three derivatives:

Code:

Cell Mapping

===========

A1 : "X0"

A2 : value for X0

A3 : "X1"

A4 : value for X1

A5 : "X2"

A6 : value for X2

A7 : "Tolerance"

A8 : value for tolerance

A9 : "F(x)"

A10: String containing expression for f(x)=0

B1 : "X"

C1 : "F(X)"  (for Generalized Secant method)

D1 : "X"

E1 : "F(X)"  (for Newton's method, used for comparison)

VBA Listing

===========

Function Fx(ByVal sFx As String, ByVal X As Double) As Double

  sFx = Replace(sFx, "EXP", "!")

  sFx = Replace(sFx, "X", "(" & CStr(X) & ")")

  sFx = Replace(sFx, "!", "EXP")

  Fx = Evaluate(sFx)

End Function

Sub genSecant3()

  Dim sFx As String, X0 As Double, X1 As Double, X2 As Double, X3 As Double, X4 As Double

  Dim Toler As Double, Df1 As Double, Df2 As Double, Df3 As Double

  Dim Row As Integer, h As Double

  Dim Fx0 As Double, Fx1 As Double, Fx2 As Double, Fx3 As Double, Diff As Double

  Range("B2:X1000").Clear

  X0 = [A2].Value

  X1 = [A4].Value

  X2 = [A6].Value

  Toler = [A8].Value

  sFx = [A10].Value

  sFx = Replace(sFx, " ", "")

  sFx = UCase(sFx)

  Fx0 = Fx(sFx, X0)

  Fx1 = Fx(sFx, X1)

  Fx2 = Fx(sFx, X2)

  Df1 = Fx0 / (X0 - X1) + Fx1 / (X1 - X0)

  Df2 = Fx2 / (X2 - X1) / (X2 - X0) + Fx1 / (X1 - X2) / (X1 - X0) + Fx0 / (X0 - X2) / (X0 - X1)

  X3 = X2 - Fx2 / (Df1 + Df2 * (X2 - X1))

  Fx3 = Fx(sFx, X3)

  Diff = X3 - X2

  Row = 2

  Cells(Row, 2) = X3

  Cells(Row, 3) = Fx3

  Row = Row + 1

  Do Until Abs(Diff) < Toler Or Row > 100

    Df1 = Fx3 / (X3 - X2) + Fx2 / (X2 - X3)

    Df2 = Fx3 / (X3 - X1) / (X3 - X2) + Fx2 / (X2 - X1) / (X2 - X3) + Fx1 / (X1 - X2) / (X1 - X3)

    Df3 = Fx3 / (X3 - X2) / (X3 - X1) / (X3 - X0) + _

          Fx2 / (X2 - X3) / (X2 - X1) / (X2 - X0) + _

          Fx1 / (X1 - X3) / (X1 - X2) / (X1 - X0) + _

          Fx0 / (X0 - X3) / (X0 - X2) / (X0 - X1)

    Diff = Fx3 / (Df1 + Df2 * (X3 - X2) + Df3 * (X3 - X2) * (X3 - X1))

    X4 = X3 - Diff

    X0 = X1

    Fx0 = Fx1

    X1 = X2

    Fx1 = Fx2

    X2 = X3

    Fx2 = Fx3

    X3 = X4

    Fx3 = Fx(sFx, X4)

    Cells(Row, 2) = X4

    Cells(Row, 3) = Fx3

    Row = Row + 1

  Loop

  ' Newton's method

  X1 = [A4].Value

  Fx1 = Fx(sFx, X1)

  Row = 2

  Diff = 2 * Toler

  Do Until Abs(Diff) < Toler

    h = 0.01 * (1 + Abs(X1))

    Diff = h * Fx1 / (Fx(sFx, X1 + h) - Fx1)

    X1 = X1 - Diff

    Fx1 = Fx(sFx, X1)

    Cells(Row, 4) = X1

    Cells(Row, 5) = Fx1

    Row = Row + 1

  Loop

End Sub

Output

======

Columns B and C show the refined guesses and their function values for the Generalized Secant Algorithm.

Columns D and E show the refined guesses and their function values for Newton's method.

Use different worksheets for each of the above methods as the input values in column A are not the same.

I wrote the above code to make it easy to ready. A programmer can optimize the above code by using arrays, matrices, and loops.

RE: Variant of the Secant Method - Albert Chan - 12-11-2020 04:34 PM

Quadratic interpolation is simply 2 linear interpolations, followed by linear interpolation.

x3 y3
x2 y2 y23
x1 y1 y13 y123

Similarly, cubic interpolation is 2 quadratic interpolations, followed by linear interpolation.

x4 y4
x3 y3
x2 y2 y234
x1 y1 y134 y1234

Using this recursive nature of polynomial interpolation, I get this neat Lua code Smile

Code:

function interp(x1,y1,x2,y2,x3,y3,x4,y4)    -- estimate y = p(0)

    if x3 then

        y2 = interp(x2,y2,x3,y3,x4,y4)

        y1 = interp(x1,y1,x3,y3,x4,y4)

    end

    local y, p = y2, x2/(x2-x1)

    if p > 0.5 then y, p = y1, x1/(x2-x1) end

    return y - p*(y2-y1)

end

Extending interp(...) is trivial, by adding (x5,y5), ...
The down-side is each bump up in polynomial degree, doubled the work.
For high degree fit, divided-difference style is better, because it can cache previous results.

OTTH, assuming we have convergence, old points are very, very far away.
Fitting those probably may not help much ...

---

For slope of improved secant, we can also consider this an interpolation problem.
Example, for a cubic fit, we generate divided-difference table:

x4 y4
x3 y3 y34
x2 y2 y24 y234
x1 y1 y14 y123 y1234

p(x) = y4 + (x-x4)*(y34 + (x-x3)*(y234 + (x-x2)*y1234)) = y4 + (x-x4)*M

p'(x) = M + (x-x4)*M' → p'(x4) = M

Basing all values relative to (x4, y4). Example X3 = x3-x4, Y3 = y3-y4 ...
M is just interpolation of individual slopes.

X3 Y3/X3
X2 Y2/X2
X1 Y1/X1 M

Again, extending this is simply adding more slopes to fit Smile

Code:

function secant(x1,y1,x2,y2,x3,y3,x4,y4)    -- estimate root of p(x)

    if not x3 then return interp(y1,x1,y2,x2) end

    if x4 then x3,y3,x4,y4 = x4,y4,x3-x4,(y3-y4)/(x3-x4) end

    x1, y1, x2, y2 = x1-x3, y1-y3, x2-x3, y2-y3

    return x3 - y3 / interp(x1,y1/x1, x2,y2/x2, x4,y4)    

end

Example:

lua> function f(x) return sin(x)/x - 1/6.5 end
lua> x1, x2 = 2.711, 2.712
lua> y1, y2 = f(x1), f(x2)
lua> x3 = secant(x1,y1, x2,y2)
lua> y3 = f(x3)

lua> interp(y1,x1, y2,x2, y3,x3) -- inverse interpolate, x = g(0)
2.711312910356982
lua> secant(x1,y1, x2,y2, x3,y3) -- improved secant, x = x3 - y3/slope
2.7113129103569826

RE: Variant of the Secant Method - Albert Chan - 12-11-2020 04:46 PM

Let's compare inverse-interpolation vs improved-secant (code from previous post)

Code:

function test(f, x0, x1, n)   -- n = 1 (linear) to 3 (cubic)

    local y0, y1 = f(x0), f(x1)

    local t0 = {y0,x0, y1,x1} -- inverse interpolate

    local t1 = {x0,y0, x1,y1} -- improved secant

    n = (n or 3)*2 + 1

    for i = 4,18,2 do         -- i = table length

        x0 = interp(unpack(t0, math.max(1,i-n)))

        x1 = secant(unpack(t1, math.max(1,i-n)))

        t0[i+1], t0[i+2] = f(x0), x0

        t1[i+1], t1[i+2] = x1, f(x1)

        print(i/2, x0, x1)

    end

end

lua> function f(x) return x^3 - 8 end

lua> test(f, 5, 4, 2) -- https://arxiv.org/pdf/2012.04248.pdf, Table 2
2 3.081967213114754 3.081967213114754
3 2.3945608933295457 2.2862188297178117
4 2.0980163342039635 2.010344209437878
5 2.0088865345029276 1.99979593345267
6 2.0001129292461193 2.0000000722313933
7 2.0000000388749792 2.000000000000015
8 2.0000000000000164 2
9 2 2

Last column is improved secant, fitting data up to quadratic (slope by linear interpolation)
Middle column is inverse-interpolation for x, at y = 0, also fitting up to quadratic.

Improved secant uses the last point as the "base", then does correction.
Inverse-interpolation is not as good, probably due to putting equal weight to its points.

Let's try again, but with up to cubic fit.

lua> test(f, 5, 4, 3)
2 3.081967213114754 3.081967213114754
3 2.3945608933295457 2.2862188297178117
4 2.082744738020821 2.034337291023909
5 2.0058425366743506 2.000576313421517
6 2.0000383915262443 2.000000166004798
7 2.0000000023200664 2.0000000000000138
8 1.9999999999999998 2
9 2 2

Based on above, anything more than a cubic fit probably not worth the trouble.
For example, lets see how Newton's method does (actual slope)

lua> function df(x) return 3*x*x end
lua> x = 4
lua> for i=2, 9 do x = x - f(x)/df(x); print(i, x) end
2 2.833333333333333
3 2.221068819684737
4 2.0212735368091126
5 2.0002231146078984
6 2.0000000248863623
7 2.0000000000000004
8 2
9 2

RE: Variant of the Secant Method - Namir - 12-12-2020 04:22 AM

I went back to my Excel VBA code and added statements to implement Halley's method. The Generalized Secant method did better than Newton's method but ws behind Halley's method in the number of steps. Of course results depend on the function tested and the initial guess.

Namir

RE: Variant of the Secant Method - ttw - 12-12-2020 02:41 PM

Methods with high convergence can have problems numerically; they also tend to have a smaller domain of convergence. Multi-point methods seem to be a bit more stable than high order methods.

RE: Variant of the Secant Method - Albert Chan - 12-12-2020 04:27 PM

(12-12-2020 02:41 PM)ttw Wrote: Methods with high convergence can have problems numerically; they also tend to have a smaller domain of convergence.

This might not always true.
In fact, recently we come across opposite case, LambertW function implementation.

Solving f(x) = x*e^x - a, with Newton's method is very unstable.

Higher convergence method, say Halley, it becomes "safe"

P.S. the next post, I had suggested another way, to suppress the exponential.
This version, Newton's method is fast, and safe.

>>> from mpmath import *
>>> a = -1e99
>>> findroot(lambda x: x + ln(x/a), 200+3j) # bad guess
mpc(real='222.55067066150301', imag='3.1275404175517207')
>>> lambertw(a)
mpc(real='222.55067066150301', imag='3.127540417551721')

RE: Variant of the Secant Method - Albert Chan - 12-11-2023 03:25 PM

(12-11-2020 04:34 PM)Albert Chan Wrote: Quadratic interpolation is simply 2 linear interpolations, followed by linear interpolation.

x3 y3
x2 y2 y23
x1 y1 y13 y123

It may seems we need 3 secant root steps for 3-points fit, in a loop, we only need 2

x2 y2
x3 y3 y23
x1 y1 y12 y123

Next iteration y23 → y12, can be reused.

Secant setup code

Code:

function newx(a,fa, b,fb) return b - fb/(fb-fa)*(b-a) end

function secant_init(f,a,b,eps,verbal)

    eps = eps or 1e-9

    local f2 = verbal and function(x) print(x); return f(x) end or f

    local fa = f2(a); if fa==0 then return nil, a end

    local fb = f2(b); if fb==0 then return nil, b end

    if signbit(eps) and abs(fa)<abs(fb) then

        a,fa, b,fb = b,fb, a,fa

        if verbal then print('-- points swapped') end

    end

    return f2, a,fa, b,fb, abs(eps)

end

Code:

function S.secant2(f,a,b,eps,verbal,c)  -- 3 pt fit

    local fa,fb,fc,ab,bc

    f, a,fa, b,fb, eps = secant_init(f,a,b,eps,verbal)

    if not f then return a,a end

    ab = newx(a,fa, b,fb)

    c = c or ab

    fc = f(c); if fc==0 then return c,c end

    while abs(b-c) > eps do

        bc, ab = ab, newx(b,fb, c,fc)

        a,fa, b,fb, c = b,fb, c,fc, newx(bc,fa, ab,fc)

        fc = f(c); if fc==0 then return c,c end

    end

    return finite(newx(b,fb, c,fc)) or (b+c)/2

end

lua> f = fn'x: x^3 - 8'
lua> x = S.secant2(f, 5, 4, 1e-9, true)
5
4
3.081967213114754
2.3945608933295457
2.0980163342039635
2.0088865345029276
2.0001129292461193
2.0000000388749792
2.0000000000000164
2

Result matched post #8

For completeness, here is S.secant3 (OP "improved" secant)

Code:

function S.secant3(f,a,b,eps,verbal,c)  -- slope extrapolated

    local fa,fb,fc

    f, a,fa, b,fb, eps = secant_init(f,a,b,eps,verbal)

    if not f then return a,a end

    c = c or newx(a,fa, b,fb)

    fc = f(c); if fc==0 then return c,c end

    while abs(b-c) > eps do

        local ca,cb = c-a,c-b

        local slope = newx((fc-fa)/ca,ca, (fc-fb)/cb,cb)

        a,fa, b,fb, c = b,fb, c,fc, c-fc/slope

        fc = f(c); if fc==0 then return c,c end

    end

    return finite(newx(b,fb, c,fc)) or (b+c)/2

end

lua> x = S.secant3(f, 5, 4, 1e-9, true)
5
4
3.081967213114754
2.2862188297178117
2.010344209437878
1.99979593345267
2.0000000722313933
2.000000000000015
2

RE: Variant of the Secant Method - Thomas Klemm - 12-16-2023 08:11 PM

These functions implement the described secant method in Python for orders \(k \in \{1, 2, 3\}\):

Code:

def f(x):

    print(x)

    return x**3 - 8

def secant_1(x_4, x_5, ε):

    f_4 = f(x_4)

    f_5 = f(x_5)

    f_45 = (f_4 - f_5) / (x_4 - x_5)

    while ε < abs(f_5):

        x_6 = x_5 - f_5 / f_45

        f_6 = f(x_6)

        f_56 = (f_5 - f_6) / (x_5 - x_6)

        x_5, f_5, f_45 = x_6, f_6, f_56

    return x_5

def secant_2(x_4, x_5, ε):

    f_4 = f(x_4)

    f_5 = f(x_5)

    f_45 = (f_4 - f_5) / (x_4 - x_5)

    x_6 = x_5 - f_5 / f_45

    f_6 = f(x_6)

    f_56 = (f_5 - f_6) / (x_5 - x_6)

    f_456 = (f_45 - f_56) / (x_4 - x_6)

    while ε < abs(f_6):

        x_7 = x_6 - f_6 / (f_56 + (x_6 - x_5) * f_456)

        f_7 = f(x_7)

        f_67 = (f_6 - f_7) / (x_6 - x_7)

        f_567 = (f_56 - f_67) / (x_5 - x_7)

        x_5, x_6 = x_6, x_7

        f_6, f_56, f_456 = f_7, f_67, f_567

    return x_6

def secant_3(x_4, x_5, ε):

    f_4 = f(x_4)

    f_5 = f(x_5)

    f_45 = (f_4 - f_5) / (x_4 - x_5)

    x_6 = x_5 - f_5 / f_45

    f_6 = f(x_6)

    f_56 = (f_5 - f_6) / (x_5 - x_6)

    f_456 = (f_45 - f_56) / (x_4 - x_6)

    x_7 = x_6 - f_6 / (f_56 + (x_6 - x_5) * f_456)

    f_7 = f(x_7)

    f_67 = (f_6 - f_7) / (x_6 - x_7)

    f_567 = (f_56 - f_67) / (x_5 - x_7)

    f_4567 = (f_456 - f_567) / (x_4 - x_7)

    while ε < abs(f_7):

        x_8 = x_7 - f_7 / (f_67 + (x_7 - x_6) * (f_567 + f_4567 * (x_7 - x_5)))

        f_8 = f(x_8)

        f_78 = (f_7 - f_8) / (x_7 - x_8)

        f_678 = (f_67 - f_78) / (x_6 - x_8)

        f_5678 = (f_567 - f_678) / (x_5 - x_8)

        x_5, x_6, x_7 = x_6, x_7, x_8

        f_7, f_67, f_567, f_4567 = f_8, f_78, f_678, f_5678

    return x_7

We can run them using the Decimal library:

Code:

from decimal import Decimal, getcontext

getcontext().prec = 34

Code:

secant_1(Decimal(5), Decimal(4), 1e-32)

5
4
3.081967213114754098360655737704918
2.519552120040923041946117994770945
2.180972989759050190092856538085551
2.037953100909517790045306042203723
2.003198489980016115117431120840932
2.000059872823468592338193344700047
2.000000095647401657566635004047756
2.000000000002863282761488899959229
2.000000000000000000136932773807772
2.000000000000000000000000000000196
2.000000000000000000000000000000000

Code:

secant_2(Decimal(5), Decimal(4), 1e-32)

5
4
3.081967213114754098360655737704918
2.286218829717811307322668037730626
2.010344209437878312641529731720143
1.999795933452669925783583536567984
2.000000072231393330599606713662298
2.000000000000015319238844912588532
2.000000000000000000000000018934481
2.000000000000000000000000000000000

Code:

secant_3(Decimal(5), Decimal(4), 1e-32)

5
4
3.081967213114754098360655737704918
2.286218829717811307322668037730626
2.034337291023909027923796138229571
2.000576313421516741692818211990178
2.000000166004797850203846960058336
2.000000000000013778794929746133383
2.000000000000000000000000000094928
2.000000000000000000000000000000000

These are the corresponding programs for the HP-42S:

Code:

00 { 76-Byte Prgm }

01▸LBL "SEC-1"

02 STO 02 # ε

03 R↓

04 STO 01 # x_5

05 R↓

06 STO 00 # x_4

07 XEQ F

08 STO 03 # f_4

09 RCL 01 # x_5

10 XEQ F

11 STO 04 # f_5

12 RCL 03

13 RCL- 04

14 RCL 00

15 RCL- 01

16 ÷

17 STO 05 # f_45 = (f_4 - f_5) ÷ (x_4 - x_5)

18▸LBL 00

19 RCL 02

20 RCL 04

21 ABS

22 X<Y?   # abs(f_5) < ε

23 GTO 01

24 RCL 01

25 RCL 04

26 RCL 05

27 ÷

28 -

29 STO 06 # x_6 = x_5 - f_5 ÷ f_45

30 XEQ F

31 X<> 04 # f_6 <> f_5

32 RCL- 04

33 RCL 06

34 X<> 01 # x_6 <> x_5

35 RCL- 01

36 ÷

37 STO 05 # f_45 = f_56 = (f_5 - f_6) ÷ (x_5 - x_6)

38 GTO 00

39▸LBL 01

40 RCL 01 # x_5

41 RTN

42▸LBL F

43 STOP

44 3

45 Y↑X

46 8

47 -

48 END

Code:

00 { 127-Byte Prgm }

01▸LBL "SEC-2"

02 STO 02 # ε

03 R↓

04 STO 01 # x_5

05 R↓

06 STO 00 # x_4

07 XEQ F

08 STO 03 # f_4

09 RCL 01 # x_5

10 XEQ F

11 STO 04 # f_5

12 RCL 03

13 RCL- 04

14 RCL 00

15 RCL- 01

16 ÷

17 STO 05 # f_45 = (f_4 - f_5) ÷ (x_4 - x_5)

18 RCL 01

19 RCL 04

20 RCL 05

21 ÷

22 -

23 STO 06 # x_6 = x_5 - f_5 ÷ f_45

24 XEQ F

25 STO 07 # f_6

26 RCL 04

27 RCL- 07

28 RCL 01

29 RCL- 06

30 ÷

31 STO 08 # f_56 = (f_5 - f_6) ÷ (x_5 - x_6)

32 RCL 05

33 RCL- 08

34 RCL 00

35 RCL- 06

36 ÷

37 STO 09 # f_456 = (f_45 - f_56) ÷ (x_4 - x_6)

38▸LBL 00

39 RCL 02

40 RCL 07

41 ABS

42 X<Y?   # abs(f_6) < ε

43 GTO 01

44 RCL 06

45 RCL 07

46 RCL 06

47 RCL- 01

48 RCL× 09

49 RCL+ 08

50 ÷

51 -

52 STO 10 # x_7 = x_6 - f_6 ÷ ((x_6 - x_5) × f_456 + f_56)

53 XEQ F

54 X<> 07 # f_7 <> f_6

55 RCL- 07

56 RCL 06

57 RCL- 10

58 ÷

59 X<> 08 # f_56 <> f_67 = (f_6 - f_7) ÷ (x_6 - x_7)

60 RCL- 08

61 RCL 01

62 RCL- 10

63 ÷

64 STO 09 # f_456 = f_567 = (f_56 - f_67) ÷ (x_5 - x_7)

65 RCL 10

66 X<> 06

67 STO 01 # x_5, x_6 = x_6, x_7

68 GTO 00

69▸LBL 01

70 RCL 06 # x_6

71 RTN

72▸LBL F

73 STOP

74 3

75 Y↑X

76 8

77 -

78 END

Code:

00 { 195-Byte Prgm }

01▸LBL "SEC-3"

02 STO 02 # ε

03 R↓

04 STO 01 # x_5

05 R↓

06 STO 00 # x_4

07 XEQ F

08 STO 03 # f_4

09 RCL 01 # x_5

10 XEQ F

11 STO 04 # f_5

12 RCL 03

13 RCL- 04

14 RCL 00

15 RCL- 01

16 ÷

17 STO 05 # f_45 = (f_4 - f_5) ÷ (x_4 - x_5)

18 RCL 01

19 RCL 04

20 RCL 05

21 ÷

22 -

23 STO 06 # x_6 = x_5 - f_5 ÷ f_45

24 XEQ F

25 STO 07 # f_6

26 RCL 04

27 RCL- 07

28 RCL 01

29 RCL- 06

30 ÷

31 STO 08 # f_56 = (f_5 - f_6) ÷ (x_5 - x_6)

32 RCL 05

33 RCL- 08

34 RCL 00

35 RCL- 06

36 ÷

37 STO 09 # f_456 = (f_45 - f_56) ÷ (x_4 - x_6)

38 RCL 06

39 RCL 07

40 RCL 06

41 RCL- 01

42 RCL× 09

43 RCL+ 08

44 ÷

45 -

46 STO 10 # x_7 = x_6 - f_6 ÷ ((x_6 - x_5) × f_456 + f_56)

47 XEQ F

48 STO 11 # f_7

49 RCL 07

50 RCL- 11

51 RCL 06

52 RCL- 10

53 ÷

54 STO 12 # f_67 = (f_6 - f_7) ÷ (x_6 - x_7)

55 RCL 08

56 RCL- 12

57 RCL 01

58 RCL- 10

59 ÷

60 STO 13 # f_567 = (f_56 - f_67) ÷ (x_5 - x_7)

61 RCL 09

62 RCL- 13

63 RCL 00

64 RCL- 10

65 ÷

66 STO 14 # f_4567 = (f_456 - f_567) ÷ (x_4 - x_7)

67▸LBL 00

68 RCL 02

69 RCL 11

70 X<Y?   # abs(f_7) < ε

71 GTO 01

72 RCL 10

73 RCL 11

74 RCL 10

75 RCL- 01

76 RCL× 14

77 RCL+ 13

78 RCL 10

79 RCL- 06

80 ×

81 RCL+ 12

82 ÷

83 -

84 STO 15 # x_8 = x_7 - f_7 ÷ (((x_7 - x_5) × f_4567 + f_567) × (x_7 - x_6) + f_67)

85 XEQ F

86 X<> 11 # f_8 <> f_7

87 RCL- 11

88 RCL 10

89 RCL- 15

90 ÷

91 X<> 12 # f_67 <> f_78 = (f_7 - f_8) ÷ (x_7 - x_8)

92 RCL- 12

93 RCL 06

94 RCL- 15

95 ÷

96 X<> 13 # f_567 <> f_678 = (f_67 - f_78) ÷ (x_6 - x_8)

97 RCL- 13

98 RCL 01

99 RCL- 15

100 ÷

101 STO 14 # f_4567 = f_5678 = (f_567 - f_678) ÷ (x_5 - x_8)

102 RCL 15

103 X<> 10

104 X<> 06

105 STO 01 # x_5, x_6, x_7 = x_6, x_7, x_8

106 GTO 00

107▸LBL 01

108 RCL 10 # x_7

109 RTN

110▸LBL F

111 STOP

112 3

113 Y↑X

114 8

115 -

116 END

From the paper: \(s_1 \doteq 1.618, \; s_2 \doteq 1.839, \; s_3 \doteq 1.928\)
It appears that the benefit of order 3 compared to order 2 is not that big.

Registers

Interestingly the 3 variants share registers and also some of the code.

SEC-1:
00: x_4
01: x_5
02: ε
03: f_4
04: f_5
05: f_45 = (f_4 - f_5) ÷ (x_4 - x_5)
06: x_6 = x_5 - f_5 ÷ f_45

SEC-2:
00: x_4
01: x_5
02: ε
03: f_4
04: f_5
05: f_45 = (f_4 - f_5) ÷ (x_4 - x_5)
06: x_6 = x_5 - f_5 ÷ f_45
07: f_6
08: f_56 = (f_5 - f_6) ÷ (x_5 - x_6)
09: f_456 = (f_45 - f_56) ÷ (x_4 - x_6)
10: x_7 = x_6 - f_6 ÷ ((x_6 - x_5) × f_456 + f_56)

SEC-3:
00: x_4
01: x_5
02: ε
03: f_4
04: f_5
05: f_45 = (f_4 - f_5) ÷ (x_4 - x_5)
06: x_6 = x_5 - f_5 ÷ f_45
07: f_6
08: f_56 = (f_5 - f_6) ÷ (x_5 - x_6)
09: f_456 = (f_45 - f_56) ÷ (x_4 - x_6)
10: x_7 = x_6 - f_6 ÷ ((x_6 - x_5) × f_456 + f_56)
11: f_7
12: f_67 = (f_6 - f_7) ÷ (x_6 - x_7)
13: f_567 = (f_56 - f_67) ÷ (x_5 - x_7)
14: f_4567 = (f_456 - f_567) ÷ (x_4 - x_7)
15: x_8 = x_7 - f_7 ÷ (((x_7 - x_5) × f_4567 + f_567) × (x_7 - x_6) + f_67)
16: f_8

Example

This is the example from the paper in table 2.
Keep hitting the R/S key to get the next result.

5 ENTER
4 ENTER
1E-32
XEQ "SEC-2"

5
4
3.081967213114754098360655737704918
2.286218829717811307322668037730626
2.010344209437878312641529731720143
1.999795933452669925783583536567984
2.000000072231393330599606713662298
2.000000000000015319238844912588532
2.000000000000000000000000018934481
2.000000000000000000000000000000000