examples/algorithms/Newton.pm

   1 use v6-alpha;
   2
   3 package Newton;
   4
   5 sub update_guess(Num $guess is rw, Num $target, Code $f, Code $fprime) {
   6     $guess += ($target - $f($guess)) / $fprime($guess);
   7 }
   8
   9 sub approxfprime(Code $f, Num $x) {
  10     my Num $delta = 0.1;
  11     return ($f($x + $delta) - $f($x - $delta))/(2 * $delta);
  12 }
  13
  14 sub newton(
  15     Num  $target,
  16     Code $f,
  17     Code :$fprime         = &approxfprime.assuming( f => $f ),
  18     Num  :$epsilon        = 0.0005,
  19     Num  :$max_iterations = 500,
  20     Bool :$verbose        = 0
  21 ) returns Num is export {
  22     my Num $guess  = $target / 2;
  23     my Int $count  = 1;
  24
  25     while (abs($f($guess) - $target) > $epsilon) {
  26
  27         if ($count > $max_iterations) {
  28             die "excessive looping";
  29         }
  30
  31         update_guess($guess, $target, $f, $fprime);
  32
  33         say "{$count++}: $guess" if $verbose;
  34     }
  35     return $guess;
  36 }
  37
  38 =head1 NAME
  39
  40 Newton - performs one dimensional Newton's method
  41
  42 =head1 SYNOPSIS
  43
  44     use v6-alpha;
  45
  46     use Newton;
  47
  48     sub f(Num $x) { return $x ** 3; }
  49
  50     say "The cube root of 7 is: {newton(7, &f, $verbose => 1)}";
  51
  52 =head1 DESCRIPTION
  53
  54 Newton's method is a simple method for finding roots, like square roots or
  55 cube roots (as in the example above).  Read a Calculus textbook for details.
  56 This implementation is meant to show the power of Perl 6 defaults and
  57 Currying in a real context.
  58
  59 The newton function must receive at least:
  60
  61     y - target value for y
  62     f - function of one variable
  63
  64 It looks for an x such that f(x) is y.
  65
  66 You may optionally supply:
  67
  68     epsilon - tolerance (how accurate the answer must be)
  69     fprime  - an exact first derivative of f (which Newton's method
  70               uses to update its guess at each iteration)
  71     verbose - a boolean which will turn on per iteration printing
  72
  73 If you omit the fprime, a second order centered difference is provided
  74 as a curried default.
  75
  76 All of the optional parameters should be supplied by name.
  77
  78 Note that the above code is not robust.  (If you need a real implementation
  79 consult a book, like Numerical Recipies.)  In particular, it suffers from
  80 convergence problems for some functions, which all Newton methods do.  All
  81 it does in such cases is die with a complaint.  Real solvers usually do
  82 a bit of algorithm switching.  When their fast method fails, they fall back
  83 on something more likely to work.
  84
  85 Further, the above method always uses a starting guess which is half
  86 the target.  This does not facilitate recovering all the roots of a
  87 function (not that Newton's method is ever really good at that).
  88
  89 =head2 The Defaults
  90
  91 Note that Perl 6 allows us to provide defaults for optional parameters.
  92 This includes providing a default fprime in the example.  For fprime, the
  93 example goes on step further and Curries approxfprime with the user supplied
  94 f.  Such things are possible in Perl 5, but simply saying .assuming is
  95 nicer.
  96
  97 =cut