plot DiracDelta(x)

[hor...@gmail.com]

Hello,
I just started using sympy for my electrical engineering studies, and already have my first question:

How is possible du plot a DiracDelta(x) ?

Cheers,
David

[hor...@gmail.com]

Here is my python notebook script

while the last plot the ipython kernel seems to be busy all the time

==============
import sympy
from sympy.plotting import plot
from sympy.abc import t,tau

u_t = 2 * DiracDelta(t)
g_t = 0.5 * ( Heaviside(t)-Heaviside(t-2) )

plot(g_t, (t, -5, 5))
==============

plot( integrate( u_t.subs(t,tau)*g_t.subs(t,t-tau), (tau,-10, 10)), (t,-10,10))
plot( integrate( u_t.subs(t,tau)*u_t.subs(t,t-tau), (tau,-10, 10)), (t,-10,10))
plot( integrate( g_t.subs(t,tau)*g_t.subs(t,t-tau), (tau,-10, 10)), (t,-10,10))

[Stefan Krastanov]

Below I try to answer the question (basically there is no way for the moment), but could you tell us how you actually expected this figure to look like?

There does not seem to be a satisfactory way to do it. After all DiracDelta is zero at all points besides x=0, so simple sampling will miss this point. If x=0 is actually sampled, it returns infinity which either will cause the code to crash or it will mask the point as a NaN so it will be invisible.

A step forward to solving this would be to have some kind of asymptote detector that changes the scale of the plot accordingly.

Another way to solve this would be to provide an appropriate api for functions to tell the plotting module that they have such degenerate points.

In any case, a lot of work is necessary for the implementation of a satisfactory solution.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Aaron Meurer]

On Sun, Jul 28, 2013 at 11:23 AM, Stefan Krastanov
<krastano...@gmail.com> wrote:
> Below I try to answer the question (basically there is no way for the
> moment), but could you tell us how you actually expected this figure to look
> like?
>
> There does not seem to be a satisfactory way to do it. After all DiracDelta
> is zero at all points besides x=0, so simple sampling will miss this point.
> If x=0 is actually sampled, it returns infinity which either will cause the
> code to crash or it will mask the point as a NaN so it will be invisible.
>
> A step forward to solving this would be to have some kind of asymptote
> detector that changes the scale of the plot accordingly.
>
> Another way to solve this would be to provide an appropriate api for
> functions to tell the plotting module that they have such degenerate points.

This would be the ideal solution, so that we could print some kind of
vertical line, like at http://en.wikipedia.org/wiki/Dirac_delta. It
would also allow plotting things like 1/x with dashed lines at that
asymptote, and even non-vertical asymptotes. I like the idea of the
API, so that we could have this be both automatic based on solve() and
extensible. This would be useful outside of plotting as well.

For now, probably the easiest way to do this would be to find some
suitable fn(x) such that fn(x) -> DiracDelta(x) as n -> oo and plot
fn(x) for n pretty large. Stefan will have to comment on what kind of
fn will be best suited for the plotting module.

Aaron Meurer

[hor...@gmail.com]

Ok, i see this isn't implemented yet. But could you give me a hint why while the last plot the kernel is busy all the time ?

import sympy
from sympy.plotting import plot
from sympy.abc import t,tau

u_t = 2 * DiracDelta(t)
g_t = 0.5 * ( Heaviside(t)-Heaviside(t-2) )

[Stefan Krastanov]

It is because the integral can not be solved symbolically by sympy so it has to be evaluated numerically _for each point t that is plotted_.

first plot: evaluate the integral symbolically (slow), calculate the value of the result at each t (fast as it is only arithmetics)

last plot: try to evaluate the integral symbolically (slow and useless as it fails), fallback to numeric integrations at each t (slow because numeric integration needs a lot of sampling)

--

[Aaron Meurer]

On Mon, Jul 29, 2013 at 3:20 AM, Stefan Krastanov
<krastano...@gmail.com> wrote:
> It is because the integral can not be solved symbolically by sympy so it has
> to be evaluated numerically _for each point t that is plotted_.
>
> first plot: evaluate the integral symbolically (slow), calculate the value
> of the result at each t (fast as it is only arithmetics)
>
> last plot: try to evaluate the integral symbolically (slow and useless as it
> fails), fallback to numeric integrations at each t (slow because numeric
> integration needs a lot of sampling)

Can we modify the numerical integration routine so that it can give a
list of points more efficiently? It seems like there would be a lot of
duplicate sampling happening. Or maybe we just need to memoize the
sampling of the integrand.

Aaron Meurer

[Stefan Krastanov]

The fallback is just to call `evalf` instead something like `lambdify`. It is always slower, but works even on the most bizarre expressions. For integrals, indeed, there are many points that are resampled with this naive solution (the algorithm becomes n^2 instead of n).

Numeric libraries like scipy provide routines for doing this in a single pass, however one provides the points to be sampled beforehand. mpmath which is used in this case does not provide this as far as I know.

Even if it is provided we will have to somehow link this routine to `Integral.evalf`, because if we just write a special case for the plotting module, it will work for `Integral(...)` but not for something more general like `exp(Integral(...))` or `x*Integral(...)`.

So yes, we can implement the single-pass algorithm, but it will require some creativity so 1) it will work without us explicitly telling the algorithm beforehand which points are to be sampled and 2) it is callable like an ordinary `evalf`.

[Stefan Krastanov]

There is also `Sum(1/x**constant, (x, 1, t))` plotted for t in [1, 10] that exhibits the same problem.

[Aaron Meurer]

The same problem as Integral sampling or the same problem as
DiracDelta (how is Sum defined for non-integer limits?).

I think in general we should expand out the computation to a symbolic
formula and use cse() to make it more efficient (what's the point of
being symbolic if we can't do cool tricks like this).

Aaron Meurer

[Stefan Krastanov]

Concerning Sum: the same problem as integral (repeated sampling).

About the noninteger limits: it seems to transform sum into an integral, but I guess this is just a coincidence/implementation detail somewhere in sympy and not anything documented.

About cse and other optimizations: Matthew has written some interesting things about using sympy to optimize numerics on his blog.

Slightly offtopic: there is a flag that can be set when plotting so that the lineplot is in a staircase pattern. It is quite useful, for instance, when making plots about approximations of integrals or plots of Sum. It is described in the plot_advanced notebook in the examples folder.

[Aaron Meurer]

On Mon, Jul 29, 2013 at 2:13 PM, Stefan Krastanov
<krastano...@gmail.com> wrote:
> Concerning Sum: the same problem as integral (repeated sampling).
>
> About the noninteger limits: it seems to transform sum into an integral, but
> I guess this is just a coincidence/implementation detail somewhere in sympy
> and not anything documented.

I think this might be part of the Karr convention we recently adopted
(before that, things were inconsistent). The important thing is that
Sum(f(x), (x, 0, 0.5)) + Sum(f(x), (x, 1.5, 2)) == Sum(f(x), (x, 0,
2)) (I hope I got that right). Raoul can give a better explanation.

Aaron Meurer