function [x, fx, a, b] = newtroot(fname, a, b, xtol, ftol, fig, start) % Newton's method for a differentiable function F mapping the real line to % the real line, with bisection safeguards % input % fname : string giving name of function that evaluates F(x) and % its derivative F'(x), with call: [f, deriv] = fname(x) % a, b: points satisfying a= b error('newtroot: a >= b') end [fa, deriva] = feval(fname, a); [fb, derivb] = feval(fname, b); if fa*fb > 0 error('newtroot: F(a), F(b) have same sign') end if start < 0 x = a; fx = fa; deriv = deriva; else x = b; fx = fb; deriv = derivb; end if fig > 0 figure(fig) hold on plot(a, fa, 'r+') plot(b, fb, 'r+') end maxit = 20; iter = 0; while b-a > xtol & abs(fx) > ftol & iter < maxit iter = iter + 1; newtonstep = -fx/deriv; % OK to divide by zero % make sure newton step doesn't go the wrong way or go past % the bisection step if (x == a & newtonstep < 0) | (x == b & newtonstep > 0) |... abs(newtonstep) > (b-a)/2 x = (a + b)/2; fprintf('bisection step to %g', x) else x = x + newtonstep; fprintf('Newton step to %20.15e', x); end [fx, deriv] = feval(fname, x); if fx*fa > 0 % update a a = x; fa = fx; else % update b b = x; fb = fx; end fprintf(', new values a= %g fa = %g b = %g fb = %g\n', a, fa, b, fb) if fig > 0 plot(x, fx, 'm*') end end