From d7a73d33b0567c9abd72151b2972d7c0eb66e6a1 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Mon, 15 Oct 2012 21:22:54 -0400
Subject: [PATCH 1/1] Add Newton's method and some tests.

---
 newtons_method.m | 43 +++++++++++++++++++++++++++++++++++++++++++
 run-tests.m      | 30 ++++++++++++++++++++++++++++++
 2 files changed, 73 insertions(+)
 create mode 100644 newtons_method.m

diff --git a/newtons_method.m b/newtons_method.m
new file mode 100644
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 aa285c3..2892ba0 100755
--- 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);
-- 
2.43.2