Numerical integration

Numerical integration
You are encouraged to solve this task according to the task description, using any language you may know.

Write functions to calculate the definite integral of a function (f(x)) using rectangular (left, right, and midpoint), trapezium, and Simpson's methods. Your functions should take in the upper and lower bounds (a and b) and the number of approximations to make in that range (n). Assume that your example already has a function that gives values for f(x).

Simpson's method is defined by the following pseudocode:

h = (b - a) / n
sum1 = sum2 = 0

loop on i from 0 to (n - 1)
 sum1 = sum1 + f(a + h * i + h / 2)

loop on i from 1 to (n - 1)
 sum2 = sum2 + f(a + h * i)

answer = (h / 6) * (f(a) + f(b) + 4 * sum1 + 2 * sum2)

What it will be doing is computing integrals for multiple quadratics (like the one shown in Wikipedia) and summing them.

Ada

This solution creates a generic package into which the function F(X) is passed during generic instantiation. The first part is the package specification. The second part is the package body.

generic
   with function F(X : Long_Float) return Long_Float;
package Integrate is
   function  Left_Rect(A, B, N : Long_Float) return Long_Float;
   function Right_Rect(A, B, N : Long_Float) return Long_Float;
   function   Mid_Rect(A, B, N : Long_Float) return Long_Float;
   function  Trapezium(A, B, N : Long_Float) return Long_Float;
   function    Simpson(A, B, N : Long_Float) return Long_Float;
end Integrate;

package body Integrate is

   ---------------
   -- Left_Rect --
   ---------------

   function Left_Rect (A, B, N : Long_Float) return Long_Float is
      H : constant Long_Float := (B - A) / N;
      Sum : Long_Float := 0.0;
      X : Long_Float := A;
   begin
      while X <= B - H loop
         Sum := Sum + (H * F(X));
         X := X + H;
      end loop;
      return Sum;
   end Left_Rect; 

   ----------------
   -- Right_Rect --
   ---------------- 

   function Right_Rect (A, B, N : Long_Float) return Long_Float is
      H : constant Long_Float := (B - A) / N;
      Sum : Long_Float := 0.0;
      X : Long_Float := A + H;
   begin
      while X <= B - H loop
         Sum := Sum + (H * F(X));
         X := X + H;
      end loop;
      return Sum;
   end Right_Rect;

   --------------
   -- Mid_Rect --
   --------------

   function Mid_Rect (A, B, N : Long_Float) return Long_Float is
      H : constant Long_Float := (B - A) / N;
      Sum : Long_Float := 0.0;
      X : Long_Float := A;
   begin
      while X <= B - H loop
         Sum := Sum + (H / 2.0) * (F(X) + F(X + H));
         X := X + H;
      end loop;
      return Sum;
   end Mid_Rect;

   ---------------
   -- Trapezium --
   --------------- 

   function Trapezium (A, B, N : Long_Float) return Long_Float is
      H : constant Long_Float := (B - A) / N;
      Sum : Long_Float := F(A) + F(B);
      X : Long_Float := 1.0;
   begin
      while X <= N - 1.0 loop
         Sum := Sum + 2.0 * F(A + X * (B - A) / N);
         X := X + 1.0;
      end loop;
      return (B - A) / (2.0 * N) * Sum;
   end Trapezium; 

   -------------
   -- Simpson --
   ------------- 

   function Simpson (A, B, N : Long_Float) return Long_Float is
      H : constant Long_Float := (B - A) / N;
      Sum1 : Long_Float := 0.0;
      Sum2 : Long_Float := 0.0;
      Limit : Integer := Integer(N) - 1;
   begin
      for I in 0..Limit loop
         Sum1 := Sum1 + F(A + H * Long_Float(I) + H / 2.0);
      end loop;
      for I in 1..Limit loop
         Sum2 := Sum2 + F(A + H * Long_Float(I));
      end loop;
      return H / 6.0 * (F(A) + F(B) + 4.0 * Sum1 + 2.0 * Sum2);
   end Simpson;

end Integrate;

ALGOL 68

MODE F = PROC(LONG REAL)LONG REAL;

###############
## left rect ##
###############

PROC left rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum:= 0;
   LONG REAL x:= a;
   WHILE x <= b - h DO
      sum := sum + (h * f(x));
      x +:= h
   OD;
   sum
END # left rect #; 

#################
## right rect  ##
#################

PROC right rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum:= 0;
   LONG REAL x:= a + h;
   WHILE x <= b - h DO
      sum := sum + (h * f(x));
      x +:= h
   OD;
   sum
