Implement 4/5th order Runge-Kutta-Fehlberg method.

src/BBS-Numerical-ode_real.adb

@@ -25,6 +25,30 @@ package body BBS.Numerical.ode_real is
  k4 := step*tf(start + step, initial + k3);
  return initial + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
  end;
+ --
+ -- 4/5th order Runge-Kutta-Fehlberg method
+ --
+ function rkf(tf : test_func; start, initial, step : f'Base; tol : out f'Base) return f'Base is
+ k1 : f'Base;
+ k2 : f'Base;
+ k3 : f'Base;
+ k4 : f'Base;
+ k5 : f'Base;
+ k6 : f'Base;
+ y1 : f'Base;
+ y2 : f'Base;
+ begin
+ k1 := step*tf(start, initial);
+ k2 := step*tf(start + step/4.0, initial + k1/4.0);
+ k3 := step*tf(start + 3.0*step/8.0, initial + 9.0*k2/32.0);
+ k4 := step*tf(start + 12.0*step/13.0, initial + (1932.0*k1 - 7200.0*k2 + 7296.0*k3)/2197.0);
+ k5 := step*tf(start + step, initial + 439.0*k1/216.0 - 8.0*k2 + 3680.0*k3/513.0 - 845.0*k4/4104.0);
+ k6 := step*tf(start + step/2.0, initial - 8.0*k1/27.0 + 2.0*k2 - 3544.0*k3/2565.0 + 1859.0*k4/4104.0 - 11.0*k5/40.0);
+ y1 := initial + 16.0*k1/135.0 + 6656.0*k3/12825.0 + 28561.0*k4/56430.0 - 9.0*k5/50.0 + 2.0*k6/55.0;
+ y2 := initial + 25.0*k1/216.0 + 1408.0*k3/2565.0 + 2197.0*k4/4104.0 - k5/5.0;
+ tol := abs(y1 - y2);
+ return y1;
+ end;
  -- --------------------------------------------------------------------------
  -- Multistep methods
  --

src/BBS-Numerical-ode_real.ads

@@ -15,6 +15,10 @@ package BBS.Numerical.ode_real is
  -- 4th order Runge-Kutta method
  --
  function rk4(tf : test_func; start, initial, step : f'Base) return f'Base;
+ --
+ -- 4/5th order Runge-Kutta-Fehlberg method
+ --
+ function rkf(tf : test_func; start, initial, step : f'Base; tol : out f'Base) return f'Base;
  -- --------------------------------------------------------------------------
  -- Multistep methods
  --

test/test.adb

@@ -9,7 +9,7 @@ with BBS.Numerical.roots_real;
 with BBS.Numerical.Statistics;
 
 procedure test is
- subtype real is Long_Float;
+ subtype real is 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);
@@ -35,6 +35,7 @@ procedure test is
  res : linreg.simple_linreg_result;
  y : real;
  yr : real;
+ yrkf : real;
  ye : real;
  ypc : real;
  t : real;
@@ -43,11 +44,12 @@ procedure test is
  y1 : real;
  y2 : real;
  y3 : real;
+ step : real;
 begin
  Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
  res := linreg.simple_linear(data);
  Ada.Text_IO.Put_Line("Linear regression result:");
- Ada.Text_IO.Put("a = ");
+ Ada.Text_IO.Put(" a = ");
  float_io.Put(res.a, 2, 3, 0);
  Ada.Text_IO.Put(", b = ");
  float_io.Put(res.b, 2, 3, 0);
@@ -75,30 +77,35 @@ begin
  Ada.Text_IO.New_Line;
  --
  Ada.Text_IO.Put_Line("Testing Differential Equations");
- Ada.Text_IO.Put_Line(" Time Euler's  Runge-Kutta Adams-Bashforth-Moulton");
+ Ada.Text_IO.Put_Line(" Time  Euler's 4th Order RK RKF Adams-Bashforth-Moulton");
  yr := 1.0;
  ye := 1.0;
+ yrkf := 1.0;
  y0 := 1.0;
  t := 0.0;
+ step := 0.05;
  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.Put(" ");
+ float_io.Put(yrkf, 2, 6, 0);
  Ada.Text_IO.Put(" ----");
  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);
+ for i in 1 .. 20 loop
+ ye := ode.euler(f2'Access, t, ye, step);
+ yr := ode.rk4(f2'Access, t, yr, step);
+ yrkf := ode.rkf(f2'Access, t, yrkf, step, tol);
  if i = 1 then
  y1 := yr;
  elsif i = 2 then
  y2 := yr;
  elsif i = 3 then
  y3 := yr;
  else
- ypc := ode.abam4(f2'Access, t, 0.1, y0, y1, y2, y3);
+ ypc := ode.abam4(f2'Access, t, step, y0, y1, y2, y3);
  y0 := y1;
  y1 := y2;
  y2 := y3;
@@ -110,13 +117,15 @@ begin
  float_io.Put(ye, 2, 6, 0);
  Ada.Text_IO.Put(" ");
  float_io.Put(yr, 2, 6, 0);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(yrkf, 2, 6, 0);
  if i > 3 then
  Ada.Text_IO.Put(" ");
  float_io.Put(ypc, 2, 6, 0);
  else
  Ada.Text_IO.Put(" ----");
  end if;
  Ada.Text_IO.New_Line;
- t := real(i)*0.1;
+ t := real(i)*step;
  end loop;
 end test;