I have to use dirac delta in a complicated integral and was hoping to see how it works with a simple case but it returns the wrong answer. Any clue what I did wrong in the following?
from sympy import DiracDelta
from scipy import integrate
def f(x):
return x*DiracDelta(x-1)
b, err = integrate.quad(f, 0, 5)
print b
This returns 0.0 while it shouldn't.
It seems sympy functions are not compatible with scipy integrate. One needs to use sympy integrate. The following gives the correct answer
from sympy import *
x = Symbol('x')
print integrate(x*DiracDelta(x-1), (x, 0, 5.0))
Although, I am not sure if sympy.integrate is as powerful and versatile as scipy.integrate.