Numerical integration

From Rosetta Code
Revision as of 18:40, 21 December 2007 by rosettacode>Mwn3d (Created page. Added Java example and explanation of Simpson's rule from my programming class.)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Task
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)

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 = a + b;
      for(int i = 1;i <= n-1;i++)
         sum += 2 * (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 - 1;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);
   }
   ...
}