Add Newton's method and some tests.

diff --git a/newtons_method.m b/newtons_method.m
new file mode 100644 (file)
index 0000000..8bdd090
--- /dev/null
+++ b/newtons_method.m
@@ -0,0 +1,43 @@
+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

diff --git a/run-tests.m b/run-tests.m
index aa285c32e38fb3f70747a75402e57ef6a8515ed2..2892ba0e84f916571960177c62bb297068d7b448 100755 (executable)
--- a/run-tests.m
+++ b/run-tests.m
@@ -72,3 +72,33 @@ expected_fp = [1.0729, 1.0821];
 unit_test_equals("Homework #3 problem #3i fixed point is correct", ...
                 expected_fp, ...
                 fixed_point_method(my_g, tol, u0));
+
+
+f = @(x) x^6 - x - 1;
+f_prime = @(x) 6*x^5 - 1;
+tol = 1/1000000;
+x0 = 2;
+expected_root = 1.1347;
+unit_test_equals("Newton's method agrees with Haskell", ...
+                expected_root, ...
+                newtons_method(f, f_prime, tol, x0));
+
+
+
+f1 = @(u) u(1)^2 + u(1)*u(2)^3 - 9;
+f2 = @(u) 3*u(1)^2*u(2) - u(2)^3 - 4;
+f = @(u) [f1(u); f2(u)];
+## The partials for the Jacobian.
+f1x = @(u) 2*u(1) + u(2)^3;
+f1y = @(u) 3*u(1)*u(2)^2;
+f2x = @(u) 6*u(1)*u(2);
+f2y = @(u) 3*u(1)^2 - 3*u(2)^2;
+## f_prime == Jacobian.
+f_prime = @(u) [ f1x(u), f1y(u); f2x(u), f2y(u) ];
+tol = 1 / 10^12;
+u0 = [1.2; 2.5];
+expected_root = [1.33635; 1.75424];
+[actual_root, iterations] = newtons_method(f, f_prime, tol, u0);
+unit_test_equals("Homework #3 problem #4 root is correct", ...
+                expected_root, ...
+                actual_root);