C++实现Monotone Cubic Interpolation时出现vector下标越界求助

算力阿峰

2026-7-1

单调三次插值代码下标越界问题排查

我写了一段实现**单调三次插值（Monotone Cubic Interpolation）**的C++代码，用给定的x、y向量生成间隔0.0004的查询点xq时，出现vector subscript out of range错误。但用x从0到5、y从0到4的测试数据时代码能正常运行，求帮忙解决下标越界问题。

原代码

#include <iostream>
#include <vector>
#include <algorithm>

// Function declarations
std::vector<double> wiki_mon_spline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& xq);

int main() {
    // Example usage
    std::vector<double> x = { 0,0.0208,0.0416,0.2,0.2208,0.2416,0.2416,0.2416,1.2416,1.2416,1.2624,1.2832,1.6416,1.6624,1.6832,1.6832,1.6832,2.6832,2.6832,2.704,2.7248,3.2832,3.304,3.3248,3.3248,3.3248,3.8248 };
    std::vector<double> y = { 0.404674,34.3471,31.8996,31.8996,6.74855,0.404674,0.404674,0.404674,0.404674,0.404674,30.3726,31.8996,31.8996,7.07452,0.404674,0.404674,0.404674,0.404674,0.404674,29.053,30.7901,30.7901,7.67308,0.404674,0.404674,0.404674,0.404674 };

    // Calculate xq in C++
    double dt = 0.0004;
    std::vector<double> xq;

    // Generate 10,000 equidistant samples
    double min_x = *std::min_element(x.begin(), x.end());
    double max_x = *std::max_element(x.begin(), x.end());
    for (double val = min_x; val <= max_x; val += dt)
    {
        xq.push_back(val);
    }

    std::cout << "Size of xq: " << xq.size() << std::endl;
    // Perform monotone cubic interpolation
    std::vector<double> result = wiki_mon_spline(x, y, xq);

    // Print the result
    for (size_t i = 0; i < xq.size(); ++i) {
        std::cout << "Interpolation at xq[" << i << "]: " << xq[i] << "     " << result[i] << std::endl;
    }

    return 0;
}

std::vector<double> wiki_mon_spline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& xq) {
    int n = x.size();

    // Calculate slope at grid points
    std::vector<double> h(n - 1);
    std::vector<double> delta(n - 1);

    for (int i = 0; i < n - 1; ++i) {
        h[i] = x[i + 1] - x[i];
        delta[i] = (y[i + 1] - y[i]) / h[i];
    }

    // Calculate slopes at interior points
    std::vector<double> m(n, 0.0);
    for (int i = 0; i < n - 1; ++i) {
        if (delta[i] * delta[i + 1] > 0) {
            double w1 = 2 * h[i + 1] + h[i];
            double w2 = h[i + 1] + 2 * h[i];
            m[i + 1] = (w1 + w2) / (w1 / delta[i] + w2 / delta[i + 1]);
        }
        else {
            m[i + 1] = 0.0; // Initialize to 0 if delta[i] and delta[i + 1] have opposite signs
        }
    }

    // Slope at end points
    m[0] = ((2 * h[0] + h[1]) * delta[0] - h[0] * delta[1]) / (h[0] + h[1]);

    if (m[0] * delta[0] <= 0) {
        m[0] = 0.0;
    }
    else if ((delta[0] * delta[1] <= 0) && (std::abs(m[0]) > std::abs(3 * delta[0]))) {
        m[0] = 3 * delta[0];
    }

    // Perform interpolation at query points xq
    std::vector<double> yq(xq.size(), 0.0);

    for (size_t i = 0; i < xq.size(); ++i) {
        // Find the interval containing xq[i]
        auto it = std::lower_bound(x.begin(), x.end(), xq[i]);
        if (it != x.end() && it != x.begin()) {
            int j = static_cast<int>(it - x.begin()) - 1;
            if (j < n - 1 && j >= 0) {
                // PCHIP interpolation formula
                double dx = xq[i] - x[j];
                double t = dx / h[j];
                double h00 = 2 * t * t * t - 3 * t * t + 1;
                double h10 = t * t * t - 2 * t * t + t;
                double h01 = -2 * t * t * t + 3 * t * t;
                double h11 = t * t * t - t * t;

                // Ensure monotonicity by considering signs
                double sign_delta = (delta[j] >= 0.0) ? 1.0 : -1.0;
                double sign_delta_next = (delta[j + 1] >= 0.0) ? 1.0 : -1.0;

                yq[i] = h00 * y[j] + sign_delta * h10 * h[j] + h01 * y[j + 1] + sign_delta_next * h11 * m[j];
            }
            else if (j < 0) {
                // Handle the case where xq[i] is beyond the first element of x
                yq[i] = y[0];
            }
            else {
                // Handle the case where xq[i] is beyond the last element of x
                yq[i] = y[n - 1];
            }
        }
    }

    return yq;
}

错误原因分析

重复x值导致除数为0：x向量存在连续重复值（如0.2416连续出现3次），计算delta[i]时h[i] = x[i+1]-x[i]为0，导致delta[i]变成NaN，后续计算斜率时引发异常，进而触发下标越界。
斜率计算数组越界：计算内部斜率的循环for (int i = 0; i < n - 1; ++i)中，访问了delta[i+1]，当i = n-2时，i+1 = n-1，但delta的大小是n-1，最大下标为n-2，直接越界访问。
查询点边界处理不全：std::lower_bound返回x.end()时（即xq[i]大于所有x值），代码未处理该情况，导致部分查询点未被赋值，后续访问时可能引发未定义行为。
插值公式错误：原代码中错误使用m[j]替代m[j+1]，且多余的符号处理不符合PCHIP插值逻辑。

