如何用scipy.interpolate.interp1d实现对数均匀间隔数据插值？

AI 小北

2026-7-27

问题描述

我有一个实值函数$G(t)$，需在对数均匀间隔的时间点$t \in [10^{-2}, 10^{-1.9}, 10^{-1.8}, ..., 10^{6}]$上插值，以得到$G$的分量关于频率$\omega = \frac{1}{t}$的近似表达式。通过以下代码定义时间列表（因numpy.logspace报错，改用numpy.linspace）：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

exponents = np.linspace(-2, 6, 81, endpoint=True) #81 values
ts = [10**(xp) for xp in exponents] #times (s)

预期结果类似下图：
$Interpolated components of $G$ depending on $\omega$$

尝试用scipy.interpolate.interp1d编写G_pr_sec函数，但调用G_pr_sec(ts, Gs)时出现错误：ValueError: A value in x_new is below the interpolation range.
函数代码如下：

def G_pr_sec(ts, Goft):
    """
    Returns G_prime and G_sec at frequencies 1/ts[i] by interpolating the array Goft

    Parameters
    ----------
    ts : 1D array
        Array of times
    Goft : 1D array
        Values of G(t)
        
    Returns
    -------
    G_pr : 1D array
        Array of G_prime(1/t)
    G_sec : 1D array
        Array of G_second(1/t)

    """
    G2 = interp1d(ts, Goft)
    
    G_pr = []
    G_sec = []
    
    for t in ts:
        res_sec = 0
        res_sec += -0.47*(G2(2*t) - G2(4*t)) + 1.674*(G2(t) - G2(2*t)) + 0.198*(G2(t/2) - G2(t))
        res_sec += 0.62*(G2(t/4) - G2(t/2)) + 0.012*(G2(t/8) - G2(t/4)) + 0.172*(G2(t/16) - G2(t/8))
        res_sec += 0.0433*(G2(t/64) - G2(t/32)) + 0.0108*(G2(t/256) - G2(t/128))
        res_sec += (0.0108/4)*(G2(t/1024) - G2(t/512)) + (0.0108/16)*(G2(t/4096) - G2(t/2048))
        res_sec += (0.0108/64)*(G2(t/(64*256)) - G2(t/(32*256))) + (0.0108/256)*(G2(t/(256**2)) - G2(t/(128*256)))
        G_sec.append(res_sec)
        
        res_pr = 0
        res_pr += G2(t) -0.142*(G2(4*t) - G2(8*t)) + 0.717*(G2(2*t) - G2(4*t)) + 0.046*(G2(t) - G2(2*t))
        res_pr += 0.099*(G2(t/2) - G2(t)) + 0.103*(G2(t/4) - G2(t/2)) + 0.001*(G2(t/8) - G2(t/4))
        res_pr += 0.00716*(G2(t/16) - G2(t/8)) + 0.000451*(G2(t/64) - G2(t/32))
        G_pr.append(res_pr)
        
    return [G_pr, G_sec]

尝试用自定义log_interp1d函数替换interp1d后仍出现相同错误：

import scipy as sp
import scipy.interpolate

def log_interp1d(xx, yy, kind='linear'):
    logx = np.log10(xx)
    logy = np.log10(yy)
    lin_interp = sp.interpolate.interp1d(logx, logy, kind=kind)
    log_interp = lambda zz: np.power(10.0, lin_interp(np.log10(zz)))
    return log_interp

希望在保持对数均匀间隔时间（及对应频率）的前提下完成插值，避免线性均匀间隔导致的对数尺度下数据分布不均问题。

解决方案

错误根源

代码中计算G_sec和G_pr时，用到了t/4096、t/(256²)这类极小时间值，远小于你定义的ts最小值$10^{-2}$，超出了插值函数的输入范围，因此触发越界错误。

解决方法

1. 扩展时间范围（优先推荐）

计算所有需要查询的t值的最小值，将ts的范围扩展到覆盖这些极小值。例如当t取最小的$10^{-2}$时，t/(256²)≈1.525e-7，因此可以将指数范围从-2扩展到-7，生成包含更小时间点的对数均匀序列：

# 扩展指数范围，覆盖所有需要查询的t值，共131个对数均匀点
exponents = np.linspace(-7, 6, 131, endpoint=True)
ts = 10**exponents  # 直接用numpy数组替代列表推导式，提升效率

扩展后，所有需要插值的t值都落在合法范围内，不会再触发越界错误。

2. 允许外插（备选方案）

若无法扩展时间范围，可以给interp1d添加fill_value="extrapolate"参数，允许函数对超出范围的值进行外插计算：

G2 = interp1d(ts, Goft, fill_value="extrapolate", kind='linear')

但注意：外插小时间（高频）区域的数据可能引入较大误差，对于流变学中松弛模量转储能/损耗模量的场景，外插结果的可靠性需要谨慎评估。

3. 优化代码效率

将原函数中的列表循环改为numpy数组操作，大幅提升运算速度：

def G_pr_sec(ts, Goft):
    """
    Returns G_prime and G_sec at frequencies 1/ts[i] by interpolating the array Goft

    Parameters
    ----------
    ts : 1D numpy array
        Array of times
    Goft : 1D numpy array
        Values of G(t)
        
    Returns
    -------
    G_pr : 1D numpy array
        Array of G_prime(1/t)
    G_sec : 1D numpy array
        Array of G_second(1/t)
    """
    # 若已扩展ts范围，无需外插；否则添加fill_value="extrapolate"
    G2 = interp1d(ts, Goft, kind='linear')
    
    ts_arr = np.asarray(ts)
    
    # 计算G_sec
    res_sec = np.zeros_like(ts_arr)
    res_sec += -0.47*(G2(2*ts_arr) - G2(4*ts_arr)) + 1.674*(G2(ts_arr) - G2(2*ts_arr)) + 0.198*(G2(ts_arr/2) - G2(ts_arr))
    res_sec += 0.62*(G2(ts_arr/4) - G2(ts_arr/2)) + 0.012*(G2(ts_arr/8) - G2(ts_arr/4)) + 0.172*(G2(ts_arr/16) - G2(ts_arr/8))
    res_sec += 0.0433*(G2(ts_arr/64) - G2(ts_arr/32)) + 0.0108*(G2(ts_arr/256) - G2(ts_arr/128))
    res_sec += (0.0108/4)*(G2(ts_arr/1024) - G2(ts_arr/512)) + (0.0108/16)*(G2(ts_arr/4096) - G2(ts_arr/2048))
    res_sec += (0.0108/64)*(G2(ts_arr/(64*256)) - G2(ts_arr/(32*256))) + (0.0108/256)*(G2(ts_arr/(256**2)) - G2(ts_arr/(128*256)))
    
    # 计算G_pr
    res_pr = np.zeros_like(ts_arr)
    res_pr += G2(ts_arr) -0.142*(G2(4*ts_arr) - G2(8*ts_arr)) + 0.717*(G2(2*ts_arr) - G2(4*ts_arr)) + 0.046*(G2(ts_arr) - G2(2*ts_arr))
    res_pr += 0.099*(G2(ts_arr/2) - G2(ts_arr)) + 0.103*(G2(ts_arr/4) - G2(ts_arr/2)) + 0.001*(G2(ts_arr/8) - G2(ts_arr/4))
    res_pr += 0.00716*(G2(ts_arr/16) - G2(ts_arr/8)) + 0.000451*(G2(ts_arr/64) - G2(ts_arr/32))
    
    return res_pr, res_sec

内容的提问来源于stack exchange，提问作者Victor G.