基于离散坐标向量求二阶导数及ODE迭代误差问题求助

已知深度坐标向量zgrid和对应深度的温度向量Tel，需要计算每个坐标点处温度对位置的二阶导数。由于无法写出温度的解析表达式，无法直接使用ForwardDiff.jl、FiniteDifferences.jl这类工具。

尝试过BSplineKit.jl与Interpolations.jl进行插值求导，但因为结果需要代入长耦合ODE迭代计算，插值方法的误差随迭代不断积累，数千次迭代后结果严重偏离预期。

当前使用的Interpolations.jl代码如下：

function Heatconductivity(κrt,Tel,Tph,zgrid)
    S=Interpolations.interpolate(zgrid,Tel,Gridded(Linear()))
    dS=only.(gradient.(Ref(S),zgrid))
    Q=dS.*kappa_linear(κrt,Tel,Tph)
    Qint=Interpolations.interpolate(zgrid,Q,Gridded(Linear()))
    dQ=only.(gradient.(Ref(Qint),zgrid))
    dQ=dQ.*-1
end

由于每个时间步都会更新温度，无法写出dTel/dz的代数表达式，现寻求可生成“精确”函数以适配ForwardDiff.jl这类工具的方法，解决误差积累问题。

方法1：构造兼容自动微分的离散二阶导数算子 #

不需要依赖插值，直接基于离散网格构造有限差分算子矩阵，并将其封装为可被ForwardDiff.jl识别的函数。这样可以避免插值带来的误差，同时让自动微分工具直接对离散温度向量求导。

针对非均匀网格的二阶中心差分，离散公式为：

d²T/dz² |_i = [ (T_{i+1} - T_i)/(z_{i+1} - z_i) - (T_i - T_{i-1})/(z_i - z_{i-1}) ] / [ (z_{i+1} - z_{i-1})/2 ]

将这个公式封装为函数，所有运算默认兼容ForwardDiff：

function second_derivative(zgrid, Tel)
    n = length(zgrid)
    d2T = similar(Tel)
    # 内部点用中心差分
    for i in 2:n-1
        dz1 = zgrid[i] - zgrid[i-1]
        dz2 = zgrid[i+1] - zgrid[i]
        dz_avg = (dz1 + dz2)/2
        dT_dz_plus = (Tel[i+1] - Tel[i])/dz2
        dT_dz_minus = (Tel[i] - Tel[i-1])/dz1
        d2T[i] = (dT_dz_plus - dT_dz_minus)/dz_avg
    end
    # 左边界i=1，用二阶单边差分
    dz1 = zgrid[2] - zgrid[1]
    dz2 = zgrid[3] - zgrid[1]
    d2T[1] = (Tel[3] - 2Tel[2] + Tel[1])/(dz1*dz2) * dz2
    # 右边界i=n，用二阶单边差分
    dz1 = zgrid[n] - zgrid[n-1]
    dz2 = zgrid[n] - zgrid[n-2]
    d2T[n] = (Tel[n] - 2Tel[n-1] + Tel[n-2])/(dz1*dz2) * dz2
    return d2T
end

同时配套实现一阶离散导数函数，修改Heatconductivity：

function first_derivative(zgrid, Tel)
    n = length(zgrid)
    dT = similar(Tel)
    for i in 2:n-1
        dz1 = zgrid[i] - zgrid[i-1]
        dz2 = zgrid[i+1] - zgrid[i]
        dT[i] = (Tel[i+1]*dz1 + Tel[i-1]*dz2) / (dz1*dz2) - Tel[i]*(dz1+dz2)/(dz1*dz2)
    end
    # 边界点用单边差分
    dT[1] = (Tel[2] - Tel[1])/(zgrid[2]-zgrid[1])
    dT[end] = (Tel[end] - Tel[end-1])/(zgrid[end]-zgrid[end-1])
    return dT
end

function Heatconductivity(κrt,Tel,Tph,zgrid)
    dS = first_derivative(zgrid, Tel)
    Q = dS .* kappa_linear(κrt, Tel, Tph)
    dQ = first_derivative(zgrid, Q)
    return -dQ
end

方法2：使用可微分的高精度插值库 #

如果仍想保留插值方式，选择支持自动微分的高阶插值实现，比如BSplineKit.jl的可微分接口。高阶B样条的插值误差远低于线性插值，且运算兼容ForwardDiff，可有效减少误差积累：

using BSplineKit, ForwardDiff

function Heatconductivity(κrt,Tel,Tph,zgrid)
    # 构造3阶B样条插值（阶数可根据需求调整）
    spl_T = interpolate(zgrid, Tel, BSplineOrder(3))
    dS = derivative.(spl_T, zgrid)
    
    Q = dS .* kappa_linear(κrt, Tel, Tph)
    spl_Q = interpolate(zgrid, Q, BSplineOrder(3))
    dQ = derivative.(spl_Q, zgrid)
    
    return -dQ
end

关键注意点 #

• 离散差分算子是无插值近似的“精确”离散形式，更适合ODE迭代场景；
• 若使用插值，优先选择3阶及以上的B样条插值，避免线性插值的低精度缺陷；
• 所有基础数组运算均兼容ForwardDiff，无需手动编写解析导数逻辑。

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