Add routines for calculating numerical derivatives.

.gitignore

@@ -7,8 +7,9 @@ lib/*
 *.ali
 
 # Test binaries
-test_ode
+test_derive
 test_integ
+test_ode
 test_poly
 test_root
 test_stats

README.md

@@ -17,6 +17,9 @@ Implemented algorithms are:
 * Adaptive Simpson's
 
 ##  Differentiation
+* 2 point method
+* 3 point methods
+* 5 point methods
 
 ##  Differential Equations
 These are basic algorithms for finding numerical solutions to differential
@@ -47,6 +50,7 @@ Implemented operations include
 * Addition and subtraction of polynomials
 * Multiplication of polynomial by a scalar
 * Multiplication of polynomial by polynomial
+* Division of polynomial by polynomial
 * Evaluation of polynomial
 * Integral of polynomial
 * Derivative of polynomial

src/BBS-Numerical-derivative_real.adb

@@ -0,0 +1,38 @@
+with Ada.Text_IO;
+package body BBS.Numerical.derivative_real is
+   package float_io is new Ada.Text_IO.Float_IO(f'Base);
+   -- --------------------------------------------------------------------------
+   --
+   --  Two point formula.  Use h > 0 for forward-difference and h < 0
+   --  for backward-difference.  Derivative calculated at point x.
+   --
+   function pt2(f1 : test_func; x, h : f'Base) return f'Base is
+   begin
+      return (f1(x + h) - f1(x))/h;
+   end;
+   --
+   --  Three point formulas.
+   --
+   function pt3a(f1 : test_func; x, h : f'Base) return f'Base is
+   begin
+      return (-3.0*f1(x) + 4.0*f1(x+h) - f1(x+2.0*h))/(2.0*h);
+   end;
+   --
+   function pt3b(f1 : test_func; x, h : f'Base) return f'Base is
+   begin
+      return (f1(x+h) - f1(x-h))/(2.0*h);
+   end;
+   --
+   --  Five point formulas.
+   --
+   function pt5a(f1 : test_func; x, h : f'Base) return f'Base is
+   begin
+      return (f1(x - 2.0*h) - 8.0*f1(x - h) + 8.0*f1(x + h) - f1(x + 2.0*h))/(12.0*h);
+   end;
+   --
+   function pt5b(f1 : test_func; x, h : f'Base) return f'Base is
+   begin
+      return (-25.0*f1(x) + 48.0*f1(x + h) - 36.0*f1(x + 2.0*h) + 16.0*f1(x + 3.0*h) - 3.0*f1(x + 4.0*h))/(12.0*h);
+   end;
+   --
+end BBS.Numerical.derivative_real;

src/BBS-Numerical-derivative_real.ads

@@ -0,0 +1,34 @@
+generic
+  type F is digits <>;
+package BBS.Numerical.derivative_real is
+   -- --------------------------------------------------------------------------
+   --  Define a type for the function to integrate.
+   --
+   type test_func is access function (x : f'Base) return f'Base;
+   -- --------------------------------------------------------------------------
+   --  WARNING:
+   --  The calculations here may involve adding small numbers to large
+   --  numbers and taking the difference of two nearly equal numbers.
+   --  The world of computers is not the world of mathematics where
+   --  numbers have infinite precision.  If you aren't careful, you
+   --  can get into a situation where (x + h) = x, or f(x) = f(x +/- h).
+   --  For example assume that the float type has 6 digits.  Then, if
+   --  x is 1,000,000 and h is 1, adding x and h is a wasted operation.
+   --
+   --  Two point formula.  Use h > 0 for forward-difference and h < 0
+   --  for backward-difference.  Derivative calculated at point x.  While
+   --  this is the basis for defining the derivative in calculus, it
+   --  generally shouldn't be used.
+   --
+   function pt2(f1 : test_func; x, h : f'Base) return f'Base;
+   --
+   --  Three point formulas.
+   --
+   function pt3a(f1 : test_func; x, h : f'Base) return f'Base;
+   function pt3b(f1 : test_func; x, h : f'Base) return f'Base;
+   --
+   --  Five point formulas.
+   --
+   function pt5a(f1 : test_func; x, h : f'Base) return f'Base;
+   function pt5b(f1 : test_func; x, h : f'Base) return f'Base;
+end BBS.Numerical.derivative_real;

test/test_derive.adb

@@ -0,0 +1,56 @@
+with Ada.Numerics.Generic_Elementary_Functions;
+with Ada.Text_IO;
+with BBS.Numerical;
+with BBS.Numerical.derivative_real;
+
+procedure test_derive is
+   subtype real is Float;
+   package derive is new BBS.Numerical.derivative_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 elem.log(x);
+   end;
+
+   procedure print(f : derive.test_func; x, h : real) is
+      y : real;
+   begin
+      float_io.Put(h, 1, 4, 0);
+      Ada.Text_IO.Put("  ");
+      y := derive.pt2(f, x, h);
+      float_io.Put(y, 1, 8, 0);
+      Ada.Text_IO.Put("  ");
+      y := derive.pt3a(f, x, h);
+      float_io.Put(y, 1, 8, 0);
+      Ada.Text_IO.Put("  ");
+      y := derive.pt3b(f, x, h);
+      float_io.Put(y, 1, 8, 0);
+      Ada.Text_IO.Put("  ");
+      y := derive.pt5a(f, x, h);
+      float_io.Put(y, 1, 8, 0);
+      Ada.Text_IO.Put("  ");
+      y := derive.pt5b(f, x, h);
+      float_io.Put(y, 1, 8, 0);
+      Ada.Text_IO.New_Line;
+   end;
+
+begin
+   Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
+   --
+   Ada.Text_IO.Put_Line("  Various derivative functions");
+   Ada.Text_IO.Put_Line("  h       pt2         pt3a        pt3b        pt5a        pt5b");