I am trying to evaluate the convolution of two continuous functions using scipy.signal.fftconvolve. The scenario of the code is as following:
I am trying to approximate the following double integral:
, i.e. in a region C_1(x',y'), representing a circle of radius 1 centered at (x', y'). This can be approximated by the following integral:
where the function K is chosen as a continuous integrable function, say, exp(-x^2-y^2), the shape of which is approximately that of a circle of radius 1. If I take a function K'(x,y)=K(-x,-y), then the integral is exactly a convolution of the two functions:
So I try to discretize these two functions into arrays and then carry out convolution.
The following code will be written in Julia and the fftconvolve function will be imported using PyCall.jl.
using PyCall
using Interpolations
r = 1
xc = -10:0.05:10
yc = -10:0.05:10
K(x, y) = exp(-(x^2+y^2)/r^2)
rho(x, y) = x^2+y^3 # Try some arbitrary function
ss = pyimport("scipy.signal") # Import scipy.signal module from Python
a = [rho(x,y) for x in xc, y in yc]
b = [K(-x,-y) for x in xc, y in yc]
c = ss.fftconvolve(a,b,mode="same") # zero-paddings beyond boundary, unimportant since rho is near zero beyond the boundary anyway
c_unscaled = interpolate(c', BSpline(Cubic(Line(OnCell()))))
# Adjoint because the array comprehension switched x and y, then interpolate the array
c_scaled = Interpolations.scale(c_unscaled, xc, yc) # Scale the interpolated function w.r.t. xc, yc
print(c_scaled(0.0,0.0)) # The result of the integral for (x', y') = (0, 0)
The result is 628.3185307178969, while the result from numerical integration is 0.785398. What is the problem here?