Add test routine for root finding and minor change to mueller to help…

… avoid NaN.

.gitignore

@@ -9,4 +9,5 @@ lib/*
 # Test binaries
 test_ode
 test_integ
+test_root
 test_stats

src/BBS-Numerical-roots_real.adb

@@ -1,16 +1,16 @@
---with Ada.Text_IO;
---with Ada.Float_Text_IO;
 with Ada.Numerics.Generic_Elementary_Functions;
+with Ada.Text_IO;
 package body BBS.Numerical.roots_real is
  package elem is new Ada.Numerics.Generic_Elementary_Functions(f);
+ package float_io is new Ada.Text_IO.Float_IO(f'Base);
  --
  -- Bisection algorithm for finding a root.
  --
- function bisection(test : test_func; lower, upper : in out f; limit : Positive; err : out errors) return f is
- val_high : f;
- val_low : f;
- val_mid : f;
- mid_limit : f;
+ function bisection(test : test_func; lower, upper : in out f'Base; limit : Positive; err : out errors) return f'Base is
+ val_high : f'Base;
+ val_low : f'Base;
+ val_mid : f'Base;
+ mid_limit : f'Base;
  begin
  val_high := test(upper);
  val_low := test(lower);
@@ -34,11 +34,11 @@ package body BBS.Numerical.roots_real is
  return (lower + upper)/2.0;
  end;
  --
- function seacant(test : test_func; lower, upper : in out f; limit : Positive; err : out errors) return f is
- val_high : f;
- val_low : f;
- val_mid : f;
- mid_limit : f;
+ function seacant(test : test_func; lower, upper : in out f'Base; limit : Positive; err : out errors) return f'Base is
+ val_high : f'Base;
+ val_low : f'Base;
+ val_mid : f'Base;
+ mid_limit : f'Base;
  begin
  val_high := test(upper);
  val_low := test(lower);
@@ -74,20 +74,20 @@ package body BBS.Numerical.roots_real is
  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;
- step1 : f := x1 - x0;
- step2 : f := x2 - x1;
- delta1 : f := (test(x1) - test(x0))/step1;
- delta2 : f := (test(x2) - test(x1))/step2;
- d_small : f := (delta2 - delta1) / (step2 + step1);
- d_big : f;
- discriminant : f;
- b : f;
- value : f;
- e : f;
- root : f;
- temp : f;
+ function mueller(test : test_func; x0, x2 : in out f'Base; limit : Positive; err : out errors) return f'Base is
+ x1 : f'Base := (x0 + x2)/2.0;
+ step1 : f'Base := x1 - x0;
+ step2 : f'Base := x2 - x1;
+ delta1 : f'Base := (test(x1) - test(x0))/step1;
+ delta2 : f'Base := (test(x2) - test(x1))/step2;
+ d_small : f'Base := (delta2 - delta1) / (step2 + step1);
+ d_big : f'Base;
+ discriminant : f'Base;
+ b  : f'Base;
+ value  : f'Base;
+ e  : f'Base;
+ root  : f'Base;
+ temp  : f'Base;
  begin
  err := none;
  for i in 0 .. limit loop
@@ -105,7 +105,7 @@ package body BBS.Numerical.roots_real is
  end if;
  value := -2.0*test(x2)/e;
  root := x2 + value;
- if value = 0.0 then
+ if (value = 0.0) or (x0 = x1) or (root = x2) then
  return root;
  end if;
  x0 := x1;

src/BBS-Numerical-roots_real.ads

@@ -13,7 +13,7 @@ package BBS.Numerical.roots_real is
  -- The lower and upper values are updated during the iterations to provide an
  -- interval containing the root.
  --
- function bisection(test : test_func; lower, upper : in out f; limit : Positive; err : out errors) return f;
+ function bisection(test : test_func; lower, upper : in out f'Base; limit : Positive; err : out errors) return f'Base;
  --
  -- If there is an odd number of roots between the lower and upper limits,
  -- the secant method will always converge. Each stage of the iteration
@@ -33,7 +33,7 @@ package BBS.Numerical.roots_real is
  -- 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;
+ function seacant(test : test_func; lower, upper : in out f'Base; limit : Positive; err : out errors) return f'Base;
  --
  -- The Mueller algorithm uses three points to model a quadratic curve and uses
  -- that to find a candidate root. Unlike the bisection method, Mueller's
@@ -51,5 +51,5 @@ package BBS.Numerical.roots_real is
  -- meaningful as upper and lower bounds for the root, except that they are
  -- both set equal to the return value if the root is exact.
  --
- function mueller(test : test_func; x0, x2 : in out f; limit : Positive; err : out errors) return f;
+ function mueller(test : test_func; x0, x2 : in out f'Base; limit : Positive; err : out errors) return f'Base;
 end;

test/test_ode.adb

@@ -50,7 +50,7 @@ begin
  Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
  --
  Ada.Text_IO.Put_Line("Testing Differential Equations");
- Ada.Text_IO.Put_Line(" Time Euler's 4th Order RK RKF Adams-Bashforth-Moulton");
+ Ada.Text_IO.Put_Line(" Time Euler's 4th Order RK  RKF RKF Tol Adams-Bashforth-Moulton");
  yr := 1.0;
  ye := 1.0;
  yrkf := 1.0;
@@ -65,6 +65,8 @@ begin
  float_io.Put(yr, 2, 6, 0);
  Ada.Text_IO.Put(" ");
  float_io.Put(yrkf, 2, 6, 0);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(0.0, 2, 6, 3);
  Ada.Text_IO.Put(" ----");
  Ada.Text_IO.New_Line;
  for i in 1 .. 20 loop
@@ -93,6 +95,8 @@ begin
  float_io.Put(yr, 2, 6, 0);
  Ada.Text_IO.Put(" ");
  float_io.Put(yrkf, 2, 6, 0);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(tol, 2, 6, 3);
  if i > 3 then
  Ada.Text_IO.Put(" ");
  float_io.Put(ypc, 2, 6, 0);
@@ -103,6 +107,7 @@ begin
  end loop;
  --
  Ada.Text_IO.Put_Line("Testing Differential Equation Systems");
+ Ada.Text_IO.Put_Line(" Time RK4-a RK4-b");
  ysys := (0.0, 0.0);
  step := 0.1;
  t := 0.0;

test/test_root.adb

@@ -0,0 +1,62 @@
+with Ada.Numerics.Generic_Elementary_Functions;
+with Ada.Text_IO;
+with BBS.Numerical;
+with BBS.Numerical.roots_real;
+
+procedure test_root is
+ subtype real is Long_Float;
+ package root is new BBS.Numerical.roots_real(real);
+ package float_io is new Ada.Text_IO.Float_IO(real);
+ package elem is new Ada.Numerics.Generic_Elementary_Functions(real);
+
+ function f1(x : real) return real is
+ begin
+ return (100.0/(x*x))*elem.sin(10.0/x);
+ end;
+
+ function f2(x : real) return real is
+ begin
+ return x*x*x + 4.0*x*x - 10.0;
+ end;
+
+ r : real;
+ l : real;
+ u : real;
+ err : root.errors;
+begin
+ Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
+ --
+ Ada.Text_IO.Put_Line("Testing root finding:");
+ l := 1.0;
+ u := 2.0;
+ r := root.bisection(f2'Access, l, u, 13, err);
+ Ada.Text_IO.Put(" Bisection gives root at");
+ float_io.Put(r, 2, 9, 0);
+ Ada.Text_IO.Put(", in range ");
+ float_io.Put(l, 2, 6, 0);
+ Ada.Text_IO.Put(" to ");
+ float_io.Put(u, 2, 6, 0);
+ Ada.Text_IO.Put_Line(", with error code " & root.errors'Image(err));
+ --
+ l := 1.0;
+ u := 2.0;
+ r := root.seacant(f2'Access, l, u, 13, err);
+ Ada.Text_IO.Put(" Seacant gives root at");
+ float_io.Put(r, 2, 9, 0);
+ Ada.Text_IO.Put(", in range ");
+ float_io.Put(l, 2, 6, 0);
+ Ada.Text_IO.Put(" to ");
+ float_io.Put(u, 2, 6, 0);
+ Ada.Text_IO.Put_Line(", with error code " & root.errors'Image(err));
+ --
+ l := 1.0;
+ u := 2.0;
+ r := root.mueller(f2'Access, l, u, 13, err);
+ Ada.Text_IO.Put(" Mueller gives root at");
+ float_io.Put(r, 2, 9, 0);
+ Ada.Text_IO.Put(", in range ");
+ float_io.Put(l, 2, 6, 0);
+ Ada.Text_IO.Put(" to ");
+ float_io.Put(u, 2, 6, 0);
+ Ada.Text_IO.Put_Line(", with error code " & root.errors'Image(err));
+end test_root;

test/test_root.gpr

@@ -0,0 +1,15 @@
+with "../BBS_Numerical.gpr";
+
+project test_root is
+ for languages use ("Ada");
+ for Source_Dirs use (".");
+ for Object_Dir use "../obj";
+ for Main use ("test_root.adb");
+ for Exec_Dir use ".";
+ package builder is
+ for switches ("Ada") use ("-s");
+ end builder;
+ package compiler is
+ for switches ("Ada") use ("-g", "-gnateE", "-gnata");
+ end compiler;
+end test_root;