Add Mueller's algorithm for root finding.

README.md

@@ -24,6 +24,7 @@ the nature of the function before selecting the root finding algorithm.
 The implemented algorithms are:
 * Bisection
 * Secant
+* Mueller
 
 ## Vectors
 

src/bbs-roots_real.adb

@@ -1,6 +1,8 @@
-with Ada.Text_IO;
-with Ada.Float_Text_IO;
+--with Ada.Text_IO;
+--with Ada.Float_Text_IO;
+with Ada.Numerics.Generic_Elementary_Functions;
 package body BBS.roots_real is
+ package elem is new Ada.Numerics.Generic_Elementary_Functions(f);
  --
  -- Bisection algorithm for finding a root.
  --
@@ -46,6 +48,10 @@ package body BBS.roots_real is
  end if;
  err := none;
  for i in 0 .. limit loop
+ if val_high = val_low then
+ err := no_solution;
+ return 0.0;
+ end if;
  mid_limit := upper - val_high*(upper - lower)/(val_high - val_low);
  val_mid := test(mid_limit);
  exit when val_mid = 0.0;
@@ -61,6 +67,49 @@ package body BBS.roots_real is
  upper := mid_limit;
  lower := mid_limit;
  end if;
+ if val_high = val_low then
+ err := no_solution;
+ return 0.0;
+ end if;
  return upper - val_high*(upper - lower)/(val_high - val_low);
  end;
+ --
+ function mueller(test : test_func; x0, x2 : in out f; limit : Positive; err : out errors) return f is
+ x1 : f := (x0 + x2)/2.0;
+ h1 : f := x1 - x0;
+ h2 : f := x2 - x1;
+ d1 : f := (test(x1) - test(x0))/h1;
+ d2 : f := (test(x2) - test(x1))/h2;
+ d_small : f := (d2 - d1) / (h2 + h1);
+ d_big : f;
+ b : f;
+ h : f;
+ e : f;
+ root : f;
+ begin
+ err := none;
+ for i in 0 .. limit loop
+ b := d2 + (h2 * d_small);
+ d_big := elem.Sqrt(b*b - (4.0*d_small*test(x2)));
+ if abs(b - d_big) < abs(b + d_big) then
+ e := b + d_big;
+ else
+ e := b - d_big;
+ end if;
+ h := -2.0*test(x2)/e;
+ root := x2 + h;
+ if h = 0.0 then
+ return root;
+ end if;
+ x0 := x1;
+ x1 := x2;
+ x2 := root;
+ h1 := x1 - x0;
+ h2 := x2 - x1;
+ d1 := (test(x1) - test(x2))/h1;
+ d2 := (test(x2) - test(x1))/h2;
+ d_small := (d2 - d1)/(h2 + h1);
+ end loop;
+ return root;
+ end;
 end;

src/bbs-roots_real.ads

@@ -29,5 +29,24 @@ package BBS.roots_real is
  -- interval containing the root. If the root is exact, the lower and upper
  -- values are equal to the returned value.
  --
+ -- This method will fail if during the process, the function values at the
+ -- upper and lower bounds are ever equal. This will cause a divide by zero
+ -- error.
+ --
  function seacant(test : test_func; lower, upper : in out f; limit : Positive; err : out errors) return f;
+ --
+ -- The Mueller algorithm uses three points to model a quadratic curve and uses
+ -- that to find a candidate root. This can potentially be used to find complex
+ -- roots, however this implementation does not.
+ --
+ -- This method will fail if the function value at the three points is equal.
+ --
+ -- In this implementation, the user provides two points and the third is
+ -- generated as the average of these two points. This keeps the call the
+ -- same as the bisection and secant functions.
+ --
+ -- Note that for this algorithm, the x0 and x2 values are not necessarily
+ -- meaningful as upper and lower bounds for the root.
+ --
+ function mueller(test : test_func; x0, x2 : in out f; limit : Positive; err : out errors) return f;
 end;