Post Reply 
Variant of the Secant Method
12-09-2020, 04:33 AM (This post was last modified: 12-09-2020 09:13 PM by ttw.)
Post: #1
Variant of the Secant Method
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
Find all posts by this user
Quote this message in a reply
12-09-2020, 05:40 AM
Post: #2
RE: Variant of the Second Method
(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.

  
All My Articles & other Materials here:  Valentin Albillo's HP Collection
 
Visit this user's website Find all posts by this user
Quote this message in a reply
12-09-2020, 04:38 PM (This post was last modified: 12-09-2020 11:21 PM by Albert Chan.)
Post: #3
RE: Variant of the Second Method
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
Find all posts by this user
Quote this message in a reply
12-09-2020, 11:51 PM
Post: #4
RE: Variant of the Secant Method
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
Find all posts by this user
Quote this message in a reply
12-10-2020, 01:54 PM (This post was last modified: 12-10-2020 02:05 PM by Namir.)
Post: #5
RE: Variant of the Secant Method
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
Find all posts by this user
Quote this message in a reply
12-10-2020, 02:56 PM (This post was last modified: 12-11-2020 01:18 PM by Namir.)
Post: #6
RE: Variant of the Secant Method
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.
Find all posts by this user
Quote this message in a reply
12-11-2020, 04:34 PM (This post was last modified: 02-07-2021 11:57 AM by Albert Chan.)
Post: #7
RE: Variant of the Secant Method
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
Find all posts by this user
Quote this message in a reply
12-11-2020, 04:46 PM (This post was last modified: 02-07-2021 11:59 AM by Albert Chan.)
Post: #8
RE: Variant of the Secant Method
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
Find all posts by this user
Quote this message in a reply
12-12-2020, 04:22 AM
Post: #9
RE: Variant of the Secant Method
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
Find all posts by this user
Quote this message in a reply
12-12-2020, 02:41 PM
Post: #10
RE: Variant of the Secant Method
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.
Find all posts by this user
Quote this message in a reply
12-12-2020, 04:27 PM
Post: #11
RE: Variant of the Secant Method
(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')
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)