Changes:

Update adaptave integration to return an estimate of the tolerance
Add package for Ordinary Differential Equations
Add Euler's method for differential equations
Add 4th order Runge-Kutta method for differential equations

src/BBS-Numerical-integration_real.adb

@@ -47,22 +47,29 @@ package body BBS.Numerical.integration_real is
  -- Integrate the provided function between the lower and upper limits using
  -- adaptive Simpson's integration.
  --
- function adapt_simpson(test : test_func; lower, upper, tol : f'Base;
- levels : in out Natural) return f'Base is
- --
- function interval(t : test_func; lower, upper, tol : f'Base;
+ function adapt_simpson(test : test_func; lower, upper : f'Base; tol : in out f'Base;
+ levels : Natural) return f'Base is
+ --
+ -- Function to recurse for adaptive Simpson's integration
+ --
+ function interval(t : test_func; lower, upper : f'Base; tol : in out f'Base;
  t_l, t_m, t_u, full : f'Base; levels : Natural) return f'Base is
  mid : constant f'Base := (upper + lower)/2.0;
  t_ml : constant f'Base := test((lower + mid)/2.0);
  t_mu : constant f'Base := test((mid + upper)/2.0);
  l : f'Base;
  r : f'Base;
+ tol1 : f'Base := tol/2.0;
+ tol2 : f'Base := tol/2.0;
  begin
  l := (upper - lower)*(t_l + 4.0*t_ml + t_m)/12.0;
  r := (upper - lower)*(t_m + 4.0*t_mu + t_u)/12.0;
  if (abs(l + r - full) > tol) and (levels > 0) then
- l := interval(test, lower, mid, tol/2.0, t_l, t_ml, t_m, l, levels - 1);
- r := interval(test, mid, upper, tol/2.0, t_m, t_mu, t_u, r, levels - 1);
+ l := interval(test, lower, mid, tol1, t_l, t_ml, t_m, l, levels - 1);
+ r := interval(test, mid, upper, tol2, t_m, t_mu, t_u, r, levels - 1);
+ tol := tol1 + tol2;
+ else
+ tol := abs(l + r - full);
  end if;
  return l + r;
  end;
@@ -76,12 +83,17 @@ package body BBS.Numerical.integration_real is
  full : constant f'Base := (upper - lower)*(t_l + 4.0*t_m + t_u)/6.0;
  l : f'Base;
  r : f'Base;
+ tol1 : f'Base := tol/2.0;
+ tol2 : f'Base := tol/2.0;
  begin
  l := (upper - lower)*(t_l + 4.0*t_ml + t_m)/12.0;
  r := (upper - lower)*(t_m + 4.0*t_mu + t_u)/12.0;
  if (abs(l + r - full) > tol) and (levels > 0) then
- l := interval(test, lower, mid, tol/2.0, t_l, t_ml, t_m, l, levels - 1);
- r := interval(test, mid, upper, tol/2.0, t_m, t_mu, t_u, r, levels - 1);
+ l := interval(test, lower, mid, tol1, t_l, t_ml, t_m, l, levels - 1);
+ r := interval(test, mid, upper, tol2, t_m, t_mu, t_u, r, levels - 1);
+ tol := tol1 + tol2;
+ else
+ tol := abs(l + r - full);
  end if;
  return l + r;
  end;

src/BBS-Numerical-integration_real.ads

@@ -19,8 +19,9 @@ package BBS.Numerical.integration_real is
  function simpson(test : test_func; lower, upper : f'Base; steps : Positive) return f'Base;
  --
  -- Integrate the provided function between the lower and upper limits using
- -- adaptive Simpson's integration.
+ -- adaptive Simpson's integration. A requested tolerance is input and
+ -- returned as an estimated value.
  --
- function adapt_simpson(test : test_func; lower, upper, tol : f'Base;
- levels : in out Natural) return f'Base;
+ function adapt_simpson(test : test_func; lower, upper : f'Base; tol : in out f'Base;
+ levels : Natural) return f'Base;
 end;

src/BBS-Numerical-ode_real.adb

@@ -0,0 +1,26 @@
+with Ada.Text_IO;
+package body BBS.Numerical.ode_real is
+ package float_io is new Ada.Text_IO.Float_IO(f'Base);
+ --
+ -- Euler's method
+ --
+ function euler(tf : test_func; start, initial, step : f'Base) return f'Base is
+ begin
+ return initial + step*tf(start, initial);
+ end;
+ --
+ -- 4th order Runge-Kutta method - single step
+ --
+ function rk4(tf : test_func; start, initial, step : f'Base) return f'Base is
+ k1 : f'Base;
+ k2 : f'Base;
+ k3 : f'Base;
+ k4 : f'Base;
+ begin
+ k1 := step*tf(start, initial);
+ k2 := step*tf(start + step/2.0, initial + k1/2.0);
+ k3 := step*tf(start + step/2.0, initial + k2/2.0);
+ k4 := step*tf(start + step, initial + k3);
+ return initial + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
+ end;
+end BBS.Numerical.ode_real;

src/BBS-Numerical-ode_real.ads

@@ -0,0 +1,17 @@
+generic
+ type F is digits <>;
+package BBS.Numerical.ode_real is
+ --
+ -- Define a type for the function to integrate.
+ --
+ type test_func is access function (t, y : f'Base) return f'Base;
+ --
+ -- Euler's method
+ --
+ function euler(tf : test_func; start, initial, step : f'Base) return f'Base;
+ --
+ -- 4th order Runge-Kutta method - single step
+ --
+ function rk4(tf : test_func; start, initial, step : f'Base) return f'Base;
+ --
+end BBS.Numerical.ode_real;

test/test.adb

@@ -3,30 +3,41 @@ with Ada.Text_IO;
 with BBS.Numerical;
 with BBS.Numerical.complex_real;
 with BBS.Numerical.integration_real;
+with BBS.Numerical.ode_real;
 with BBS.Numerical.regression;
 with BBS.Numerical.roots_real;
 with BBS.Numerical.Statistics;
 
 procedure test is
- package linreg is new BBS.Numerical.regression(Long_Float);
- package integ is new BBS.Numerical.integration_real(Long_Float);
- package float_io is new Ada.Text_IO.Float_IO(Long_Float);
- package elem is new Ada.Numerics.Generic_Elementary_Functions(Long_Float);
+ subtype real is Long_Float;
+ package linreg is new BBS.Numerical.regression(real);
+ package integ is new BBS.Numerical.integration_real(real);
+ package ode is new BBS.Numerical.ode_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 : Long_Float) return Long_Float is
+ function f1(x : real) return real is
  begin
  return (100.0/(x*x))*elem.sin(10.0/x);
  end;
 
+ function f2(t, y : real) return real is
+ begin
+ return -y + t + 1.0;
+ end;
+
  data : linreg.data_array :=
  ((x => 1.0, y => 1.0),
  (x => 2.0, y => 1.0),
  (x => 3.0, y => 2.0),
  (x => 4.0, y => 2.0),
  (x => 5.0, y => 4.0));
  res : linreg.simple_linreg_result;
- y : Long_Float;
- i : Positive;
+ y : real;
+ yr : real;
+ ye : real;
+ t : real;
+ tol : real;
 begin
  Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
  res := linreg.simple_linear(data);
@@ -40,18 +51,46 @@ begin
  Ada.Text_IO.Put(", variance = ");
  float_io.Put(res.variance, 2, 3, 0);
  Ada.Text_IO.New_Line;
+ --
  Ada.Text_IO.Put_Line("Testing integration:");
  y := integ.trapezoid(f1'Access, 1.0, 3.0, 10);
- Ada.Text_IO.Put("Trapazoid rule gives ");
+ Ada.Text_IO.Put(" Trapazoid rule gives ");
  float_io.Put(y, 2, 3, 0);
  Ada.Text_IO.New_Line;
  y := integ.simpson(f1'Access, 1.0, 3.0, 10);
- Ada.Text_IO.Put("Simpson rule gives ");
- float_io.Put(y, 2, 3, 0);
+ Ada.Text_IO.Put(" Simpson rule gives ");
+ float_io.Put(y, 2, 6, 0);
  Ada.Text_IO.New_Line;
- i := 10;
- y := integ.adapt_simpson(f1'Access, 1.0, 3.0, 1.0e-6, i);
- Ada.Text_IO.Put("Adaptive Simpson gives ");
+ tol := 1.0e-6;
+ y := integ.adapt_simpson(f1'Access, 1.0, 3.0, tol, 8);
+ Ada.Text_IO.Put(" Adaptive Simpson gives ");
  float_io.Put(y, 2, 6, 0);
+ Ada.Text_IO.Put(" with estimated tolerance ");
+ float_io.Put(tol, 2, 6, 0);
+ Ada.Text_IO.New_Line;
+ --
+ Ada.Text_IO.Put_Line("Testing Differential Equations");
+ Ada.Text_IO.Put_Line(" Time Euler's Runge-Kutta");
+ yr := 1.0;
+ ye := 1.0;
+ t := 0.0;
+ Ada.Text_IO.Put(" ");
+ float_io.Put(0.0, 2, 2, 0);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(ye, 2, 6, 0);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(yr, 2, 6, 0);
  Ada.Text_IO.New_Line;
+ for i in 1 .. 10 loop
+ ye := ode.euler(f2'Access, t, ye, 0.1);
+ yr := ode.rk4(f2'Access, t, yr, 0.1);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(t, 2, 2, 0);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(ye, 2, 6, 0);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(yr, 2, 6, 0);
+ Ada.Text_IO.New_Line;
+ t := real(i)*0.1;
+ end loop;
 end test;