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 the wiki) 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;

BASIC

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 * (b - a) / n))
     NEXT i
     trap = (b - a) / (2 * n) * 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

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 += h * (f(x));
      return 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 += h * (f(x));
      return 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 += (h / 2) * (f(x) + f(x + h));
      return 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-1;i++)
         sum += 2 * f((a + i * (b - a) / n));
      return (b - a) / (2 * n) * 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 - 1;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
}