END # right rect #;

###############
## mid rect  ##
###############

PROC mid rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum:= 0;
   LONG REAL x:= a;
   WHILE x <= b - h DO
      sum := sum + (h / 2) * (f(x) + f(x + h));
      x +:= h
   OD;
   sum
END # mid rect #;

###############
## trapezium ##
############### 

PROC trapezium = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum:= f(a) + f(b);
   LONG REAL x:= 1;
   WHILE x <= n - 1 DO
      sum := sum + 2 * f(a + x * h );
      x +:= 1
   OD;
   (b - a) / (2 * n) * sum
END # trapezium #; 

#############
## simpson ##
############# 

PROC simpson = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum1:= 0;
   LONG REAL sum2:= 0;
   INT limit:= n - 1;
   FOR i FROM 0 TO limit DO
      sum1 := sum1 + f(a + h * LONG REAL(i) + h / 2)
   OD;
   FOR i FROM 1 TO limit DO
      sum2 +:= f(a + h * LONG REAL(i))
   OD;
   h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END # simpson #;
SKIP

BASIC

Works with: QuickBasic version 4.5

Translation of: Java

FUNCTION leftRect(a, b, n)
     h = (b - a) / n
     sum = 0
     FOR x = a TO b - h STEP h
        sum = sum + h * (f(x))
     NEXT x
     leftRect = sum
END FUNCTION

FUNCTION rightRect(a, b, n)
     h = (b - a) / n
     sum = 0
     FOR x = a + h TO b - h STEP h
        sum = sum + h * (f(x))
     NEXT x
     rightRect = sum
END FUNCTION

FUNCTION midRect(a, b, n)
     h = (b - a) / n
     sum = 0
     FOR x = a TO b - h STEP h
        sum = sum + (h / 2) * (f(x) + f(x + h))
     NEXT x
     midRect = sum
END FUNCTION

FUNCTION trap(a, b, n)
     h = (b - a) / n
     sum = f(a) + f(b)
     FOR i = 1 TO n-1
        sum = sum + 2 * f((a + i * h))
     NEXT i
     trap = h / 2 * sum
END FUNCTION

FUNCTION simpson(a, b, n)
     h = (b - a) / n
     sum1 = 0
     sum2 = 0

     FOR i = 0 TO n-1
        sum1 = sum + f(a + h * i + h / 2)
     NEXT i

     FOR i = 1 TO n - 1
        sum2 = sum2 + f(a + h * i)
     NEXT i

     simpson = h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END FUNCTION

C++

Due to their similarity, it makes sense to make the integration method a policy.

// the integration routine
template<typename Method, typename F, typename Float>
 double integrate(F f, Float a, Float b, int steps, Method m)
{
  double s = 0;
  double h = (b-a)/steps;
  for (int i = 0; i < steps; ++i)
    s += m(f, a + h*i, h);
  return h*s;
}

// methods
class rectangular
{
public:
  enum position_type { left, middle, right };
  rectangular(position_type pos): position(pos) {}
  template<typename F, typename Float>
   double operator()(F f, Float x, Float h)
  {
    switch(position)
    {
    case left:
      return f(x);
    case middle:
      return f(x+h/2);
    case right:
      return f(x+h);
    }
  }
private:
  position_type position;
};

class trapezium
{
public:
  template<typename F, typename Float>
   double operator()(F f, Float x, Float h)
  {
    return (f(x) + f(x+h))/2;
  }
};

class simpson
{
public:
  template<typename F, typename Float>
   double operator()(F f, Float x, Float h)
  {
    return (f(x) + 4*f(x+h/2) + f(x+h))/6;
  }
};

// sample usage
double f(double x) { return x*x; )

// inside a function somewhere:
double rl = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::left));
double rm = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::middle));
double rr = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::right));
double t  = integrate(f, 0.0, 1.0, 10, trapezium());
double s  = integrate(f, 0.0, 1.0, 10, simpson());

D

The function f(x) is Func1/Func2 in this program.

module integrate ;
import std.stdio ;
import std.math ;

alias real function(real) realFn ;

