Add 4th order Runge-Kutta for systems of differential equations.

README.md

@@ -19,6 +19,15 @@ Implemented algorithms are:
 ##  Differentiation
 
 ##  Differential Equations
+These are basic algorithms for finding numerical solutions to differential
+equations.
+
+the implemented algorithms are:
+* Euler's method
+* 4th order Runge-Kutta
+* 4/5th order Runge-Kutta-Fehlberg
+* 4th order Adams-Bashforth/Adams-Moulton
+* 4th order Runge-Kutta for systems of differential equations
 
 ##  Root Finding
 These are used to find roots (zero crossings) of functions.  There are a number

src/BBS-Numerical-ode_real.adb

@@ -74,5 +74,44 @@ package body BBS.Numerical.ode_real is
    begin
       return am4(tf, start, step, i1, i2, i3, i4);
    end;
+   -- --------------------------------------------------------------------------
+   --  Systems of differential equation methods
+   --
+   --  4th order Runge-Kutta method
+   --
+   function rk4s(tf : functs; start, initial : params; step : f'Base) return params is
+      w  : params(tf'First .. tf'Last);
+      k1 : params(tf'First .. tf'Last);
+      k2 : params(tf'First .. tf'Last);
+      k3 : params(tf'First .. tf'Last);
+      k4 : params(tf'First .. tf'Last);
+      res : params(tf'First .. tf'Last);
+   begin
+      for i in tf'Range loop
+         k1(i) := step*tf(i)(start(i), initial);
+      end loop;
+      for i in tf'Range loop
+         w(i) := initial(i) + k1(i)/2.0;
+      end loop;
+      for i in tf'Range loop
+         k2(i) := step*tf(i)(start(i) + step/2.0, w);
+      end loop;
+      for i in tf'Range loop
+         w(i) := initial(i) + k2(i)/2.0;
+      end loop;
+      for i in tf'Range loop
+         k3(i) := step*tf(i)(start(i) + step/2.0, w);
+      end loop;
+      for i in tf'Range loop
+         w(i) := initial(i) + k3(i);
+      end loop;
+      for i in tf'Range loop
+         k4(i) := step*tf(i)(start(i) + step, w);
+      end loop;
+      for i in tf'Range loop
+         res(i) := initial(i) + (k1(i) + 2.0*k2(i) + 2.0*k3(i) + k4(i))/6.0;
+      end loop;
+      return res;
+   end;
    --
 end BBS.Numerical.ode_real;

src/BBS-Numerical-ode_real.ads

@@ -5,10 +5,14 @@ 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;
+   type params is array (integer range <>) of f'Base;
+   type sys_func is access function (t : f'Base; y : params) return f'Base;
+   type functs is array (integer range <>) of sys_func;
    -- --------------------------------------------------------------------------
    --  Single step methods
    --
-   --  Euler's method
+   --  Euler's method (It is not recommended to use this.  It was included only
+   --  for illustration and debugging purposes.)
    --
    function euler(tf : test_func; start, initial, step : f'Base) return f'Base;
    --
@@ -27,4 +31,14 @@ package BBS.Numerical.ode_real is
    function ab4(tf : test_func; start, step : f'Base; i0, i1, i2, i3 : f'Base) return f'Base;
    function am4(tf : test_func; start, step : f'Base; i0, i1, i2, i3 : f'Base) return f'Base;
    function abam4(tf : test_func; start, step : f'Base; i0, i1, i2, i3 : f'Base) return f'Base;
+   -- --------------------------------------------------------------------------
+   --  Systems of differential equation methods
+   --
+   --  4th order Runge-Kutta method
+   --
+   --  Note that the arrays for tf, start, and initial must have the same bounds.
+   --
+   function rk4s(tf : functs; start, initial : params; step : f'Base) return params
+      with pre => (tf'First = start'First) and (tf'First = initial'First) and
+                  (tf'Last = start'Last) and (tf'Last = initial'Last);
 end BBS.Numerical.ode_real;

test/test.adb

@@ -26,6 +26,19 @@ procedure test is
       return -y + t + 1.0;
    end;
 
+   type param_array is array(1 .. 2) of real;
+   type sys_func is access function (t : real; y : ode.params) return real;
+
+   function f1sys(t : real; y : ode.params) return real is
+   begin
+      return -4.0*y(1) + 3.0*y(2) + 6.0;
+   end;
+
+   function f2sys(t : real; y : ode.params) return real is
+   begin
+      return -2.4*y(1) + 1.6*y(2) + 3.6;
+   end;
+
    data : linreg.data_array :=
       ((x => 1.0, y => 1.0),
        (x => 2.0, y => 1.0),
@@ -45,6 +58,10 @@ procedure test is
    y2  : real;
    y3  : real;
    step : real;
+   tsys : ode.params(1..2);
+   ysys : ode.params(1..2);
+   rsys : ode.params(1..2);
+--   func : constant func_array := (f1sys'Access, f2sys'Access);
 begin
    Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
    res := linreg.simple_linear(data);
@@ -111,6 +128,7 @@ begin
          y2 := y3;
          y3 := ypc;
       end if;
+      t := real(i)*step;
       Ada.Text_IO.Put("  ");
       float_io.Put(t, 2, 2, 0);
       Ada.Text_IO.Put("  ");
@@ -126,6 +144,31 @@ begin
          Ada.Text_IO.Put("       ----");
       end if;
       Ada.Text_IO.New_Line;
+   end loop;
+   --
+   Ada.Text_IO.Put_Line("Testing Differential Equation Systems");
+   ysys := (0.0, 0.0);
+   tsys := (0.0, 0.0);
+   step := 0.1;
+   t := 0.0;
+   Ada.Text_IO.Put("  ");
+   float_io.Put(0.0, 2, 2, 0);
+   Ada.Text_IO.Put("  ");
+   float_io.Put(ysys(1), 2, 6, 0);
+   Ada.Text_IO.Put("  ");
+   float_io.Put(ysys(2), 2, 6, 0);
+   Ada.Text_IO.New_Line;
+   for i in 1 .. 10 loop
+      rsys := ode.rk4s((1 => f1sys'Access, 2 => f2sys'Access),
+               tsys, ysys, step);
+      ysys := rsys;
       t := real(i)*step;
+      Ada.Text_IO.Put("  ");
+      float_io.Put(t, 2, 2, 0);
+      Ada.Text_IO.Put("  ");
+      float_io.Put(ysys(1), 2, 6, 0);
+      Ada.Text_IO.Put("  ");
+      float_io.Put(ysys(2), 2, 6, 0);
+      Ada.Text_IO.New_Line;
    end loop;
 end test;