MATLAB - Secant method produces NaN

Question

I'm writing a secant method in MATLAB, which I want to iterate through exactly n times.

function y = secantmeth(f,xn_2,xn_1,n)

xn = (xn_2*f(xn_1) - xn_1*f(xn_2))/(f(xn_1) - f(xn_2));
k = 0;
while (k < n)
    k = k + 1;
    xn_2 = xn_1;
    xn_1 = xn;
    xn = (xn_2*f(xn_1) - xn_1*f(xn_2))/(f(xn_1) - f(xn_2));
end
y = xn;
end

I believe the method works for small values of n, but even something like n = 9 produces NaN. My guess is that the quantity f(xn_1) - f(xn_2) is approximately zero, which causes this error. How can I prevent this?

Examples:

Input 1

eqn = @(x)(x^2 + x -9)

secantmeth(eqn,2,3,5)

Input 2

eqn = @(x)(x^2 + x - 9)

secantmeth(eqn, 2, 3, 9)

Output 1

2.7321

Output 2

NaN

Can you give us an example input you use to this function? I mean `f, xn_2, xn_1, n` — Karls, Mar 20 '19 at 14:04
In Newtons method this problem also exists. The common way of solving it is before computing the update on `xn`, you check the divisor. If its close to zero, you use a different method, such as bisection, and keep going. — Ander Biguri, Mar 20 '19 at 14:06
Use the formula with the inverse slope, `xn = xn_1 - fn_1*(xn_1-xn_2)/(fn_1-fn_2)`, somewhere it is proven that it is more numerically stable than the symmetric midpoint formula. — Lutz Lehmann, Mar 20 '19 at 14:44

gnovice · Accepted Answer · 2019-03-20T14:51:09.090

The value for xn will be NaN when xn_2 and xn_1 are exactly equal, which results in a 0/0 condition. You need to have an additional check in your while loop condition to see if xn_1 and x_n are equal (or, better yet, within some small tolerance of one another), thus suggesting that the loop has converged on a solution and can't iterate any further:

...
while (k < n) && (xn_1 ~= xn)
  k = k + 1;
  xn_2 = xn_1;
  xn_1 = xn;
  xn = (xn_2*f(xn_1) - xn_1*f(xn_2))/(f(xn_1) - f(xn_2));
end
...

As Ander mentions in a comment, you could then continue with a different method after your while loop if you want to try and get a more accurate approximation:

...
if (xn_1 == xn)  % Previous loop couldn't iterate any further
  % Try some new method
end
...

And again, I would suggest reading through this question to understand some of the pitfalls of floating-point comparison (i.e. == and ~= aren't usually the best operators to use for floating-point numbers).

MATLAB - Secant method produces NaN

1 Answers1