Comparing convolutions in Mathematica and Python

Question

I'm comparing the results of convolution in Python (using sympy's symbolic variables) and Mathematica with its Convolve function.

In Python, my MWE is

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

which calculates a convolution using symbolic variable x. The resulting function (before lambdifying) is

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

There is no agreement between analytic (fogA, Mathematica) and sympy (fogS, Python):

My Mathematica code is:

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

where

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

and graphs for fogS and fogA are the same as for Python.

Why is there such a large disagreement between the analytic and sympy solution? I suspect the problem lies with sympy. Another Pythonic method is to convolve two arrays which seems to agree with the analytic solution.

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

(Note: this is a MWE. The actual f and g I want to convolve are much more complicated than the Lorentzians defined in this post.)

Mathematica's Convolve function has a bug. See this question and this one explaining the bug. — polfosol ఠ_ఠ

norok2 norok2 · Accepted Answer · 2019-07-01T11:11:17

I do not think this is a reasonable way of using scipy + sympy. I am actually quite surprised that you get a result from lambdify at all.

What you should be doing, instead of using scipy.signal.fftconvolve(), is to use a symbolic definition of the convolution, e.g.:

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

adapted from here.

Comparing convolutions in Mathematica and Python

1 Answers