小能豆

Python numpy.convolve 积分极限

py

我有一个卷积积分,其形式与本文描述的类似:(链接

我需要从负无穷到 t 进行积分,而不是从 0 到 t 进行积分。但是,使用 我无法做到这一点,numpy.convolve因为它总是返回从负无穷到正无穷的结果。使用scipy.integrate.quad会非常慢,因为我必须循环遍历每个t,并且它仅适用于具有解析表达式的被积函数。

有没有办法指定 的下限和上限numpy.convolve?非常感谢。

以下是代码(抱歉无法在这里输入 LaTeX 方程式):

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

def gaussian(tau, mu, sigma):

    pdf = 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (tau - mu)**2 / (2 * sigma**2))
    return pdf

def gaussian_deriv(tau, mu, sigma):
    """
    derivative of gaussian function
    """
    pdf = -(tau - mu)/sigma**2 * gaussian(tau, mu, sigma)
    return pdf

def integral_kernel(tau):
    return np.cbrt(1/tau)

def integrand(tau, t, mu, sigma):

    return gaussian_deriv(tau, mu, sigma) * integral_kernel(t - tau + 1E-28)

tau = np.linspace(-7, 7, 1000)
dtau = tau[1] - tau[0]
lower_lim = tau[0]
g_deriv = gaussian_deriv(tau, mu=0, sigma=1)

result = np.zeros_like(tau)
for idx, t in np.ndenumerate(tau):
    result[idx], err = integrate.quad(integrand, lower_lim, t, args=(t, mu, sigma), points=[t])

result_convolve = np.convolve(g_deriv, integral_kernel(tau), mode='same') * dtau

fig, ax = plt.subplots(2,1, figsize=(10, 6))
ax[0].plot(tau, result, 'r-', label='scipy quad')
ax[1].plot(tau, result_convolve, '.', label='numpy convolve')
ax[0].legend()
ax[1].legend()
plt.show()

阅读 25

收藏
2024-12-24

共1个答案

小能豆

如果我理解正确的话,你想要积分

int{-inf…T} f(tau) g(t-tau) dtau

写为 H(t) = 1/2 (sign(t) + 1),这等同于

int{-inf…inf} f(tau) g(t-tau) [1 - H(tau-T)] dtau

或者

int{-inf…inf} f^(tau) g(t-tau) dtau = f^*g(t)

其中 f^(tau) = f(tau) * [1 - H(tau-T)]。

因此您需要做的就是将 T 右侧的 f 清零以获得 f^,然后计算常规卷积 f^*g。

更新

如果 t 和 T 是相同的数字,那就更简单了:

int{-inf…t} f(tau)g(t-tau) dtau

在这种情况下,你可以将积分变量移动 t,这样你就会得到:

int{-inf…0} f(tau+t}g(-tau) dtau

接下来用 -tau 替换 tau

int{0…inf} g(tau) f(t-tau) dtau

现在您或多或少地按照上面的方法操作,只是 T 现在是固定的并且等于零,并且您不是将 T 的右侧清零,而是将 T 的左侧清零。

2024-12-24