Optimal way to convolute continuous functions in python

user_19577021 0
Viewed 214

I am trying to numerically compute in python integrals of the form

2D convolution

To that aim, I first define two discrete sets of x and t values, let's say

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]

declare symbolically that the function g(x,t) is equal to 0 if 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])

I have then tried to run

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

Yet, on basic test functions such as

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

I don't find an agreement between the the numerical integration and the convolution using scipy and I would like to have a fairly fast way of computing these integrals, especially when I only have access to discretised representations of the functions rather than their symbolic one.

convolution scipy integration python
1 Answers

According to your code, I assume you want to conduct convolution on two function g and h that are non-zero only on [a, b]*[m,n].

Of course you can use signal.fftconvolve to compute the convolution. The key is don't forget the transformation between the indices inside discretF and the real coordinates. Here I use interpolation to compute for arbitrary (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

Note: you can extend my codes to compute convolution on a meshgrid defined by x and y, where x and y are 1D array. (In my code, x and y are float now)

You can use the following code to explore the "agreement" between "the numerical integration" and "the convolution using scipy" (and also, the correctness of compute_f function above)：

# 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))

Which gives:

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
user_97587900
Related