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]

象征性地声明函数 g(x,t) 等于 0 如果 t<0 and discretise the two functions to integrate as

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 进行数值积分和卷积之间的一致性，我希望有一种相当快的方法来计算这些积分，特别是当我只能访问函数的离散表示而不是它们的符号表示时.

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

是一维数组。 （在我的代码中，

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