Root-Bracketing Quartile Algorithm

by

Namir C. Shammas

Numerical analysis offers a rich set of algorithms to solve for real roots of single nonlinear functions. The simplest class of root-seeking algorithms is the root-bracketing algorithms. These algorithms work by taking an interval that contains a root and zooming in on that root. The simplest, and slowest, root-seeking algorithm is the Bisection algorithm. This method works by selecting the middle of the interval and determining which half of the interval to discard. Each iteration reduces the root-hosting interval by half until the interval becomes small enough to contain the root within an accepted tolerance value.

The pseudo code for the Bisection method is:

Given: An interval [A, B] that contains a root for function f(X), such that f(A)*f(B) < 0.

1. Select the middle point M as (A + B)/2
2. If f(M) * f(A) > 0 then set A = M else set B = M
3. If |A - B| > Tolerance then resume at step 1
4. Root is M

The Bisection method is the slowest root-seeking method. However, since the method uses the sign of the targeted function, convergence is guaranteed if the initial conditions are true.

Another root-seeking algorithm is the method of False Position. It is more efficient than the Bisection method. The False Position algorithm looks at the function values at A and B as weights used to calculate the next guess. The following equation shows how to calculate the new guess:

M = (f(B) * A - f(A) * B) / (f(B) - f(A))

The algorithm performs a function sign test, just like with the Bisection algorithm, to determine which interval end is replaced by the new guess M. However, this method does not shrink the interval as one might expect. To stop iteration, an implementation of the method of False Position must compare the most two recent values M. Iteration resumes as long as the absolute difference of these values exceeds a tolerance value.

The new Quartile algorithm adopts a middle strategy regarding the use of the function values at A and B. The method compares the absolute values of f(A) and f(B) to determine how close to A or B the next guess for the root will be located.

Here is the pseudo code for Quartile algorithm:

Given: An interval [A, B] that contains a root for function f(X), such that f(A)*f(B) < 0.

1. If |f(A)| < |f(B)| then
   1. M = A + alpha * (B - A)
   2. Else
   3. M = A + (1 - alpha) * (B - A)
2. If f(M) * f(A) > 0 then set A = M else set B = M
3. If |A - B| > Tolerance then resume at step 1
4. Root is M

The coefficient alpha can be in the range of 0.25 to 0.30. When the value for alpha is 0.5, the Quartile algorithm behaves just like the Bisection algorithm. Thus, conceptually, the Bisection method is a special case of the Quartile algorithm.

To compare the Bisection, Quartile, and False Position algorithms, consider the function f(x) = e^x – 3x² = 0 and tolerance = 1E-7. This function has roots near -0.459, 0.910, and 3.733.  The following table shows the initial values for A and B, and the number of function calls for the Bisection, Quartile, and False Position algorithms. Comparing the Bisection and Quartile methods, the table shows that the Quartile algorithm requires fewer function calls to obtain a root in all the cases listed. Comparing the Quartile and False Position methods, the table shows the Quartile algorithm is no worse than the False Position in 6 (out of 10) cases.

Option Explicit

Function F(X As Double) As Double
  F = Exp(X) - 3 * X ^ 2
End Function

Sub test()
  Const TOLER = 0.0000001
  Const ALPHA = 0.27
  Dim A As Double, B As Double, M As Double, LastM As Double
  Dim FA As Double, FB As Double, FM As Double
  Dim ColA As Integer, ColB As Integer, ColM As Integer
  Dim Row As Integer, Count As Integer
  
  Range("C2:Z1000").Clear
  
  Row = 2
  Do Until Cells(Row, 1) = ""
  
    ' Perform the Bisection Method
    ColA = 4
    ColB = 5
    ColM = 6
    
    A = Cells(Row, 1)
    B = Cells(Row, 2)
    
    FA = F(A)
    FB = F(B)
    Count = 2
    Do
      M = (A + B) / 2
      FM = F(M)
      Count = Count + 1
      If FM * FA > 0 Then
        A = M
        FA = FM
      Else
        B = M
        FB = FM
      End If
    Loop Until Abs(A - B) < TOLER
  
    Cells(Row, ColA - 1) = Count
    Cells(Row, ColA) = A
    Cells(Row, ColB) = B
    Cells(Row, ColM) = M
    Cells(Row, ColM + 1) = FM
    
    ' Perform the Quartile Method
    ColA = ColM + 3
    ColB = ColA + 1
    ColM = ColB + 1
    
    A = Cells(Row, 1)
    B = Cells(Row, 2)
    
    FA = F(A)
    FB = F(B)
    Count = 2
    Do
      If Abs(FA) < Abs(FB) Then
        M = A + ALPHA * (B - A)
      Else
        M = A + (1 - ALPHA) * (B - A)
      End If
      FM = F(M)
      Count = Count + 1
      If FM * FA > 0 Then
        A = M
        FA = FM
      Else
        B = M
        FB = FM
      End If
    Loop Until Abs(A - B) < TOLER
   
    Cells(Row, ColA - 1) = Count
    Cells(Row, ColA) = A
    Cells(Row, ColB) = B
    Cells(Row, ColM) = M
    Cells(Row, ColM + 1) = FM
    
    ' Perform the False Position Method
    ColA = ColM + 3
    ColB = ColA + 1
    ColM = ColB + 1
      
    A = Cells(Row, 1)
    B = Cells(Row, 2)
    
    FA = F(A)
    FB = F(B)
    M = A
    Count = 2
    Do
      LastM = M
      M = (FB * A - FA * B) / (FB - FA)
      FM = F(M)
      Count = Count + 1
      If FM * FA > 0 Then
        A = M
        FA = FM
      Else
        B = M
        FB = FM
      End If
    Loop Until (Abs(M - LastM) < TOLER) Or (FA = FB)
   
    Cells(Row, ColA - 1) = Count
    Cells(Row, ColA) = A
    Cells(Row, ColB) = B
    Cells(Row, ColM) = M
    Cells(Row, ColM + 1) = FM
   
    Row = Row + 1
  Loop
  
End Sub