均值保持二次插值的不稳定实现问题问询

我一直在寻找高效的均值保持二次插值计算方法，但始终没有找到合适的方案。于是我用Python构建了一个线性方程组系统，可即便输入完全平坦的数组，该系统仍存在数值不稳定的问题。

模型设定 #

每个数据分箱宽度为1，分箱值等于二次多项式在该分箱上的积分，即：

Int(a_i+b_ix+c_ix^2,{dx,x-1/2,x+1/2})=data[x,i, x=i]=a_i+b_ix+c_i*(x^2+1/12)

同时设置了相邻分箱间的两个匹配条件：

• 数值匹配：

a_i+b_i*(x+1/2)+c_i*(x+1/2)^2=a_{i+1}+b_{i+1}*(x+1/2)+c_{i+1}*(x+1/2)^2

• 斜率匹配条件（对应多项式导数在分箱边界相等）

样条节点绑定在±1/2处，第一个分箱的负端自由，最后一个分箱定义边界条件。对于长度为n的数据集，我构建了3n个方程和3n个变量，np.linalg.solve()可以处理该系统，但所有输入数据集的结果均不稳定，输入[1,2,3]时直接失效，方程组根本无法正确求解。

我期望integro_quadratic_spline_orig函数在输入integro_quadratic_spline(I)的输出后能返回原数组I，但实际结果偏差极大。

实现代码 #

import numpy as np
import scipy as sp

def integro_quadratic_spline(I):
    # 区间数量
    n = len(I)

    # 初始化系数数组
    a = np.zeros(n)
    b = np.zeros(n)
    c = np.zeros(n)

    # 构建方程组矩阵与结果向量
    A = sp.sparse.csr_array((3*n, 3*n))
    B = np.zeros(3*n)

    # 填充方程组：积分条件、斜率连续、数值连续
    for i in range(0,n-1):
        # 积分匹配条件：分箱积分等于原数据值
        A[3*(i),3*(i):3*(i)+6]=[1,(i),(i)**2+1/12,0,0,0]
        # 斜率连续条件：当前分箱右边界斜率等于下一分箱左边界斜率
        A[3*(i)+1,3*(i):3*(i)+6]=[0,1,2*(i+1/2),0,-1,-2*(i+1/2)]
        # 数值连续条件：当前分箱右边界值等于下一分箱左边界值
        A[3*(i)+2,3*(i):3*(i)+6]=[1,i+1/2,(i+1/2)**2,-1,-(i+1/2),-(i+1/2)**2]
        B[3*(i)]=I[i]
        B[3*(i)+1]=0
        B[3*(i)+2]=0
    # 处理最后一个分箱的边界条件
    i=i+1    
    A[3*(i),3*(i):3*(i)+3]=[1,(i),(i)**2+1/12]
    A[3*(i)+1,3*(i):3*(i)+3]=[0,1,2*(i+1/2)]
    A[3*(i)+2,3*(i):3*(i)+3]=[1,i+1/2,(i+1/2)**2]
    B[3*(i)]=I[i]
    B[3*(i)+1]=0
    B[3*(i)+2]=I[i]

    # 求解方程组
    coeffs = sp.sparse.linalg.spsolve(A, B)
    print(np.allclose(np.dot(A.toarray(), coeffs), B))

    # 提取各分箱的多项式系数
    for i in range(n):
        a[i] = coeffs[3*i]
        b[i] = coeffs[3*i+1]
        c[i] = coeffs[3*i+2]

    return a, b, c


def integro_quadratic_spline_return(coeffs, mag):
    original_length = len(coeffs[0])
    new_length=original_length*mag
    newplot=np.empty(new_length)
    # 根据系数生成高分辨率插值结果
    for i in range(original_length):
        for j in range(int(mag)):
            newplot[mag*i+j]=coeffs[0][i]+coeffs[1][i]*(i+j/mag-1/2)+coeffs[2][i]*(i+j/mag-1/2)**2

    return newplot

def integro_quadratic_spline_orig(coeffs):
    original_length = len(coeffs[0])
    newplot=np.empty(original_length)
    # 尝试从系数还原原分箱值
    for i in range(original_length):
        newplot[i]=coeffs[0][i]+coeffs[1][i]*(i)+coeffs[2][i]*((i)**2-12)

    return newplot

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