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 ##
24 if (columns(x0) > 1)
25 root = NA;
26 return;
27 end
29 iterations = 0;
30 xn = x0;
31 f_val = f(x0);
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
42 root = xn;
43 end