Numerical integration

Revision as of 21:34, 21 December 2007 by rosettacode>Mwn3d (→‎{{header|Java}}: Fixed trapezoid method again...)

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).

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

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)

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
}