10 views (last 30 days)

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 Cleve Moler

% 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

% Inverse quadratic interpolation

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.

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.