class Sigma{
  int n ;
  real a, b , h ;
  realFn fn ; 
  string desc ;
  enum Method {LEFT = 0, RGHT, MIDD, TRAP, SIMP} ;
  static string[5] methodName = 
  ["LeftRect ", "RightRect", "MidRect  ", "Trapezium", "Simpson  "] ;
  static Sigma opCall() { 
    Sigma s = new Sigma() ; 
    return s ; 
  }
  static Sigma opCall(realFn f, int n, real a, real b) { 
    Sigma s = new Sigma(f, n, a, b) ; 
    return s ; 
  } 
  static real opCall(realFn f, int n, real a, real b, Method m) { 
    return Sigma(f,n,a,b).getSum(m) ;
  }
  private this() {} ;
  this(realFn f, int n, real a, real b) {
    setFunction(f) ;
    setStep(n) ;
    setRange(a,b) ;   
  } 
  Sigma opCall(Method m) { 
    return doSum(m) ; 
  }
  Sigma setFunction(realFn f) { 
    this.fn = f ; 
    return this ; 
  }
  Sigma setRange(real a, real b) { 
    this.a = a ; this.b = b ; 
    setInterval() ;
    return this ; 
  }
  Sigma setStep(int n) { 
    this.n = n ; 
    setInterval() ;
    return this ; 
  }
  Sigma setDesc(string d) { 
    this.desc = d ; 
    return this ; 
  }
  private void setInterval() { 
    this.h = (b - a) / n ; 
  }
  private real partSum(int i, Method m) {
    real x = a + h * i ;
    switch(m) {
      case Method.LEFT:
        return fn(x) ;
      case Method.RGHT:
        return fn(x + h) ;
      case Method.MIDD:
        return fn(x + h/2) ;
      case Method.TRAP:
        return (fn(x) + fn(x + h))/2 ;
      default:
    }
    //case SIMPSON:
    return (fn(x) + 4 * fn(x + h/2) + fn(x + h))/6 ;
  } 
  real getSum(Method m) {
    real sum = 0 ;
    for(int i = 0; i < n ; i++)
      sum += partSum(i, m) ;
    return sum * h ;
  } 
  Sigma doSum(Method m) {
    writefln("%10s = %9.6f", methodName[m], getSum(m)) ; 
    return this ; 
  }
  Sigma showSetting() {
    writefln("\n%s\nA = %9.6f, B = %9.6f, N = %s, h = %s", desc, a, b, n, h) ;
    return this ;
  } 
  Sigma doLeft() { return doSum(Method.LEFT) ; }
  Sigma doRight() { return doSum(Method.RGHT) ; }
  Sigma doMid() { return doSum(Method.MIDD) ; }
  Sigma doTrap() { return doSum(Method.TRAP) ; }
  Sigma doSimp() { return doSum(Method.SIMP) ; }
  Sigma doAll() {
    showSetting() ;
    doLeft() ; doRight() ; doMid() ; doTrap() ; doSimp() ;
    return this ; 
  }
}
real Func1(real x) {
  return cos(x) + sin(x)  ;
}
real Func2(real x) {
  return 2.0L/(1 + 4*x*x)  ;
}
void main() {
  // use as a re-usable and manageable object
  Sigma s = Sigma(&Func1, 10, -PI/2, PI/2).showSetting()
  .doLeft().doRight().doMid().doTrap()(Sigma.Method.SIMP) ;     
  s.setFunction(&Func2).setStep(4).setRange(-1.0L,2.0L) ; 
  s.setDesc("Function(x) = 2 / (1 + 4x^2)").doAll() ;
  // use as a single function call 
  writefln("\nLeftRect Integral of FUNC2 =") ;
  writefln("%12.9f (%3dstep)\n%12.9f (%3dstep)\n%12.9f (%3dstep).",
    Sigma(&Func2,   8, -1.0L, 2.0L, Sigma.Method.LEFT),  8,
    Sigma(&Func2,  64, -1.0L, 2.0L, Sigma.Method.LEFT), 64,
    Sigma(&Func2, 512, -1.0L, 2.0L, Sigma.Method.LEFT),512) ;
  writefln("\nSimpson Integral of FUNC2 =") ;
  writefln("%12.9f (%3dstep).",
    Sigma(&Func2, 512, -1.0L, 2.0L, Sigma.Method.SIMP),512) ;
}

Parts of the output:

Function(x) = 2 / (1 + 4x^2)
A = -1.000000, B =  2.000000, N = 4, h = 0.75
 LeftRect  =  2.456897
 RightRect =  2.245132
 MidRect   =  2.496091
 Trapezium =  2.351014
 Simpson   =  2.447732

• One can verifies his result with this applet

Forth

fvariable step

defer method ( fn F: x -- fn[x] )

