Updates to mueller's method for root finding. Minor updates to polyno…

…mial

division.  Add utility function to print polynomials.

src/BBS-Numerical-polynomial_real.adb

@@ -95,13 +95,14 @@ package body BBS.Numerical.polynomial_real is
  -- u/v => q, r
  --
  procedure divide(u, v : poly; q : out poly; r : out poly) is
- temp : poly(u'Range) := u;
+ temp : poly := u;
+ v_last : constant f'Base := v(v'Last);
  begin
  q := (others => 0.0);
  r := (others => 0.0);
  --
  for k in reverse 0 .. u'Last - v'Last loop
- q(k) := temp(v'Last + k)/v(v'Last);
+ q(k) := temp(v'Last + k)/v_last;
  for j in k .. v'Last + k - 1 loop
  temp(j) := temp(j) - q(k)*v(j-k);
  end loop;
@@ -149,6 +150,19 @@ package body BBS.Numerical.polynomial_real is
  return limit;
  end;
  --
+ -- Utility print procedure for debugging
+ --
+ procedure print(p : in poly; fore, aft, exp : Natural) is
+ begin
+ for i in reverse p'Range loop
+ if p(i) >= 0.0 then
+ Ada.Text_IO.Put("+");
+ end if;
+ float_io.put(p(i), fore, aft, exp);
+ Ada.Text_IO.Put("*X^" & Natural'Image(i));
+ end loop;
+ end;
+ --
  -- Basic calculus
  --
  function integrate(p : poly; c : f'Base) return poly is

src/BBS-Numerical-polynomial_real.ads

@@ -50,6 +50,10 @@ package BBS.Numerical.polynomial_real is
  function order(p : poly) return Natural
  with pre => (p'First = 0);
  --
+ -- Utility print procedure for debugging
+ --
+ procedure print(p : poly; fore, aft, exp : Natural);
+ --
  -- Basic calculus
  --
  function integrate(p : poly; c : f'Base) return poly

src/BBS-Numerical-roots_complex.adb

@@ -21,7 +21,7 @@ package body BBS.Numerical.roots_complex is
  root : cmplx.complex;
  temp : cmplx.complex;
  nTwo : constant cmplx.complex := (-2.0)*cmplx.one;
- atemp  : ada_cmplx.Complex;
+ atemp : ada_cmplx.Complex;
  begin
  err := none;
  for i in 0 .. limit loop

src/BBS-Numerical-roots_complex.ads

@@ -23,6 +23,10 @@ package BBS.Numerical.roots_complex 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.
  --
+ --
+ -- Note that the sucess may be sensitive to the choice of x0 and x2. If you
+ -- know that a root exists and get a no_solution error, try different values.
+ --
  function mueller(test : test_func; x0, x2 : in out cmplx.complex;
  limit : in out Positive; err : out errors) return cmplx.complex;
 end BBS.Numerical.roots_complex;

src/BBS-Numerical-roots_real.adb

@@ -86,16 +86,18 @@ package body BBS.Numerical.roots_real is
  b : f'Base;
  value : f'Base;
  e : f'Base;
- root : f'Base;
+ root : f'Base := x2;
  temp : f'Base;
+ epsilon : constant f'Base := 1.0e-10; -- Adjust as needed. Maybe make a parameter
  begin
  err := none;
- for i in 0 .. limit loop
+ for i in 1 .. limit loop
  b := delta2 + (step2 * d_small);
  discriminant := b*b - (4.0*d_small*test(x2));
  if discriminant < 0.0 then
  err := no_solution;
- return 0.0;
+ limit := i;
+ return root;
  end if;
  d_big := elem.Sqrt(discriminant);
  if abs(b - d_big) < abs(b + d_big) then
@@ -105,7 +107,7 @@ package body BBS.Numerical.roots_real is
  end if;
  value := -2.0*test(x2)/e;
  root := x2 + value;
- if (value = 0.0) or (x0 = x1) or (root = x2) then
+ if (abs(value) <= epsilon) or (x0 = x1) or (root = x2) then
  limit := i;
  return root;
  end if;

src/BBS-Numerical-roots_real.ads

@@ -51,5 +51,8 @@ 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.
  --
+ -- Note that the sucess may be sensitive to the choice of x0 and x2. If you
+ -- know that a root exists and get a no_solution error, try different values.
+ --
  function mueller(test : test_func; x0, x2 : in out f'Base; limit : in out Positive; err : out errors) return f'Base;
 end;

test/test_poly.adb

@@ -11,17 +11,6 @@ procedure test_poly is
  package root is new BBS.Numerical.roots_real(real);
  package float_io is new Ada.Text_IO.Float_IO(real);
 
- procedure put_poly(p : poly.poly) is
- begin
- for i in reverse p'Range loop
- if p(i) >= 0.0 then
- Ada.Text_IO.Put("+");
- end if;
- float_io.put(p(i), 1, 6, 0);
- Ada.Text_IO.Put("*X^" & Natural'Image(i));
- end loop;
- end;
-
  p1 : poly.poly := (-1.0, 2.0, 3.0);
  p2 : poly.poly := (3.0, 2.0, 1.0);
  p3 : poly.poly(0 .. 2);
@@ -31,43 +20,63 @@ procedure test_poly is
  p7 : poly.poly(0 .. 2);
  p8 : poly.poly(0 .. 2);
  p9 : poly.poly(0 .. 1);
+ b1 : poly.poly := (1.0, 1.0);
+ b2 : poly.poly := (2.0, 1.0);
+ b3 : poly.poly := (3.0, 1.0);
+ b4 : poly.poly := (4.0, 1.0);
+ b5 : poly.poly(0 .. 1);
+ b6 : poly.poly(0 .. 1);
+ d0 : poly.poly(0 .. 4);
+ d1 : poly.poly(0 .. 3);
+ d2 : poly.poly(0 .. 2);
+ d3 : poly.poly(0 .. 1);
  x : real;
  r : real;
  l : real;
  u : real;
  err : root.errors;
  iter : Positive;
 
- function test(x : real) return real is
+ function t0(x : real) return real is
+ begin
+ return poly.evaluate(d0, x);
+ end;
+
+ function t1(x : real) return real is
  begin
- return poly.evaluate(p4, x);
+ return poly.evaluate(d1, x);
+ end;
+
+ function t2(x : real) return real is
+ begin
+ return poly.evaluate(d2, x);
  end;
 begin
  Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
  Ada.Text_IO.Put_Line("Basic polynomial operations");
  Ada.Text_IO.Put(" p1 = ");
- put_poly(p1);
+ poly.print(p1, 1, 2, 0);
  Ada.Text_IO.New_Line;
  Ada.Text_IO.Put(" p2 = ");
- put_poly(p2);
+ poly.print(p2, 1, 2, 0);
  Ada.Text_IO.New_Line;
- p3 := p1 + p2;
- Ada.Text_IO.Put(" p3 = p1+p2 = ");
- put_poly(p3);
+ p3 := p1 - p2;
+ Ada.Text_IO.Put(" p3 = p1-p2 = ");
+ poly.print(p3, 1, 2, 0);
  Ada.Text_IO.New_Line;
  p4 := p1*p2;
  Ada.Text_IO.Put(" p4 = p1*p2 = ");
- put_poly(p4);
+ poly.print(p4, 1, 2, 0);
  Ada.Text_IO.New_Line;
  p5 := (-1.0)*p2;
  Ada.Text_IO.Put(" p5 = -p2 = ");
- put_poly(p5);
+ poly.print(p5, 1, 2, 0);
  Ada.Text_IO.New_Line;
  poly.divide(p4, p1, p8, p9);
  Ada.Text_IO.Put(" p8 = p4/p1 = ");
- put_poly(p8);
+ poly.print(p8, 1, 2, 0);
  Ada.Text_IO.Put(", remainder p9 = ");
- put_poly(p9);
+ poly.print(p9, 1, 2, 0);
  Ada.Text_IO.New_Line;
  --
  Ada.Text_IO.Put_Line("Evaluations of polynomials");
@@ -88,29 +97,72 @@ begin
  Ada.Text_IO.New_Line;
  end loop;
  --
- Ada.Text_IO.Put_Line("Find a root of P4");
- l := 0.0;
- u := 1.0;
+ d0 := b1*b2*b3*b4;
+ Ada.Text_IO.Put("Find roots of ");
+ poly.print(d0, 1, 2, 0);
+ Ada.Text_IO.New_Line;
+ l := -5.1;
+ u := -3.1;
+ iter := 20;
+ r := root.mueller(t0'Access, l, u, iter, err);
+ Ada.Text_IO.Put(" After " & Positive'image(iter) & " iterations, 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));
+ b5 := (0 => -r, 1 => 1.0);
+ poly.divide(d0, b5, d1, b6);
+ Ada.Text_IO.Put("Find roots of ");
+ poly.print(d1, 1, 2, 0);
+ Ada.Text_IO.New_Line;
+ l := -4.1;
+ u := -0.9;
+ iter := 13;
+ r := root.mueller(t1'Access, l, u, iter, err);
+ Ada.Text_IO.Put(" After " & Positive'image(iter) & " iterations, 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));
+ b5 := (0 => -r, 1 => 1.0);
+ poly.divide(d1, b5, d2, b6);
+ Ada.Text_IO.Put("Find roots of ");
+ poly.print(d2, 1, 2, 0);
+ Ada.Text_IO.New_Line;
+ l := -5.1;
+ u := 0.0;
  iter := 13;
- r := root.mueller(test'Access, l, u, iter, err);
+ r := root.mueller(t2'Access, l, u, iter, err);
  Ada.Text_IO.Put(" After " & Positive'image(iter) & " iterations, 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));
+ b5 := (0 => -r, 1 => 1.0);
+ poly.divide(d2, b5, d3, b6);
+ Ada.Text_IO.Put(" Last root at ");
+ poly.print(d3, 1, 2, 0);
+ Ada.Text_IO.New_Line;
+ Ada.Text_IO.Put(" Remainder ");
+ poly.print(b6, 1, 2, 0);
+ Ada.Text_IO.New_Line;
  --
  Ada.Text_IO.Put_Line("Integrals and derivatives");
  p6 := poly.integrate(p1, 1.0);
  Ada.Text_IO.Put(" P1 = ");
- put_poly(p1);
+ poly.print(p1, 1, 2, 0);
  Ada.Text_IO.New_Line;
  Ada.Text_IO.Put(" p6 = Integral of p1 = ");
- put_poly(p6);
+ poly.print(p6, 1, 2, 0);
  Ada.Text_IO.New_Line;
  p7 := poly.derivative(p6);
  Ada.Text_IO.Put(" p7 = Derivative of P6 = ");
- put_poly(p7);
+ poly.print(p7, 1, 2, 0);
  Ada.Text_IO.New_Line;
 end test_poly;