Add composite simpson's rule for integration.

README.md

@@ -13,6 +13,7 @@ for some of the operations, if that's important to you.
 ## Integration
 Implemented algorithms are:
 * Composite trapezoid
+* Composite simpson
 
 ## Differentiation
 

src/bbs-integration_real.adb

@@ -1,11 +1,11 @@
 package body BBS.integration_real is
  --
  function trapezoid(test : test_func; lower, upper : f; steps : Positive) return f is
+ step_size : constant f := (upper - lower)/f(steps);
  accumulator : f := 0.0;
  last_func : f := test(lower);
  next_func : f;
  base : f := lower;
- step_size : f := (upper - lower)/f(steps);
  begin
  for i in 0 .. steps - 1 loop
  base := base + step_size;
@@ -16,4 +16,22 @@ package body BBS.integration_real is
  return accumulator;
  end;
  --
+ function simpson(test : test_func; lower, upper : f; steps : Positive) return f is
+ step_size : constant f := (upper - lower)/f(2*steps);
+ ends : constant f := test(lower) + test(upper);
+ sum1 : f := 0.0;
+ sum2 : f := 0.0;
+ base : f := lower;
+ begin
+ for i in 1 .. steps - 1 loop
+ base := base + step_size;
+ sum1 := sum1 + test(base);
+ base := base + step_size;
+ sum2 := sum2+ test(base);
+ end loop;
+ base := base + step_size;
+ sum1 := sum1 + test(base);
+ return step_size*(ends + 2.0*sum2 + 4.0*sum1)/3.0;
+ end;
+ --
 end;

src/bbs-integration_real.ads

@@ -7,4 +7,12 @@ package BBS.integration_real is
  -- the composite trapezoid method with the specified number of steps.
  --
  function trapezoid(test : test_func; lower, upper : f; steps : Positive) return f;
+ --
+ -- Integrate the provided function between the lower and upper limits using
+ -- the composite simpson's rule with the specified number of steps. Note
+ -- that since simpson's rule evaluates the function in the midpoint of a
+ -- segment, the effective number of steps is doubled.
+ --
+ function simpson(test : test_func; lower, upper : f; steps : Positive) return f;
+ --
 end;