1

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

enter image description here

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

Medulla Oblongata
  • 3,771
  • 8
  • 36
  • 75
  • _Mathematica_'s `Convolve` function has a bug. See [this question](https://mathematica.stackexchange.com/q/134670/37848) and [this one explaining the bug](https://math.stackexchange.com/q/2082959/301977). – polfosol ఠ_ఠ Jul 03 '19 at 08:42

1 Answers1

1

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.

norok2
  • 25,683
  • 4
  • 73
  • 99
  • Thanks, that's an interesting solution. In terms of speed, your function evaluates the convolution integral itself so is this the most efficient method? I specifically chose to use scipy's `fftconvolve` instead of numpy's `convolve` because I'll eventually be dealing with large arrays (that have up to millions of elements). – Medulla Oblongata Jul 01 '19 at 11:18
  • @MedullaOblongata Once the convolution is computed, it should no longer contain information of the integral (which is the time consuming part). However, this may result in an extremely long expression. Whether such expression, after lambdifying, outperforms a numerical fftconvolve approach or not remains to be proved. However, when plugging in this in your code I get an error from `sympy` which points toward a bug in `sympy`: https://github.com/sympy/sympy/issues/7999 ... bad luck! – norok2 Jul 01 '19 at 11:44
  • Do you also get a `raise PolynomialDivisionFailed(f, g, K)` error? I use sympy v.1.3 and the error message is `PolynomialDivisionFailed: couldn't reduce degree in a polynomial division algorithm when dividing [...] by [...]. This can happen when it's not possible to detect zero in the coefficient domain. The domain of computation is RR[w,_t]. Zero detection is guaranteed in this coefficient domain. This may indicate a bug in SymPy or the domain is user defined and doesn't implement zero detection properly.` – Medulla Oblongata Jul 02 '19 at 07:52
  • You can get rid of that error by defining `x = symbols('x', real=True)`, but you get another one related to the above issue. You can get rid of them all by defining `a` to `d` as `symbols`, at which point I am waiting since yesterday for SymPy to find a solution... perhaps you should try to use Mathematica results instead. – norok2 Jul 02 '19 at 08:04