Plotting 3d in python with DiracDelta function

I am trying to plot a two variable function A simple function x*DiracDelta(y)+y*DiracDelta(x) code will work with simple function without DiracDelta But when I enter the delta function, it gives an error? Can you help me find the problem with this code?

from sympy import symbols, DiracDelta
from sympy.plotting import plot_parametric
from sympy import symbols, DiracDelta
from sympy.plotting import plot3d

x, y = symbols('x y')
f = x*DiracDelta(y)+y*DiracDelta(x)

plot3d(f, (x, -1, 1), (y, -1, 1))

Thanks!

try to plot a 3d function It works when I have function without Delta function but it give me error when put delta function

About DiracDelta in general

The Dirac distribution is not a function in the mathematical sense.

The sympy documentation for sympy.DiracDelta states it explicitly:

DiracDelta is not an ordinary function. It can be rigorously defined either as a distribution or as a measure.

DiracDelta only makes sense in definite integrals, [...] care must be taken to not treat it as a real function. [...] if DiracDelta is treated too much like a function, it is easy to get wrong or nonsensical results.

So, you can't plot it like a function, because it isn't a function.

In simple words, the Dirac distribution behaves like a function that is 0 almost everywhere, but has an integral of 1 in any interval that includes 0.

Intuitively, if it could be plotted like a function, the Dirac distribution would have a plot that looks kinda like this:

Except the width of the peak should be 0, and the area under the peak should be exactly 1. This of course would be self-contradictory for a function, so no plot will be truly satisfactory. Again, the Dirac distribution is not a function.

About your f = x*DiracDelta(y)+y*DiracDelta(x)

Again, f is not a function in the mathematical sense. It will behave like a function that is 0 almost everywhere, but with an integral equal to (x2**2 - x1**2)/2 on regions that intersect segment [x1, x2] on the x-axis, and an integral equal to (y2**2 - y1**2)/2 on regions that intersect segment [y1, y2] on the y-axis.

Instead of plotting f, which won't work since f is not a function, I suggest approximating DiracDelta by a "peak" function like the one I drew above, and plot an approximation of f using this approximation of DiracDelta. For instance:

# approximating DiracDelta by a piecewise-affine peak function
from sympy import Piecewise

def dn(x, n):
    return Piecewise(
        (0,               x <= -1/n),
        (n + n**2 * x,    x <= 0),
        (n - n**2 * x,    x <= 1/n),
        (0,               True)
    )


from sympy import symbols, DiracDelta
x, y = symbols('x y')
f = x * DiracDelta(y) + y * DiracDelta(x)

n = 20
fn = x * dn(y, n) + y * dn(x, n)

# now plot fn instead of f