DiracDelta not giving correct result

Question

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.

Agreed, the answer should be one, the area under the delta function being zero and concentrated around x=1 , where the multipler x is not changing anything. — roadrunner66, Apr 20 '16 at 21:30

score 5 · Accepted Answer · edited Apr 21 '16 at 08:25

5

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.

edited Apr 21 '16 at 08:25

user

5,370
8
47
75

answered Apr 20 '16 at 21:57

HuShu

501
6
19

1

**Symbolic** integration (done by sympy) and **numeric** integration (done by scipy) are very different things. Sticking an abstract object like DiracDelta into a numeric integration routine is like sticking fingers into a power outlet; the thing just doesn't belong there. – Apr 20 '16 at 22:02
1

Ah, beat me to it. Basically, scipy's integration is a *numerical* integration, and doesn't play nicely with sympy in general, even more so when we have something as peculiar as the delta "function". – DSM Apr 20 '16 at 22:02
@DSM Your answer was better, please undelete. – Apr 20 '16 at 22:02

score 2 · Answer 2 · answered Apr 21 '16 at 19:57

HuShu's answer is correct. I'll add that the Dirac δ function is a symbolic method of representing function evaluation as an integral. It's useful as a symbolic abstraction, but if you only care about numeric evaluation, just do the function evaluation. That is instead of

b
⌠
⎮ f(x)⋅DiracDelta(x - 1) dx
⌡
a

just use

f(1) if a <= 1 <= b else 0

DiracDelta not giving correct result

2 Answers2

Linked