function [root, iterations] = newtons_method(f, f_prime, epsilon, x0) ## Find a root of the function `f` with initial guess x0. ## ## INPUTS: ## ## * ``f`` - The function whose root we seek. Must return a column ## vector or scalar. ## ## * ``f_prime`` - The derivative or Jacobian of ``f``. ## ## * ``epsilon`` - We stop when the value of `f` becomes less ## than epsilon. ## ## * ``x0`` - An initial guess. Either 1d or a column vector. ## ## OUTPUTS: ## ## * ``root`` - The root that we found. ## ## * ``iterations`` - The number of iterations that we performed ## during the search. ## if (columns(x0) > 1) root = NA; return; end iterations = 0; xn = x0; f_val = f(x0); while (norm(f_val, Inf) > epsilon) ## This uses the modified algorithm (2.11.5) in Atkinson. ## Should work in 1d, too. delta = f_prime(xn) \ f_val; xn = xn - delta; f_val = f(xn); iterations = iterations + 1; end root = xn; end