3
$\begingroup$

I have coded a simple program for a root finding problem using Halley's method. Here is the code:

function [x,n] = root_finder_Halley(f0,f1,f2,x)

i = true;
n = 0;
while i
n = n + 1;

x(n+1) = x(n) - ( 2*f0(x(n))*f1(x(n)) )/...
    ( 2*f1(x(n))^2 - f0(x(n))*f2(x(n)) );

dx(n) = abs(x(n+1) - x(n));
i = dx(n) > eps(x(n+1));
end

where eps(x) is the machine epsilon at the value x. I need to carefully make a statement for the loop's convergence. The goal is to achieve the best precision allowed by the computer, so I did not pre-define a desired tolerance value, instead, I've set a dynamic tolerance equal to eps(x(n+1)). This means that the loop should continue until the difference of result at two consecutive loops is no more than their epsilon. This program works very good, but sometimes it goes forever. an example of this program fail is given below:

x =
         0.049278137931035
        0.0821484600250533
        0.0853747083523107
        0.0853759293572672
        0.0853759293572688
        0.0853759293572672
        0.0853759293572688
dx =
        0.0328703220940183
       0.00322624832725742
      1.22100495650312e-06
      1.60982338570648e-15
      1.62370117351429e-15
      1.60982338570648e-15
      1.62370117351429e-15
eps(x) =
      6.93889390390723e-18
      1.38777878078145e-17
      1.38777878078145e-17
      1.38777878078145e-17
      1.38777878078145e-17
      1.38777878078145e-17
      1.38777878078145e-17

As can be seen, dx oscillates between two values greater than eps(x). I even tried the difference of difference (dx(n+1)-dx(n)) but it was oscillating too. Is there a better convergence criterion (other than predefined tolerance) for this issue?

$\endgroup$
3
  • 1
    $\begingroup$ This is normal. The functions f0, f1, f2 have evaluation errors that combine with those of the formula, so getting a loop in the last 2 or 3 digits is common. If the step size reduction stalls out, not even reducing by a factor of 2, you could switch to a simpler method with less computational effort like the secant method. $\endgroup$ Commented Feb 9, 2023 at 15:46
  • $\begingroup$ Thank you dr. Lehmann, I guess I should bring the program in an if-statement and define that if the number of iterations exceed a specified number, then switch to secant. However I want to try my best to somehow "detect" stall of step size reduction using eps. $\endgroup$ Commented Feb 9, 2023 at 17:49
  • $\begingroup$ @LutzLehmann Please don't use comments to answer the question. $\endgroup$ Commented Feb 9, 2023 at 19:57

1 Answer 1

3
$\begingroup$

You are using the wrong measure for the achievable accuracy. $|x|\mu$ is in some way a lower bound for the error, the smallest increment that would give a different value, so that any smaller increments leave the value of $x$ unchanged.

However, what you actually get in the numerical picture is a thick line of floating-point values if you zoom in sufficiently. This thick line will give an interval of $x$ values from where the lower bound of the line is above zero to where the upper bound is below zero. The radius of this interval can be computed from an estimate of the thickness of the line and the slope of it, the derivative of the function. So what is needed here is the interval of values that the numerical evaluation of $f(x)$ can take even if the exact evaluation has the value zero.

To get a function value near zero for arguments close to a root requires, except in some trivial cases, that unrelated terms in the function expression cancel. This quite naturally leads to catastrophic cancellation in the most extreme way, where the error level of the terms is preserved even if the exact value of the composition is zero.

A rough estimate of this error can be $|f|(x)\mu$ where $|f|$ has the same computational algorithm as $f$, only that all sums and differences are replaced by sums of absolute values. For a polynomial $f(x)=c_0+c_1x+...+c_nx^n$ this would give $|f|(x)=|c_0|+|c_1||x|+...+|c_n||x|^n$. This works reasonably well, see https://math.stackexchange.com/questions/3021532/is-it-possible-for-the-bisection-method-to-provide-fake-zeros/3021704#3021704. The actual error propagation will look a little more complex.

So the error bound should be changed to eps(x)=mu*abs(fabs(x)/f1(x)), with fabs as described above or some better estimate for the scale factor in the error of the floating-point evaluation of $f$. This will give the radius of an interval where $f$ changes sign, and where the numerical values left and right from the interval boundaries will most likely strictly above or below zero. So if the length of the Newton step falls below this radius, it is most likely that the iteration point is inside such an interval around a root or close to it.

One could guess that the exact $f$ is almost linear at the scale eps(x), so computing the secant root for the points x-10*eps(x), x+10*eps(x) is probably as good as it gets in getting close to the root. Or perform linear regression for an arithmetic sequence at this scale centered at $x$, taking the root of the regression line as last root approximation.

$\endgroup$
4
  • $\begingroup$ As a sidenote, in my experience, dichotomy is the best algorithm in terms of achievable precision on the result, since it only evaluates $f$ without further operations. However, efficiency is not on its side... Could be used to further refine the solution though. $\endgroup$
    – Laurent90
    Commented Feb 10, 2023 at 19:35
  • 2
    $\begingroup$ @Laurent90 Do you mean the bisection method? Then yes, it will reduce the interval to zero, but the result will be one of the many sign changes in this small interval that the numerical evaluation of $f$ provides, see the linked answer. $\endgroup$ Commented Feb 10, 2023 at 20:11
  • $\begingroup$ @LutzLehmann Let $x(n+1) = x(n) + H(x(n))$ be the recurrence relation of any algorithm where $H(x(n))$ is the iteration function. According to your approximate error measurement, can we generalize the tolerance to be $\mu |H|(x(n))$ ? $\endgroup$ Commented Feb 16, 2023 at 6:50
  • 1
    $\begingroup$ Yes. If the difference of the step falls below that point, the further iteration is mostly random noise. If you can decompose $H$ into a quotient of products (denominator may be trivial), and you know which factor in the numerator has the zero, then the other factors are unlikely to have catastrophic cancellations, their relative errors do not contribute to the error bound in a major war. So only the zero factor contributes to the error level, needs to be analyzed for cancellation, the others can be evaluated as they are. $\endgroup$ Commented Feb 16, 2023 at 9:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.