如何实现可高效计算数值原函数的`NumericalAntiderivative`类？

Let's build this class step by step—we'll focus on precomputing integral data during initialization for fast lookups later, which aligns perfectly with your usage example. Here's a practical, robust implementation approach:

Core Design Idea #

The key is to precompute cumulative integral values across the interval [minXofInterest, maxXofInterest] when the class is initialized. We'll sample the function at regular intervals, calculate the integral from minXofInterest to each sample point, and store these values. When you call ad(x) later, we'll use interpolation to quickly estimate the antiderivative at x.

Full Implementation Code #

#include <vector>
#include <algorithm>
#include <stdexcept>
#include <cmath>
#include <functional>

class NumericalAntiderivative {
private:
    std::function<double(double)> func;
    double min_x;
    double max_x;
    size_t num_samples;
    std::vector<double> x_samples;
    std::vector<double> cumulative_integrals;

    // Helper: Compute cumulative integrals using Simpson's Rule (higher precision than trapezoidal)
    void compute_cumulative_integrals() {
        cumulative_integrals.resize(num_samples, 0.0);
        const double step = (max_x - min_x) / (num_samples - 1);

        // Simpson's Rule requires an even number of intervals (odd number of samples)
        if (num_samples % 2 == 0) {
            num_samples++;
            x_samples.resize(num_samples);
            for (size_t i = 0; i < num_samples; ++i) {
                x_samples[i] = min_x + i * step;
            }
        }

        // Base case: integral from min_x to min_x is 0
        cumulative_integrals[0] = 0.0;

        // Calculate integrals in pairs of intervals (Simpson's Rule requirement)
        for (size_t i = 1; i < num_samples; i += 2) {
            const double x0 = x_samples[i-1];
            const double x1 = x_samples[i];
            const double x2 = x_samples[i+1];
            const double f0 = func(x0);
            const double f1 = func(x1);
            const double f2 = func(x2);

            // Simpson's Rule formula for [x0, x2]
            const double interval_integral = (x2 - x0) / 6.0 * (f0 + 4*f1 + f2);
            cumulative_integrals[i+1] = cumulative_integrals[i-1] + interval_integral;

            // Fill in the middle point (even index) with linear interpolation
            if (i+1 < num_samples - 1) {
                cumulative_integrals[i] = cumulative_integrals[i-1] + (interval_integral / 2.0);
            }
        }
    }

public:
    // Constructor: takes target function, interval bounds, and optional sample count
    NumericalAntiderivative(std::function<double(double)> f, double minX, double maxX, size_t samples = 10001) 
        : func(std::move(f)), min_x(minX), max_x(maxX), num_samples(samples) {
        if (min_x >= max_x) {
            throw std::invalid_argument("minX must be strictly less than maxX");
        }
        if (num_samples < 3) {
            throw std::invalid_argument("At least 3 samples are required for Simpson's Rule");
        }

        // Initialize evenly spaced sample points
        x_samples.reserve(num_samples);
        const double step = (max_x - min_x) / (num_samples - 1);
        for (size_t i = 0; i < num_samples; ++i) {
            x_samples.push_back(min_x + i * step);
        }

        // Run precomputation (the "slow initialization" step)
        compute_cumulative_integrals();
    }

    // Operator to get antiderivative value at a given x
    double operator()(double x) const {
        // Handle out-of-bounds with linear extrapolation (adjust to throw if preferred)
        if (x <= min_x) {
            const double slope = (cumulative_integrals[1] - cumulative_integrals[0]) / (x_samples[1] - x_samples[0]);
            return cumulative_integrals[0] + slope * (x - min_x);
        }
        if (x >= max_x) {
            const size_t last = num_samples - 1;
            const size_t second_last = last - 1;
            const double slope = (cumulative_integrals[last] - cumulative_integrals[second_last]) / (x_samples[last] - x_samples[second_last]);
            return cumulative_integrals[last] + slope * (x - max_x);
        }

        // Find the interval containing x using binary search
        auto it = std::upper_bound(x_samples.begin(), x_samples.end(), x);
        const size_t idx = std::distance(x_samples.begin(), it) - 1;

        // Linear interpolation between adjacent precomputed points
        const double x0 = x_samples[idx];
        const double x1 = x_samples[idx+1];
        const double y0 = cumulative_integrals[idx];
        const double y1 = cumulative_integrals[idx+1];

        return y0 + (y1 - y0) * (x - x0) / (x1 - x0);
    }
};

Key Details Explained #

• Precomputation: The constructor uses Simpson's Rule (more accurate than the trapezoidal rule for smooth functions) to calculate cumulative integrals. This is the slow initialization step you noted, but it only runs once.
• Fast Lookups: When you call ad(x), we use binary search to find the relevant interval and linear interpolation to get the antiderivative value—this is nearly instantaneous.
• Out-of-Bounds Handling: The code uses linear extrapolation for values outside your target interval. If you prefer strict bounds checking, replace this with a throw statement.
• Sample Count: The default 10001 samples balance precision and memory usage. Tweak this number based on how accurate you need the results to be.

Usage Example (Matching Your Code) #

#include <iostream>
#include <cmath>

int main() {
    auto f = [](double x){ return std::sin(x); };
    try {
        NumericalAntiderivative ad(f, 0.0, 10.0); // Initialize over [0, 10]
        std::cout << "Integral of f from 0 to 5 is " << (ad(5.0) - ad(0.0)) << std::endl;
        // Expected value: ~1 - cos(5) ≈ 1.28366
    } catch (const std::exception& e) {
        std::cerr << "Error: " << e.what() << std::endl;
        return 1;
    }
    return 0;
}

Notes for Improvement #

• If your function has discontinuities, add breakpoints in the interval during initialization to sample around those points.
• For very large intervals, use adaptive sampling (more samples in regions where the function changes rapidly) instead of uniform sampling.
• For higher precision, replace linear interpolation with cubic spline interpolation (requires adding a spline calculation step during initialization).

内容的提问来源于stack exchange，提问作者Museful