# Zeros and Roots function from "Numerical Computing with Matlab" textbook

10 views (last 30 days)
Ryan M on 6 Oct 2017
Answered: Patrick Hew on 25 Jun 2020
Hi, I'm trying to answer question 4.7 from the Clever Moler MATLAB textbook, which can be found here (chapter 4): https://au.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/zeros.pdf
Does anyone know why this is? Also out of curiosity, does anyone know what sort of stopping criteria the function fzerotx uses? Note that the actual zero between the interval 0 and pi is 2.4048.
Here is the code for the function fzerotx. More information on the code can be found in section 4.7 of the link.
function b = fzerotx(F,ab,varargin)
%FZEROTX Textbook version of FZERO.
% x = fzerotx(F,[a,b]) tries to find a zero of F(x) between a and b.
% F(a) and F(b) must have opposite signs. fzerotx returns one
% end point of a small subinterval of [a,b] where F changes sign.
% Arguments beyond the first two, fzerotx(F,[a,b],p1,p2,...),
% are passed on, F(x,p1,p2,..).
%
% Examples:
% fzerotx(@sin,[1,4])
% F = @(x) sin(x); fzerotx(F,[1,4])
% Copyright 2014 The MathWorks, Inc.
% Initialize.
a = ab(1);
b = ab(2);
fa = F(a,varargin{:});
fb = F(b,varargin{:});
if sign(fa) == sign(fb)
error('Function must change sign on the interval')
end
c = a;
fc = fa;
d = b - c;
e = d;
% Main loop, exit from middle of the loop
while fb ~= 0
% The three current points, a, b, and c, satisfy:
% f(x) changes sign between a and b.
% abs(f(b)) <= abs(f(a)).
% c = previous b, so c might = a.
% The next point is chosen from
% Bisection point, (a+b)/2.
% Secant point determined by b and c.
% Inverse quadratic interpolation point determined
% by a, b, and c if they are distinct.
if sign(fa) == sign(fb)
a = c; fa = fc;
d = b - c; e = d;
end
if abs(fa) < abs(fb)
c = b; b = a; a = c;
fc = fb; fb = fa; fa = fc;
end
% Convergence test and possible exit
m = 0.5*(a - b);
tol = 2.0*eps*max(abs(b),1.0);
if (abs(m) <= tol) || (fb == 0.0)
break
end
% Choose bisection or interpolation
if (abs(e) < tol) || (abs(fc) <= abs(fb))
% Bisection
d = m;
e = m;
else
% Interpolation
s = fb/fc;
if (a == c)
% Linear interpolation (secant)
p = 2.0*m*s;
q = 1.0 - s;
else
q = fc/fa;
r = fb/fa;
p = s*(2.0*m*q*(q - r) - (b - c)*(r - 1.0));
q = (q - 1.0)*(r - 1.0)*(s - 1.0);
end;
if p > 0
q = -q;
else
p = -p;
end;
% Is interpolated point acceptable
if (2.0*p < 3.0*m*q - abs(tol*q)) && (p < abs(0.5*e*q))
e = d;
d = p/q;
else
d = m;
e = m;
end;
end
% Next point
c = b;
fc = fb;
if abs(d) > tol
b = b + d;
else
b = b - sign(b-a)*tol;
end
fb = F(b,varargin{:});
end

Patrick Hew on 25 Jun 2020
The 'trick' is that the call
z = fzerotx(@besselj,[0 pi],0)
doesn't invoke besselj in the way that you'd expect. What happens is that the line
fb = F(b,varargin{:});
is executed as
fb = besselj(pi,0);
which is equivalent to
fb = 0;
Consequently fzerotx falls through its main loop (while fb ~= 0) and returns b = pi.
A correct call is
z = fzerotx(@(z) besselj(0,z), [0 pi])
which returns z = 2.4048 as expected.