# AdaptiveQuadrature.java

Below is the syntax highlighted version of AdaptiveQuadrature.java from §9.3 Symbolic Methods.

```
/******************************************************************************
*  Compilation:  javac AdaptiveQuadrature.java
*  Execution:    java AdaptiveQuadrature a b
*
*  Numerically integrate the function in the interval [a, b] using
*  the extrapolated Simpson's rule in an adaptive recursive algorithm.
*
*  See https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/quad.pdf
*
*  % java AdaptiveQuadrature -3 3          // true answer = 0.9973002040...
*  0.9973002039366737
*  0.9973002036280225
*
*  % java AdaptiveQuadrature 0 10000000    // true answer = 1/2
*  1.9947114020071635
*  0.500000121928651
*
*
*  Note that even with 1 million subintervals, the trapezoid rule
*  gets the integral of f(x) = exp(-x^2) / sqrt(2 pi) from 0 to 10,000,000
*  completely wrong.
*  In contrast, the adaptive quadrature rule is accurate to 8 significant
*  digits. In this example, both methods evaluate f() at roughly 1 million
*  points.
*
*  Note: many function calls can be saved by computing f(a), f(b), and f(c)
*        only once per recursive call (instead of twice) and adding extra
*        arguments to adapative() and passing the already computed values
*        that will get reused in the next recursive call.
*
******************************************************************************/

public class AdaptiveQuadrature {
private final static double EPSILON = 1E-6;

/***************************************************************************
* 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);
}

// trapezoid rule
static double trapezoid(double a, double b, int n) {
double h = (b - a) / n;
double sum = 0.5 *  h * (f(a) + f(b));
for (int k = 1; k < n; k++)
sum = sum + h * f(a + h*k);
return sum;
}

// adaptive quadrature
public static double adaptive(double a, double b) {
double h = b - a;
double c = (a + b) / 2.0;
double d = (a + c) / 2.0;
double e = (b + c) / 2.0;
double Q1 = h/6  * (f(a) + 4*f(c) + f(b));
double Q2 = h/12 * (f(a) + 4*f(d) + 2*f(c) + 4*f(e) + f(b));
if (Math.abs(Q2 - Q1) <= EPSILON)
return Q2 + (Q2 - Q1) / 15;
else
return adaptive(a, c) + adaptive(c, b);
}

// sample client program
public static void main(String[] args) {
double a = Double.parseDouble(args);
double b = Double.parseDouble(args);
StdOut.println(trapezoid(a, b, 1000000));
StdOut.println(adaptive(a, b));
}

}
```

Copyright © 2000–2017, Robert Sedgewick and Kevin Wayne.
Last updated: Fri Aug 31 21:08:05 EDT 2018.