在Python中卷积连续函数的最佳方式

Question

我正在尝试在Python中数值计算形如以下的积分

为此，我首先定义了两个离散的x和t值集合，假设如下：

x_samples = np.linspace(-10, 10, 100)
t_samples = np.linspace(0, 1, 100)
dx = x_samples[1]-x_samples[0]
dt = t_samples[1]-t_samples[0]

声明符号化地，如果t<0，则函数g(x,t)等于0，并将这两个函数离散化以进行积分。

discretG = g(x_samples[None, :], t_samples[:, None])
discretH = h(x_samples[None, :], t_samples[:, None])

我尝试运行了

discretF = signal.fftconvolve(discretG, discretH, mode='full') * dx * dt 

然而，在基本测试函数方面，例如

g(x,t) = lambda x,t: np.exp(-np.abs(x))+t
h(x,t) = lambda x,t: np.exp(-np.abs(x))-t

我发现使用scipy进行数值积分和卷积之间存在不一致，我希望能够找到一种相当快速的计算这些积分的方法，特别是当我只能获得函数的离散表示而非符号表示时。

- Alexandre Krajenbrink

---

你能否编辑你的帖子，附上数值积分结果的示例？ - Yohann L.

你的问题似乎不太清晰。你指的“我在使用scipy进行数值积分和卷积时没有找到一致性”的意思是什么？ - sadegh arefizadeh

Answer 1

根据你的代码，我假设你想在两个函数 g 和 h 上进行卷积，这两个函数只在 [a,b]*[m,n] 上非零。

当然，你可以使用 signal.fftconvolve 来计算卷积。关键是不要忘记在 discretF 内部的指数和实际坐标之间的转换。在这里，我使用插值来计算任意的 (x,t)。

import numpy as np
from scipy import signal, interpolate

a = -1
b = 2
m = -10
n = 15

samples_num = 1000
x_eval_index = 200
t_eval_index = 300

x_samples = np.linspace(a, b, samples_num)
t_samples = np.linspace(m, n, samples_num)
dx = x_samples[1]-x_samples[0]
dt = t_samples[1]-t_samples[0]

g = lambda x,t: np.exp(-np.abs(x))+t
h = lambda x,t: np.exp(-np.abs(x))-t

discretG = g(x_samples[None, :], t_samples[:, None])
discretH = h(x_samples[None, :], t_samples[:, None])

discretF = signal.fftconvolve(discretG, discretH, mode='full')


def compute_f(x, t):
    if x < 2*a or x > 2*b or t < 2*m or t > 2*n:
        return 0
    # use interpolation t get data on new point
    x_samples_for_conv = np.linspace(2*a, 2*b, 2*samples_num-1)
    t_samples_for_conv = np.linspace(2*m, 2*n, 2*samples_num-1)
    f = interpolate.RectBivariateSpline(x_samples_for_conv, t_samples_for_conv, discretF.T)
    return f(x, t)[0, 0] * dx * dt

注意：您可以扩展我的代码，以计算由x和y定义的网格上的卷积，其中x和y是1D数组。（在我的代码中，x和y现在是浮点数）

---

您可以使用以下代码来探索“数值积分”和“使用scipy的卷积”的“一致性”（以及上述compute_f函数的正确性）：

# how the convolve work
# for 1D f[i]=sigma_{j} g[j]h[i-j]
sum = 0
for y_idx, y in enumerate(x_samples[0:]):
    for s_idx, s in enumerate(t_samples[0:]):
        if x_eval_index - y_idx < 0 or t_eval_index - s_idx < 0:
            continue
        if t_eval_index - s_idx >= len(x_samples[0:]) or x_eval_index - y_idx >= len(t_samples[0:]):
            continue
        sum += discretG[t_eval_index - s_idx, x_eval_index - y_idx] * discretH[s_idx, y_idx] * dx * dt
print("Do discrete convolution manually, I get: %f" % sum)
print("Do discrete convolution using scipy, I get: %f" % (discretF[t_eval_index, x_eval_index] * dx * dt))


# numerical integral
# the x_val and t_val
# take 1D convolution as example, function defined on [a, b], and index of your samples range from [0, samples_num-1]
# after convolution, function defined on [2a, 2b], index of your samples range from [0, 2*samples_num-2]
dx_prime = (b-a) / (samples_num-1)
dt_prime = (n-m) / (samples_num-1)
x_eval = 2*a + x_eval_index * dx_prime
t_eval = 2*m + t_eval_index * dt_prime


sum = 0
for y in x_samples[:]:
    for s in t_samples[:]:
        if x_eval - y < a or x_eval - y > b:
            continue
        if t_eval - s < m or t_eval - s > n:
            continue
        if y < a or y >= b:
            continue
        if s < m or s >= n:
            continue
        sum += g(x_eval - y, t_eval - s) * h(y, s) * dx * dt
print("Do numerical integration, I get: %f" % sum)
print("The convolution result of 'compute_f' is: %f" % compute_f(x_eval, t_eval))

这将会给出：

Do discrete convolution manually, I get: -154.771369
Do discrete convolution using scipy, I get: -154.771369
Do numerical integration, I get: -154.771369
The convolution result of 'compute_f' is: -154.771369

- hellohawaii