yuygfgg的小站

Simpson's Rule is based on approximating the curve with a parabola (i.e., a quadratic polynomial), instead of a straight line.

Mathematical Proof

Simpson's Rule approximates the integrand using a quadratic polynomial. Suppose we have three points on the interval ([a, b]): (a), (m) (the midpoint), and (b). We approximate the function values at these three points with a quadratic polynomial (P(x)):

$P(x) = Ax^2 + Bx + C$

Through interpolation, we can ensure that (P(a) = f(a)), (P(m) = f(m)), and (P(b) = f(b)). This quadratic polynomial captures the function's changes over the interval better than a straight line.

Interpolation Calculation

First, we write three equations:

Starting from the first equation:

$Aa^2 + Ba + C = f(a)$

Setting (x = a), we get:

$Aa^2 + Ba + C = f(a)$

Similarly, setting (x = b), we have:

$Ab^2 + Bb + C = f(b)$

Setting (x = \frac{a+b}{2}), we obtain:

$A\left(\frac{a+b}{2}\right)^2 + B\left(\frac{a+b}{2}\right) + C = f\left(\frac{a+b}{2}\right)$

Now, we have three equations:

$Aa^2 + Ba + C = f(a) \quad \text{(1)}$

$Ab^2 + Bb + C = f(b) \quad \text{(2)}$

$A\left(\frac{a+b}{2}\right)^2 + B\left(\frac{a+b}{2}\right) + C = f\left(\frac{a+b}{2}\right) \quad \text{(3)}$

Eliminating (C) from equations (1) and (2):

$Aa^2 + Ba = f(a) - C \quad \text{(4)}$

$Ab^2 + Bb = f(b) - C \quad \text{(5)}$

Subtracting (5) from (4), we have:

$A(a^2 - b^2) + B(a - b) = f(a) - f(b)$

Since (a^2 - b^2 = (a+b)(a-b)):

$A(a+b)(a-b) + B(a-b) = f(a) - f(b)$

Factoring out (a-b):

$\left(A(a+b) + B\right)(a-b) = f(a) - f(b)$

Thus:

$A(a+b) + B = \frac{f(a) - f(b)}{a-b} \quad \text{(6)}$

Substituting (6) into (3):

$A\left(\frac{a+b}{2}\right)^2 + \left(\frac{f(a) - f(b)}{a-b}\right)\left(\frac{a+b}{2}\right) + C = f\left(\frac{a+b}{2}\right)$

Simplifying:

$A\frac{a^2 + 2ab + b^2}{4} + \frac{(f(a) - f(b))(a+b)}{2(a-b)} + C = f\left(\frac{a+b}{2}\right)$

$4A(a^2 + 2ab + b^2) + 8\frac{(f(a) - f(b))(a+b)}{a-b} + 4C = 4f\left(\frac{a+b}{2}\right)$

$4A(a^2 + 2ab + b^2) = 4f\left(\frac{a+b}{2}\right) - 8\frac{(f(a) - f(b))(a+b)}{a-b} - 4C$

$A = \frac{4f\left(\frac{a+b}{2}\right) - 8\frac{(f(a) - f(b))(a+b)}{a-b} - 4C}{4(a^2 + 2ab + b^2)}$

$B = \frac{f(a) - f(b)}{a-b} - A(a+b)$

$C = f(a) - Aa^2 - Ba$

Thus, (A), (B), and (C) can be expressed using (a), (b), (C), and (f(x)) as follows:

$A = \frac{4f\left(\frac{a+b}{2}\right) - 8\frac{(f(a) - f(b))(a+b)}{a-b} - 4C}{4(a^2 + 2ab + b^2)}$

$B = \frac{f(a) - f(b)}{a-b} - A(a+b)$

$C = f(a) - Aa^2 - Ba$

Thus, we can write out the specific expression for the interpolation function (P(x)):

Substituting (A), (B), and (C) completely into the equation (Ax^2 + Bx + C), we get:

$\left(\frac{2f(a) + 2f(b) - 4f\left(\frac{a+b}{2}\right)}{(a-b)^2}\right)x^2 + \left(\frac{-af(a) - 3af(b) + 4a f\left(\frac{a+b}{2}\right) - 3bf(a) - bf(b) + 4b f\left(\frac{a+b}{2}\right)}{(a-b)^2}\right)x + \frac{a^2 f(b) + ab f(a) + ab f(b) - 4ab f\left(\frac{a+b}{2}\right) + b^2 f(a)}{(a-b)^2}$

Integration

We integrate the expression we just obtained.

Here we calculate the definite integral of the function from (a) to (b).

from sympy import symbols, Function, integrate, simplify

# Define symbols
x, a, b, c = symbols('x a b c')
f = Function('f')

# Define coefficients A, B, C
A = 2*f(a)/(a**2 - 2*a*b + b**2) + 2*f(b)/(a**2 - 2*a*b + b**2) - 4*f(a/2 + b/2)/(a**2 - 2*a*b + b**2)
B = -a*f(a)/(a**2 - 2*a*b + b**2) - 3*a*f(b)/(a**2 - 2*a*b + b**2) + 4*a*f(a/2 + b/2)/(a**2 - 2*a*b + b**2) - 3*b*f(a)/(a**2 - 2*a*b + b**2) - b*f(b)/(a**2 - 2*a*b + b**2) + 4*b*f(a/2 + b/2)/(a**2 - 2*a*b + b**2)
C = a**2*f(b)/(a**2 - 2*a*b + b**2) + a*b*f(a)/(a**2 - 2*a*b + b**2) + a*b*f(b)/(a**2 - 2*a*b + b**2) - 4*a*b*f(a/2 + b/2)/(a**2 - 2*a*b + b**2) + b**2*f(a)/(a**2 - 2*a*b + b**2)

