adaptive Simpson integration over an 1D interval
| 63 | |
| 64 | /// adaptive Simpson integration over an 1D interval |
| 65 | inline double adaptiveSimpson(const std::function<double (double)> &f, double x0, double x1, double eps = 1e-6, int depth = 6) { |
| 66 | int count = 0; |
| 67 | /* Define an recursive lambda function for integration over subintervals */ |
| 68 | std::function<double (double, double, double, double, double, double, double, double, int)> integrate = |
| 69 | [&](double a, double b, double c, double fa, double fb, double fc, double I, double eps, int depth) { |
| 70 | /* Evaluate the function at two intermediate points */ |
| 71 | double d = 0.5 * (a + b), e = 0.5 * (b + c), fd = f(d), fe = f(e); |
| 72 | |
| 73 | /* Simpson integration over each subinterval */ |
| 74 | double h = c-a, |
| 75 | I0 = (1.0/12.0) * h * (fa + 4.0*fd + fb), |
| 76 | I1 = (1.0/12.0) * h * (fb + 4.0*fe + fc), |
| 77 | Ip = I0+I1; |
| 78 | ++count; |
| 79 | |
| 80 | /* Stopping criterion from J.N. Lyness (1969) |
| 81 | "Notes on the adaptive Simpson quadrature routine" */ |
| 82 | if (depth <= 0 || std::abs(Ip-I) < 15.0*eps) { |
| 83 | // Richardson extrapolation |
| 84 | return Ip + (1.0/15.0) * (Ip-I); |
| 85 | } |
| 86 | |
| 87 | return integrate(a, d, b, fa, fd, fb, I0, 0.5*eps, depth-1) + |
| 88 | integrate(b, e, c, fb, fe, fc, I1, 0.5*eps, depth-1); |
| 89 | }; |
| 90 | double a = x0, b = 0.5 * (x0+x1), c = x1; |
| 91 | double fa = f(a), fb = f(b), fc = f(c); |
| 92 | double I = (c-a) * (1.0/6.0) * (fa+4.0*fb+fc); |
| 93 | return integrate(a, b, c, fa, fb, fc, I, eps, depth); |
| 94 | } |
| 95 | |
| 96 | /// Nested adaptive Simpson integration over a 2D rectangle |
| 97 | inline double adaptiveSimpson2D(const std::function<double (double, double)> &f, double x0, double y0, |
no test coverage detected