AIOJ

Definition of the definite integral

For a function $f(x)$ on the interval $[l,r]$ , its definite integral $\int_{l}^{r}f(x)\mathrm{d}x$ is the area of the region enclosed by $f(x)$ over $[l,r]$ and the $x$ -axis (the part above the $x$ -axis is positive, while the part below the $x$ -axis is negative).

In many situations, we need to compute an approximate value of an integral efficiently and accurately. The Simpson's method introduced below is one such numerical integration method.

Simpson's method

The idea is to divide the interval of integration into several smaller intervals and apply the integral formula for a quadratic function to each one.

Integral formula for a quadratic function (Simpson's formula)

For a quadratic function $f(x)=ax^2+bx+c$ , we have:

\int_l^r f(x) {\mathrm d}x = \frac{(r-l)(f(l)+f(r)+4 f(\frac{l+r}{2}))}{6}

Derivation: For a quadratic function $f(x)=ax^2+bx+c$ ; Integrating gives $F(x)=\int_0^x f(x) {\mathrm d}x = \frac{a}{3}x^3+\frac{b}{2}x^2+cx+D$ where D is a constant, so

\begin{aligned} \int_l^r f(x) {\mathrm d}x &= F(r)-F(l) \\ &= \frac{a}{3}(r^3-l^3)+\frac{b}{2}(r^2-l^2)+c(r-l) \\ &=(r-l)(\frac{a}{3}(l^2+r^2+lr)+\frac{b}{2}(l+r)+c) \\ &=\frac{r-l}{6}(2al^2+2ar^2+2alr+3bl+3br+6c)\\ &=\frac{r-l}{6}((al^2+bl+c)+(ar^2+br+c)+4(a(\frac{l+r}{2})^2+b(\frac{l+r}{2})+c)) \\ &=\frac{r-l}{6}(f(l)+f(r)+4f(\frac{l+r}{2})) \end{aligned}

Based on Simpson's formula, we first introduce the ordinary Simpson integration method.

Ordinary Simpson's method

This method was published in a paper by Thomas Simpson in 1743.

Description

Given a natural number $n$ , divide the interval $[l, r]$ into $2n$ equal-length intervals $x$ .

$x_i = l + i h, ~~ i = 0 \ldots 2n,$ $h = \frac {r-l} {2n}.$

We can compute the integral over each small interval $[x_ {2i-2}, x_ {2i}]$ , $i = 1\ldots n$ , and sum the values over all intervals to obtain the total integral.

For an interval $[x_ {2i-2}, x_ {2i}]$ , $i = 1\ldots n$ , choose its three points $(x_ {2i-2}, x_ {2i-1}, x_ {2i})$ . They determine a parabola and hence a function $P(x)$ , which exists and is unique. Computing the integral of the original function over this interval is then reduced to computing the integral of the new quadratic function $P(x)$ over the same interval. We can approximate it using Simpson's formula.

$\int_{x_ {2i-2}} ^ {x_ {2i}} f (x) ~dx \approx \int_{x_ {2i-2}} ^ {x_ {2i}} P (x) ~dx = \left(f(x_{2i-2}) + 4f(x_{2i-1})+(f(x_{2i})\right)\frac {h} {3}$

Summing the intervals gives the following conclusion:

$\int_l ^ r f (x) dx \approx \left(f (x_0) + 4 f (x_1) + 2 f (x_2) + 4f(x_3) + 2 f(x_4) + \ldots + 4 f(x_{2N-1}) + f(x_{2N}) \right)\frac {h} {3}$

Error

We state directly that the error of the ordinary Simpson's method is:

-\tfrac{1}{90} \left(\tfrac{r-l}{2}\right)^5 f^{(4)}(\xi)

Here $\xi$ is some value in the interval $[l,r]$ .

Implementation

C++

constexpr int N = 1000 * 1000;

double simpson_integration(double a, double b) {
  double h = (b - a) / N;
  double s = f(a) + f(b);
  for (int i = 1; i <= N - 1; ++i) {
    double x = a + h * i;
    s += f(x) * ((i & 1) ? 4 : 2);
  }
  s *= h / 3;
  return s;
}

Python

N = 1000 * 1000

def simpson_integration(a, b):
    h = (b - a) / N
    s = f(a) + f(b)
    for i in range(1, N):
        x = a + h * i
        if i & 1:
            s = s + f(x) * 4
        else:
            s = s + f(x) * 2
    s = s * (h / 3)
    return s

Adaptive Simpson's method

The ordinary method is inevitably limited by $n$ when accuracy must be guaranteed, so we should look for a more suitable method.

The only remaining question is how to divide the interval. Too few segments lead to a large error, while too many reduce time efficiency. We need to find a balance between accuracy and efficiency.

Consider a portion of the graph that is already very close to a quadratic function. Substituting it directly into the formula gives a highly accurate integral, so there is no need to subdivide that portion further.

This gives the following subdivision method: each time, judge how similar the current segment is to a quadratic function. If it is sufficiently similar, substitute it directly into the formula; otherwise, split it into left and right segments and solve them recursively.

There is only one problem left: how do we determine whether each segment is similar to a quadratic function?

Substitute the current segment directly into the formula to obtain its integral, then split the segment at its midpoint and apply the formula directly to the two resulting segments. If the integral over the current segment differs only slightly from the sum of the two split integrals, we can regard the current segment as very similar to a quadratic function and stop the recursion.

This is the idea behind adaptive Simpson's method. When making divide-and-conquer decisions, in addition to checking the required accuracy, one generally also forces a minimum number of iterations.

Reference code follows: