Я пытаюсь оценить свертку двух непрерывных функций, используя scipy.signal.fftconvolve. Сценарий кода выглядит следующим образом:
Я пытаюсь приблизить следующий двойной интеграл:
, то есть вобласть C_1 (x ', y'), представляющая круг радиуса 1 с центром в (x ', y'). Это может быть аппроксимировано следующим интегралом:
, где функция K выбрана как непрерывная интегрируемая функция, скажем, exp(-x^2-y^2), формаиз которых приблизительно равен окружности радиуса 1. Если я возьму функцию K '(x, y) = K (-x, -y), то интеграл точно является сверткой двух функций:
Поэтому я пытаюсь разделить эти две функции на массивы, а затем выполнить свертку.
Следующий код будет написан на Юлии и *Функция 1028 * будет импортирована с помощью 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)
Результат - 628.3185307178969, а результат численного интегрирования - 0.785398. В чем здесь проблема?