newtons_method.m

   1 function [root, iterations] = newtons_method(f, f_prime, epsilon, x0)
   2   % Find a root of the function `f` with initial guess x0.
   3   %
   4   % INPUTS:
   5   %
   6   %   * ``f`` - The function whose root we seek. Must return a column
   7   %     vector or scalar.
   8   %
   9   %   * ``f_prime`` - The derivative or Jacobian of ``f``.
  10   %
  11   %   * ``epsilon`` - We stop when the value of `f` becomes less
  12   %     than epsilon.
  13   %
  14   %   * ``x0`` - An initial guess. Either 1d or a column vector.
  15   %
  16   % OUTPUTS:
  17   %
  18   %   * ``root`` - The root that we found.
  19   %
  20   %   * ``iterations`` - The number of iterations that we performed
  21   %     during the search.
  22   %
  23
  24   if (columns(x0) > 1)
  25     root = NA;
  26     return;
  27   end
  28
  29   iterations = 0;
  30   xn = x0;
  31   f_val = f(x0);
  32
  33   while (norm(f_val, Inf) > epsilon)
  34     % This uses the modified algorithm (2.11.5) in Atkinson.
  35     % Should work in 1d, too.
  36     delta = f_prime(xn) \ f_val;
  37     xn = xn - delta;
  38     f_val = f(xn);
  39     iterations = iterations + 1;
  40   end
  41
  42   root = xn;
  43 end