From: Michael Orlitzky Date: Tue, 16 Oct 2012 01:22:54 +0000 (-0400) Subject: Add Newton's method and some tests. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=d7a73d33b0567c9abd72151b2972d7c0eb66e6a1;p=octave.git Add Newton's method and some tests. --- 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);