: left                    execute ;
: right  step f@       f+ execute ;
: mid    step f@ 2e f/ f+ execute ;
: trap
  dup fdup  left
      fswap right f+  2e f/ ;
: simpson
  dup fdup  left
  dup fover mid 4e f* f+
      fswap right f+  6e f/ ;
     
: set-step ( n F: a b -- n F: a )
  fover f- dup 0 d>f f/ step f! ;

: integrate ( xt n F: a b -- F: sigma )
  set-step
  0e
  0 do
    dup fover method f+
    fswap step f@ f+ fswap
  loop
  drop fnip
  step f@ f* ;

\ testing similar to the D example
: test
  ' is method  ' 4 -1e 2e integrate f. ;

: fn1  fsincos f+ ;
: fn2  fdup f* 4e f* 1e f+ 2e fswap f/ ;

7 set-precision
test left    fn2  \ 2.456897
test right   fn2  \ 2.245132
test mid     fn2  \ 2.496091
test trap    fn2  \ 2.351014
test simpson fn2  \ 2.447732

Fortran

In ISO Fortran 95 and later if function f() is not already defined to be "elemental", define an elemental wrapper function around it to allow for array-based initialization:

  ELEMENTAL FUNCTION ELEMF(X)
     REAL :: ELEMF, X
     ELEMF = F(X)
  END FUNCTION ELEMF

Use Array Initializers, Pointers, Array invocation of Elemental functions, Elemental array-array and array-scalar arithmetic, and the SUM intrinsic function

  PROGRAM INTEGRATE
     INTEGER, PARAMETER :: N = 20          ! or whatever
     REAL, PARAMETER :: A = 0.0, B = 15.0  ! or whatever
     REAL, PARAMETER :: H = (B - A) / N
     REAL, PARAMETER, DIMENSION(0:2*N) :: XPOINTS = (/ (A + H*i/(N*2), i = 0, 2*N) /) 
     REAL, DIMENSION(0:2*N), TARGET :: FPOINTS     ! gather up all the f(x) values needed for all methods once and distribute via pointers
     REAL, DIMENSION(:), POINTER :: FLEFT, FMID, FRIGHT
     REAL :: LEFTRECT, MIDRECT, RIGHTRECT, TRAPEZOID, SIMPSON
     
     FPOINTS = ELEMF(XPOINTS)    ! elemental function wrapper for F runs once for each element of XPOINTS
     FLEFT =>  FPOINTS(0 : 2*N-2 : 2)
     FMID =>   FPOINTS(1 : 2*N-1 : 2)
     FRIGHT => FPOINTS(2 : 2*N : 2)
     
     ! Left rectangular rule integral
     LEFTRECT  = H * SUM(FLEFT)
     
     ! Middle rectangular rule integral
     MIDRECT   = H * SUM(FMID)
     
     ! Right rectangular rule integral
     RIGHTRECT = H * SUM(FRIGHT)
     
     ! Trapezoid rule integral
     TRAPEZOID = H / 2 * SUM(FLEFT + FRIGHT)
     
     ! Simpson's rule integral (coming soon)...
     SIMPSON   = H / 6 * SUM(FLEFT + FRIGHT + 4*FMID)
  END PROGRAM INTEGRATE

Haskell

Different approach from most of the other examples: First, the function f might be expensive to calculate, and so it should not be evaluated several times. So, ideally, we want to have positions x and weights w for each method and then just calculate the approximation of the integral by

approx f xs ws = sum [w * f x | (x,w) <- zip xs ws]

Second, let's to generalize all integration methods into one scheme. The methods can all be characterized by the coefficients vs they use in a particular interval. These will be fractions, and for terseness, we extract the denominator as an extra argument v.

Now there are the closed formulas (which include the endpoints) and the open formulas (which exclude them). Let's do the open formulas first, because then the coefficients don't overlap:

integrateOpen :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a
integrateOpen v vs f a b n = approx f xs ws * h / v where
  m = fromIntegral (length vs) * n
  h = (b-a) / fromIntegral m
  ws = concat $ replicate n vs
  c = a + h/2
  xs = [c + h * fromIntegral i | i <- [0..m-1]]

Similarly for the closed formulas, but we need an additional function overlap which sums the coefficients overlapping at the interior interval boundaries:

integrateClosed :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a
integrateClosed v vs f a b n = approx f xs ws * h / v where
  m = fromIntegral (length vs - 1) * n
  h = (b-a) / fromIntegral m
  ws = overlap n vs
  xs = [a + h * fromIntegral i | i <- [0..m]]