修复后的代码

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>

// Function declarations
std::vector<double> wiki_mon_spline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& xq);

int main() {
    // Example usage
    std::vector<double> x = { 0,0.0208,0.0416,0.2,0.2208,0.2416,0.2416,0.2416,1.2416,1.2416,1.2624,1.2832,1.6416,1.6624,1.6832,1.6832,1.6832,2.6832,2.6832,2.704,2.7248,3.2832,3.304,3.3248,3.3248,3.3248,3.8248 };
    std::vector<double> y = { 0.404674,34.3471,31.8996,31.8996,6.74855,0.404674,0.404674,0.404674,0.404674,0.404674,30.3726,31.8996,31.8996,7.07452,0.404674,0.404674,0.404674,0.404674,0.404674,29.053,30.7901,30.7901,7.67308,0.404674,0.404674,0.404674,0.404674 };

    // Calculate xq in C++
    double dt = 0.0004;
    std::vector<double> xq;

    // Generate equidistant samples
    double min_x = *std::min_element(x.begin(), x.end());
    double max_x = *std::max_element(x.begin(), x.end());
    for (double val = min_x; val <= max_x; val += dt)
    {
        xq.push_back(val);
    }

    std::cout << "Size of xq: " << xq.size() << std::endl;
    // Perform monotone cubic interpolation
    std::vector<double> result = wiki_mon_spline(x, y, xq);

    // Print the result (optional, comment out for large xq)
    // for (size_t i = 0; i < xq.size(); ++i) {
    //     std::cout << "Interpolation at xq[" << i << "]: " << xq[i] << "     " << result[i] << std::endl;
    // }

    return 0;
}

std::vector<double> wiki_mon_spline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& xq) {
    int n = x.size();
    if (n < 2) {
        return std::vector<double>(xq.size(), n == 1 ? y[0] : 0.0);
    }

    // Calculate slope at grid points
    std::vector<double> h(n - 1);
    std::vector<double> delta(n - 1);

    for (int i = 0; i < n - 1; ++i) {
        h[i] = x[i + 1] - x[i];
        // Handle duplicate x values to avoid division by zero
        if (std::abs(h[i]) < 1e-10) {
            delta[i] = 0.0;
        } else {
            delta[i] = (y[i + 1] - y[i]) / h[i];
        }
    }

    // Calculate slopes at interior points
    std::vector<double> m(n, 0.0);
    // 修正循环边界：避免访问delta[i+1]越界
    for (int i = 0; i < n - 2; ++i) {
        if (std::abs(delta[i]) > 1e-10 && std::abs(delta[i + 1]) > 1e-10 && delta[i] * delta[i + 1] > 0) {
            double w1 = 2 * h[i + 1] + h[i];
            double w2 = h[i + 1] + 2 * h[i];
            m[i + 1] = (w1 + w2) / (w1 / delta[i] + w2 / delta[i + 1]);
        } else {
            m[i + 1] = 0.0;
        }
    }

    // Slope at first point
    if (std::abs(delta[0]) > 1e-10) {
        if (n == 2) {
            m[0] = delta[0];
        } else {
            m[0] = ((2 * h[0] + h[1]) * delta[0] - h[0] * delta[1]) / (h[0] + h[1]);
            if (m[0] * delta[0] <= 0) {
                m[0] = 0.0;
            } else if ((delta[0] * delta[1] <= 0) && (std::abs(m[0]) > std::abs(3 * delta[0]))) {
                m[0] = 3 * delta[0];
            }
        }
    }

    // Slope at last point
    if (std::abs(delta[n-2]) > 1e-10) {
        if (n == 2) {
            m[n-1] = delta[n-2];
        } else {
            m[n-1] = ((2 * h[n-2] + h[n-3]) * delta[n-2] - h[n-2] * delta[n-3]) / (h[n-2] + h[n-3]);
            if (m[n-1] * delta[n-2] <= 0) {
                m[n-1] = 0.0;
            } else if ((delta[n-2] * delta[n-3] <= 0) && (std::abs(m[n-1]) > std::abs(3 * delta[n-2]))) {
                m[n-1] = 3 * delta[n-2];
            }
        }
    }

    // Perform interpolation at query points xq
    std::vector<double> yq(xq.size(), 0.0);

    for (size_t i = 0; i < xq.size(); ++i) {
        double xi = xq[i];
        // Find the interval containing xi
        auto it = std::lower_bound(x.begin(), x.end(), xi);
        
        if (it == x.begin()) {
            // xi <= x[0], use first point's value
            yq[i] = y[0];
        } else if (it == x.end()) {
            // xi >= x.back(), use last point's value
            yq[i] = y.back();
        } else {
            int j = static_cast<int>(it - x.begin()) - 1;
            // Ensure j is within valid range
            if (j < 0) j = 0;
            if (j >= n-1) j = n-2;

            double dx = xi - x[j];
            // Handle duplicate x in interval
            if (std::abs(h[j]) < 1e-10) {
                yq[i] = y[j];
                continue;
            }
            double t = dx / h[j];
            double h00 = 2 * t * t * t - 3 * t * t + 1;
            double h10 = t * t * t - 2 * t * t + t;
            double h01 = -2 * t * t * t + 3 * t * t;
            double h11 = t * t * t - t * t;

            // Correct PCHIP interpolation formula
            yq[i] = h00 * y[j] + h10 * m[j] * h[j] + h01 * y[j+1] + h11 * m[j+1] * h[j];
        }
    }

    return yq;
}

C++实现Monotone Cubic Interpolation时出现vector下标越界求助

原代码 #

错误原因分析 #

修复后的代码 #

原代码

错误原因分析

修复后的代码