Numerical Integration

Essential Algorithms: A Practical Approach to Computer Algorithms
The most straightforward uses Newton-Cotes formulas, which use a series of polynomials to approximate the function. The two most basic kinds of Newton-Cotes formulas are the rectangle rule and the trapezoid rule.

http://introcs.cs.princeton.edu/java/93integration/

Midpoint rule. Goal: given continuous function f(x) of one variable, compute ∫ f(x) dx over interval from a to b. Estimate integral by M = (b-a) * f(c), where c = (a + b)/2.

Trapezoidal rule. Goal: given continuous function f(x) of one variable, compute ∫ f(x) dx over interval from a to b.TrapezoidalRule.java numerically integrates a function of one variable using the trapezoidal rule. We can estimate the integral of f(x) from a to b using the formula T = (b-a)/2 (f(a) + f(b)). Breaking the interval from a to b up into N equally spaced intervals (and combining common terms) we obtain the formula:

where the interval [a, b] is broken up into N subintervals of uniform size h = (b - a) / N. Under certain technical conditions, if N is large then the formula above is a good estimate of the integral.

Simpson's rule. The trapezoidal rule is rarely used to integrate in practice. For smooth f, the midpoint rule is approximately twice as accurate as the trapezoidal rule, and the errors have different signs. By combining the two expressions, we obtain a more accurate estimate of f: S = 2/3*M + 1/3*T. This combination is known as Simpson's 1/3 rule. S = (b-a)/6 (f(a) + 4f(c) + f(b)), where c = (a + b)/2. Breaking the interval from a to b up into N equally spaced intervals (and combining common terms) we obtain the formula:

where the interval [a, b] is broken up into N subintervals of uniform size h = (b - a) / (N - 1) and the 2/3 and 4/3 coefficients alternate throughout the interior. Here are some nice animations of numerical quadrature. Under certain technical conditions, if N is large then the formula above is a good estimate of the integral. The program SimpsonsRule.java numerically integrates x^4 log (x + sqrt(x^2 + 1)) from a = 0 to b = 2.

http://introcs.cs.princeton.edu/java/93integration/TrapezoidalRule.java.html

public class TrapezoidalRule {

  /**********************************************************************
   * Standard normal distribution density function.
   * Replace with any sufficiently smooth function.
   **********************************************************************/
   static double f(double x) {
      return Math.exp(- x * x / 2) / Math.sqrt(2 * Math.PI);
   }

  /**********************************************************************
   * Integrate f from a to b using the trapezoidal rule.
   * Increase N for more precision.
   **********************************************************************/
   static double integrate(double a, double b, int N) {
      double h = (b - a) / N;              // step size
      double sum = 0.5 * (f(a) + f(b));    // area
      for (int i = 1; i < N; i++) {
         double x = a + h * i;
         sum = sum + f(x);
      }

      return sum * h;
   }



   // sample client program
   public static void main(String[] args) { 
      double a = Double.parseDouble(args[0]);
      double b = Double.parseDouble(args[1]);
      System.out.println(integrate(a, b, 1000));
   }

}

Essential Algorithms: A Practical Approach to Computer Algorithms

Float: UseRectangleRule(Float: function(), Float: xmin, Float: xmax, Integer: num_intervals) // Calculate the width of a rectangle. Float: dx = (xmax - xmin) / num_intervals // Add up the rectangles' areas. Float: total_area = 0 Float: x = xmin For i = 1 To num_intervals total_area = total_area + dx * function(x) x = x + dx Next i Return total_area End UseRectangleRule

Float: UseTrapezoidRule(Float: function(), Float: xmin, Float: xmax, Integer: num_intervals) // Calculate the width of a trapezoid. Float: dx = (xmax - xmin) / num_intervals // Add up the trapezoids' areas. Float: total_area = 0 Float: x = xmin For i = 1 To num_intervals total_area = total_area + dx * (function(x) + function(x + dx)) / 2 x = x + dx Next i Return total_area End UseTrapezoidRule

ADAPTIVE QUADRATURE
When calculating the area of a slice, this program first uses a single trapezoid to approximate its area. It then breaks the slice into two pieces and uses two smaller trapezoids to calculate their areas. If the difference between the larger trapezoid's area and the sum of the areas of the smaller trapezoids is more than a certain percentage, the program divides the slice into two pieces and calculates the areas of the pieces in the same way.
// Integrate by using an adaptive midpoint trapezoid rule. Float: IntegrateAdaptiveMidpoint(Float: function(), Float: xmin, Float: xmax, Integer: num_intervals, Float: max_slice_error) // Calculate the width of the initial trapezoids. Float: dx = (xmax - xmin) / num_intervals double total = 0 // Add up the trapezoids' areas. Float: total_area = 0 Float: x = xmin For i = 1 To num_intervals // Add this slice's area. total_area = total_area + SliceArea(function, x, x + dx, max_slice_error) x = x + dx Next i Return total_area End IntegrateAdaptiveMidpoint // Return the area for this slice. Float: SliceArea(Float: function(),Float: x1, Float: x2, Float: max_slice_error) // Calculate the function at the endpoints and the midpoint. Float: y1 = function(x1) Float: y2 = function(x2) Float: xm = (x1 + x2) / 2 Float: ym = function(xm) // Calculate the area for the large slice and two subslices. Float: area12 = (x2 - x1) * (y1 + y2) / 2.0 Float: area1m = (xm - x1) * (y1 + ym) / 2.0 Float: aream2 = (x2 - xm) * (ym + y2) / 2.0 Float: area1m2 = area1m + aream2 // See how close we are. Float: error = (area1m2 - area12) / area12 // See if this is small enough. If (Abs(error) < max_slice_error) Then Return area1m2 // The error is too big. Divide the slice and try again. Return SliceArea(function, x1, xm, max_slice_error) + SliceArea(function, xm, x2, max_slice_error) End SliceArea
http://introcs.cs.princeton.edu/java/93integration/AdaptiveQuadrature.java.html
http://rosettacode.org/wiki/Numerical_integration#Java