# Define the equation Ax^2 + Bx + C
expr = A*x**2 + B*x + C

# Calculate the definite integral
integral_result = integrate(expr, (x, a, b))

simp = simplify(integral_result)

# Print the result
print("Definite integral result:")
print(simp)

The output is as follows:

-a*f(a)/6 - a*f(b)/6 - 2*a*f(a/2 + b/2)/3 + b*f(a)/6 + b*f(b)/6 + 2*b*f(a/2 + b/2)/3

Which simplifies to:

$-\frac{a f(a)}{6} - \frac{a f(b)}{6} - \frac{2a f\left(\frac{a}{2} + \frac{b}{2}\right)}{3} + \frac{b f(a)}{6} + \frac{b f(b)}{6} + \frac{2b f\left(\frac{a}{2} + \frac{b}{2}\right)}{3}$

Continuing to simplify:

We define (h = \frac{b-a}{2}), (m = \frac{a+b}{2})

$-\frac{af(a)}{6} - \frac{af(b)}{6} - \frac{2af\left(\frac{a}{2} + \frac{b}{2}\right)}{3} + \frac{bf(a)}{6} + \frac{bf(b)}{6} + \frac{2bf\left(\frac{a}{2} + \frac{b}{2}\right)}{3}$

We already know:

$\frac{a}{2} + \frac{b}{2} = m$

So:

$-\frac{2af(m)}{3} + \frac{2bf(m)}{3}$

Can be simplified to:

$\frac{2(b-a)f(m)}{3} = \frac{4h f(m)}{3}$

Now, the remaining part is:

$-\frac{af(a)}{6} - \frac{af(b)}{6} + \frac{bf(a)}{6} + \frac{bf(b)}{6}$

This can be rearranged as:

$\frac{-a f(a) + b f(a) - a f(b) + b f(b)}{6}$

Substituting (a = m - h) and (b = m + h):

$\frac{-(m-h) f(a) + (m+h) f(a) - (m-h) f(b) + (m+h) f(b)}{6}$

$\frac{2h f(a) + 2h f(b)}{6} = \frac{h f(a) + h f(b)}{3}$

Combining these two parts:

$\frac{h f(a) + h f(b)}{3} + \frac{4h f(m)}{3} = \frac{h(f(a) + 4f(m) + f(b))}{3}$

This is the expression we need:

$\frac{h}{3} [f(a) + 4f(m) + f(b)]$

This is the final expression for Simpson's Rule.

By taking more points (a, b), we can accurately approximate any given function with a quadratic function, thereby using numerical methods to calculate the definite integral of any function.

Code Verification

In the code, we calculate $4f(m)$ separately and calculate $f(a)$ and $f(b)$ separately. Since all $f(a)$ and $f(b)$ values, except for the first and last, act as the end of one interval and the start of another, they need to be multiplied by 2.

To run the code properly, we need to install the following dependencies: libomp for parallel computing and clang++ for compiling.

#include <iostream>
#include <cmath>
#include <stdexcept>
#include <functional>
#include <omp.h>

// Define the function f(x)
long double f(long double x) {
    return std::tan(x) + 114514 * std::pow(x, 1.0235) + 0.6 * std::sin(1 / std::pow(2 * x + 1, std::tan(std::pow(x, 0.14))));
}

// Use Simpson's Rule to approximate the definite integral
long double simpsons_rule(std::function<long double(long double)> f, long double a, long double b, long long n) {
    if (n % 2 != 0) {
        throw std::invalid_argument("n must be an even integer.");
    }

    long double h = (b - a) / n;
    long double total = f(a) + f(b);
    long double sum1 = 0.0;
    long double sum2 = 0.0;

    #pragma omp parallel for reduction(+:sum1)
    for (long long i = 1; i < n; i += 2) {
        sum1 += 4 * f(a + i * h);
    }

    #pragma omp parallel for reduction(+:sum2)
    for (long long i = 2; i < n; i += 2) {
        sum2 += 2 * f(a + i * h);
    }

    total += sum1 + sum2;

    return total * (h / 3);
}

int main() {
    try {
        long double a = 0.1; // Lower limit of integration
        long double b = 1;   // Upper limit of integration
        long long n = 10000000000LL; // Number of subintervals (must be an even number)
        
        long double simpsons_result = simpsons_rule(f, a, b, n);
        std::cout << "Simpson's Rule Result: " << simpsons_result << std::endl;
    } catch (const std::exception &e) {
        std::cerr << e.what() << std::endl;
    }
}

We can compile and run with the following command:

clang++ -Ofast -march=native -Xpreprocessor -fopenmp -I$(brew --prefix libomp)/include -L$(brew --prefix libomp)/lib -lomp  -mfpu=neon -mfloat-abi=hard -o integrate ./*.cpp -lm -std=c++20 -funroll-loops && ./integrate