Numerical Integration

Definition of Definite Integral

Simply put, the definite integral \(\int_{l}^{r}f(x)\mathrm{d}x\) of a function \(f(x)\) over the interval \([l,r]\) refers to the area of the region bounded by \(f(x)\) and the x-axis over \([l,r]\) (where the part above the x-axis is positive, and the part below is negative).

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

Simpson's Method

The idea of this method is to divide the integration interval into several small segments, and apply the integral formula for quadratic functions to each segment.

Quadratic Function Integral Formula (Simpson's Rule)

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. Then

\[ \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 this Simpson's rule, we first introduce a common Simpson's integration method.

Common 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 subintervals \(x\).

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

We can then compute the integral value for each small interval \([x_ {2i-2}, x_ {2i}]\), \(i = 1\ldots n\), and sum the integrals of all intervals to obtain the total integral.

For an interval \([x_ {2i-2}, x_ {2i}]\), \(i = 1\ldots n\), selecting three points \((x_ {2i-2}, x_ {2i-1}, x_ {2i})\) uniquely defines a parabola, giving a function \(P(x)\). Computing the integral of the original function over this interval becomes computing the integral of the new quadratic function \(P(x)\) over this segment. Thus, we can use Simpson's rule to approximate it.

\(\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 over all segments yields 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 directly state the conclusion: the error of the common Simpson's method is:

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

where \(\xi\) is some value in the interval \([l,r]\).

Implementation

C++Python

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;
}

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 common method's accuracy is inevitably limited by \(n\) in terms of time. We should find a more suitable approach.

The only remaining problem is how to perform segmentation. Too few segments lead to large computational errors, while too many segments reduce time efficiency. We need to find a balance between accuracy and efficiency.

We consider this: if a segment of the graph is already very close to a quadratic function, directly applying the formula to compute the integral yields high precision, and there is no need to further subdivide this segment.

Thus, we have the following segmentation method: each time, judge how closely the current segment resembles a quadratic function. If it is sufficiently similar, directly substitute into the formula for computation; otherwise, split the current segment into left and right halves and solve recursively.

Now only one question remains: how to judge whether each segment is similar to a quadratic function?

We directly substitute the current segment into the formula to compute its integral. Then we split the current segment at its midpoint into two segments and directly substitute these two segments into the formula to compute their integrals. If the integral of the current segment and the sum of the integrals of the two split segments differ very little, we can consider the current segment sufficiently similar to a quadratic function and stop recursive subdivision.

The above is the idea of the adaptive Simpson's method. When performing the divide-and-conquer judgment, besides checking if the precision is correct, a minimum number of iterations is generally enforced.

Reference code is as follows:

C++Python

double simpson(double l, double r) {
  double mid = (l + r) / 2;
  return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;  // Simpson's rule
}

double asr(double l, double r, double eps, double ans, int step) {
  double mid = (l + r) / 2;
  double fl = simpson(l, mid), fr = simpson(mid, r);
  if (abs(fl + fr - ans) <= 15 * eps && step < 0)
    return fl + fr + (fl + fr - ans) / 15;  // If sufficiently similar, return directly
  return asr(l, mid, eps / 2, fl, step - 1) +
         asr(mid, r, eps / 2, fr, step - 1);  // Otherwise, split into two segments and solve recursively
}

double calc(double l, double r, double eps) {
  return asr(l, r, eps, simpson(l, r), 12);
}

def simpson(l, r):
    mid = (l + r) / 2
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6  # Simpson's rule


def asr(l, r, eps, ans, step):
    mid = (l + r) / 2
    fl = simpson(l, mid)
    fr = simpson(mid, r)
    if abs(fl + fr - ans) <= 15 * eps and step < 0:
        return fl + fr + (fl + fr - ans) / 15  # If sufficiently similar, return directly
    return asr(l, mid, eps / 2, fl, step - 1) + asr(
        mid, r, eps / 2, fr, step - 1
    )  # Otherwise, split into two segments and solve recursively


def calc(l, r, eps):
    return asr(l, r, eps, simpson(l, r), 12)

Exercises

References

https://doi.org/10.1145/321526.321537: This article discusses improvements to the adaptive Simpson's method, detailing the origin and advantages of the constant 15 in the code above.

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用