overlap :: Num a => Int -> [a] -> [a]
overlap n []  = []
overlap n (x:xs) = x : inter n xs where
  inter 1 ys     = ys
  inter n []     = x : inter (n-1) xs
  inter n [y]    = (x+y) : inter (n-1) xs
  inter n (y:ys) = y : inter n ys

And now we can just define

intLeftRect  = integrateClosed  1 [1,0]     
intRightRect = integrateClosed  1 [0,1]     
intMidRect   = integrateOpen    1 [1]  
intTrapezium = integrateClosed  2 [1,1]   
intSimpson   = integrateClosed  3 [1,4,1] 

or, as easily, some additional schemes:

intMilne     = integrateClosed 45 [14,64,24,64,14]
intOpen1     = integrateOpen    2 [3,3]
intOpen2     = integrateOpen    3 [8,-4,8]

Some examples:

*Main> intLeftRect  (\x -> x*x) 0 1 10  
0.2850000000000001
*Main> intRightRect (\x -> x*x) 0 1 10  
0.38500000000000006
*Main> intMidRect   (\x -> x*x) 0 1 10 
0.3325
*Main> intTrapezium (\x -> x*x) 0 1 10 
0.3350000000000001
*Main> intSimpson   (\x -> x*x) 0 1 10 
0.3333333333333334

Java

The function in this example is assumed to be f(double x).

public class Integrate{
   public static double leftRect(double a, double b, double n){
      double h = (b - a) / n;
      double sum = 0;
      for(double x = a;x <= b - h;x += h)
         sum += f(x);
      return h * sum;
   }

   public static double rightRect(double a, double b, double n){
      double h = (b - a) / n;
      double sum = 0;
      for(double x = a + h;x <= b - h;x += h)
         sum += f(x);
      return h * sum;
   }

   public static double midRect(double a, double b, double n){
      double h = (b - a) / n;
      double sum = 0;
      for(double x = a;x <= b - h;x += h)
         sum += (f(x) + f(x + h));
      return (h / 2) * sum;
   }

   public static double trap(double a, double b, double n){
      double h = (b - a) / n;
      double sum = f(a) + f(b);
      for(int i = 1;i < n;i++)
         sum += f(a + i * h);
      return  h * sum;
   }

   public static double simpson(double a, double b, double n){
      double h = (b - a) / n;
      double sum1 = 0;
      double sum2 = 0;

      for(int i = 0;i < n;i++)
         sum1 += f(a + h * i + h / 2);

      for(int i = 1;i < n;i++)
         sum2 += f(a + h * i);

      return h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2);
   }
   //assume f(double x) elsewhere in the class
}

OCaml

let integrate f a b steps meth =
  let h = (b -. a) /. float_of_int steps in
  let rec helper i s =
    if i >= steps then s
    else helper (succ i) (s +. meth f (a +. h *. float_of_int i) h)
  in
  h *. helper 0 0.

let left_rect f x _ =
  f x

let mid_rect f x h =
  f (x +. h /. 2.)

let right_rect f x h =
  f (x +. h)

let trapezium f x h =
  (f x +. f (x +. h)) /. 2.

let simpson f x h =
  (f x +. 4. *. f (x +. h /. 2.) +. f (x +. h)) /. 6.

let square x = x *. x

let rl = integrate square 0. 1. 10 left_rect
let rm = integrate square 0. 1. 10 mid_rect
let rr = integrate square 0. 1. 10 right_rect
let t = integrate square 0. 1. 10 trapezium
let s = integrate square 0. 1. 10 simpson

Scheme

(define (integrate f a b steps meth)
  (define h (/ (- b a) steps))
  (* h
     (let loop ((i 0) (s 0))
       (if (>= i steps)
           s
           (loop (+ i 1) (+ s (meth f (+ a (* h i)) h)))))))

(define (left-rect f x h) (f x))
(define (mid-rect f x h) (f (+ x (/ h 2))))
(define (right-rect f x h) (f (+ x h)))
(define (trapezium f x h) (/ (+ (f x) (f (+ x h))) 2))
(define (simpson f x h) (/ (+ (f x) (* 4 (f (+ x (/ h 2)))) (f (+ x h))) 6))

(define (square x) (* x x))

(define rl (integrate square 0 1 10 left-rect))
(define rm (integrate square 0 1 10 mid-rect))
(define rr (integrate square 0 1 10 right-rect))
(define t (integrate square 0 1 10 trapezium))
(define s (integrate square 0 1 10 simpson))