Сравнение сверток в Mathematica и Python - PullRequest
Новая
       54

Сравнение сверток в Mathematica и Python

0 голосов
/ 01 июля 2019

Я сравниваю результаты свертки в Python ( с использованием символических переменных sympy ) и Mathematica с его функцией Convolve .

В Python мой MWE равен 

from numpy import linspace, pi
from numpy.random import randn
from scipy.signal import fftconvolve
import matplotlib.pyplot as plt
from sympy import symbols
from sympy.utilities.lambdify import lambdify

a = 0.43
b = 0.41
c = 0.65
d = 0.71

x = symbols('x')
f = 2*b / ((x-a)**2 + b**2)
g = 2*d / ((x-c)**2 + d**2)
fog = fftconvolve(f,g,mode='same')
fog_fun = lambdify(x,fog,'numpy') # returns a numpy-ready function
x = linspace(-20,20,int(1e3))
dx = x[1]-x[0]
fogS = fog_fun(x)

fogA = 4*pi*(b+d)/((x-a-c)**2+(b+d)**2) # correct analytic solution

plt.figure()
plt.plot(x,fogA,lw=2,label='analytic')
plt.plot(x,fogS,lw=2,label='sympy')
plt.grid()
plt.legend(loc='best')
plt.show()

, который вычисляет свертку, используя символическую переменную x. Результирующая функция (до лямбда) 

fog = 1.1644/(((x - 0.65)**2 + 0.5041)*((x - 0.43)**2 + 0.1681))

Нет никакого соглашения между analytic (fogA, Mathematica) и sympy (fogS, Python):

enter image description here

Мой код Mathematica: 

a = 0.43; b = 0.41; c = 0.65; d = 0.71;
fogA = FullSimplify[Convolve[2*b/((t-a)^2+b^2),2*d/((t-c)^2+d^2), t, x]];
fogS = 1.1644/(((x - 0.65)^2 + 0.5041)*((x - 0.43)^2 + 0.1681));

, где 

fogA = (17.683+x*(-30.4006+14.0743*x))/(3.04149+x*(-7.9428+x*(8.3428+x*(-4.32+1.*x))))

и графики для fogS и fogA такие же, как для Python.

Почему между решениями analytic и sympy так много разногласий? Я подозреваю, что проблема заключается в sympy. Другой метод Pythonic состоит в том, чтобы объединить два массива , которые, кажется, согласуются с решением analytic. 

f = 2*b / ((x-a)**2 + b**2)
g = 2*d / ((x-c)**2 + d**2)
fogN = fftconvolve(f,g,mode='same')*dx # numeric

(Примечание: это MWE. Фактические f и g, которые я хочу свернуть, гораздо сложнее, чем лоренцы, определенные в этом посте.)

1 Ответ

1 голос
/ 01 июля 2019

Я не думаю, что это разумный способ использования scipy + sympy.Я на самом деле очень удивлен, что вы получаете результат от lambdify.

Что вы должны делать вместо использования scipy.signal.fftconvolve(), это использовать символическое определение свертки, например: 

from sympy import oo, Symbol, integrate

def convolve(f, g, t, lower=-oo, upper=oo):
    tau = Symbol('__very_unlikely_name__', real=True)
    return integrate(
        f.subs(t, tau) * g.subs(t, t - tau), (tau, lower, upper))

адаптировано с здесь .

...