用于 ODE 的常规时间数据的快速表插值

计算科学 颂 插值 麻木的
2021-12-16 04:11:45

我正在使用scipy.integrate.odeint通过集成来模拟具有已知输入信号的系统的反应。下面的简化代码说明了我在做什么。它通过首先生成数据并在积分函数中对数据进行插值来模拟许多一阶系统对正弦输入的响应。

我想通过利用数据表的时间点是规则间隔的事实来加速模拟，并且我们正在为所有输入变量在表的相同行之间进行插值。分析表明，我正在处理的实际代码确实在numpy.interp.

那么确切的问题是我是否应该推出自己的插值函数来处理这种特殊情况，或者是否存在一个现成的解决方案来有效地处理这种情况。我已经搜索过这个，但空手而归。

import numpy
import scipy.integrate

import matplotlib.pyplot as plt

Ntime = 400
Ninputs = 4
time = numpy.arange(Ntime)

# Generate some out-of-phase sinusoids for input data
DATA = numpy.array([numpy.sin(time/10 - i) for i in range(Ninputs)]).T

# Time constants
taus = numpy.ones(Ninputs)*10

def intfun(x, t):
    # I would like to speed up this lookup:
    inputs = numpy.array([numpy.interp(t, time, DATA[:, i]) 
                          for i in range(Ninputs)])

    return -(1/taus)*x + inputs  # first order

x0 = numpy.ones(Ninputs)

x = scipy.integrate.odeint(intfun, x0, time)

plt.plot(time, x)
plt.show()
1个回答

我在您的方法中看到以下效率低下：

  • 对于 , 的每次调用，每次调用时intfun都numpy.interp必须从零开始搜索DATA用于插值的点。虽然这使用了编译的二进制搜索，但它需要一些时间。

  • 对于 的每次调用intfun，您都有相当大的开销numpy.interp用于插值多个数据点。

  • 您根本不会利用等距的数据。

  • 您必须遍历输入维度 ( Ninputs)，因为numpy.interp仅适用于一维数据。

所以，你可以做什么？

  1. 为避免搜索，您可以利用常规采样并仅使用数组的相关部分进行插值：

    def intfun(x,t):
    j = int(t)
    inputs = numpy.array([
                numpy.interp( t, time[j:j+1], DATA[j:j+1,i] )
                for i in range(Ninputs)
                ])
    return -(1/taus)*x + inputs

    根据经验，这会略微提高速度。

  2. 您可以利用您对常规采样的了解并自己进行插值。这样您就可以摆脱开销，numpy.interp并且可以进行多维插值：

    def intfun(x,t):
    j = int(t)
    t -= j
    inputs = DATA[j]*(1-t) + DATA[j+1]*t
    return -(1/taus)*x + inputs

    根据经验，这提供了比上述更好的速度提升。

  3. 您使用固定步长进行积分（而不是像 那样的自适应积分器odeint）。这允许您以矢量化（因此有效）的方式预先执行插值。但是，您可能需要更多的步骤并以这种方式浪费时间。此外，您还必须找到合适的步长，并且必须为 Python 找到或实现一个固定步长的积分器。

  4. 您切换到具有静态类型的编译语言。这样，您就可以摆脱建议 2 中插值操作的 Python 开销，并且可以同时使用自适应步长。

其它你可能感兴趣的问题
上一篇如何有效地处理 2D 和 3D 模型 下一篇高效计算乙_=一个BX=A当给出的LU 分解时一种A