Re: [sympy] Convolution of piecewise functions

[David Joyner]

On Tue, Feb 5, 2013 at 8:46 AM, Francesco Paparella
<francesco...@unile.it> wrote:
> Is there any way to teach sympy to yield a (useful and correct) result when
> convolving two piecewise functions?


I wrote a python program for Sage to compute the convolution of two piecewise
defined functions. If you want to convert it to Sympy, let me know and I'll
email you a link to it.


> Here is an example (SymPy 0.7.2-git)
>

...

>

[Aaron Meurer]

On Tue, Feb 5, 2013 at 6:46 AM, Francesco Paparella

<francesco...@unile.it> wrote:
> Is there any way to teach sympy to yield a (useful and correct) result when
> convolving two piecewise functions?

> Here is an example (SymPy 0.7.2-git)
>

> --------------------------------------------------------------------------------------------------------------------------------------
> In [1]: f = Piecewise((0, x<=-1), (Rational(1,2), x<1), (0, True))
> In [2]: integrate(f.subs(x, x-y, simultaneous=True)*f(x), y)
> Out[2]:
> ⎧ 0 for x ≤ -1
> ⎪
> ⎪⎧0 for x - y ≤ -1 for x < 1
> ⎪⎪
> ⎪⎪y
> ⎨⎨─ for x - y < 1
> ⎪⎪4
> ⎪⎪
> ⎪⎩0 otherwise
> ⎪
> ⎩ 0 otherwise
> ---------------------------------------------------------------------------------------------------------------------------------------
> Note that I have used f.subs(x, x-y, simultaneous=True) because
> piecewise-defined functions seem to
> have memory of the name of the dummy variable used to define them.

Why do you have to do this? It seems to be the same without simultaneous=True.

Probably this would require something like
http://code.google.com/p/sympy/issues/detail?id=2128. Also, we don't
currently have any kinds of algorithms to simplify the piecewise
conditions.

One idea is that in that issue it is noted that integration of
Heaviside functions works. So maybe you could try that.
Unfortunately. Piecewise.rewrite(Heaviside) is not implemented yet, so
you'll have to do it manually. For example, your function is
Heaviside(x + 1)/2 - Heaviside(x - 1)/2 (except at the points -1 and
1, because of the way Heaviside is defined at 0). Integrating the
convolution of that gives a complex piecewise expression, which is
based on whether |y| is less than or greater than |x + 1| and/or |x -
1|. It's actually more complex than it looks, because it gives some
(probably) singular G-functions in the cases where |y| = |x + 1| or |x
- 1|.

This might actually not be useful to you, since things get complex
quite fast, but you can see if it leads you anywhere.

Aaron Meurer

>
> If I am not mistaken, the correct and useful result should be
>
> (2-abs(x))/4, abs(x)<2
> 0, otherwise
>
>
> --
> 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?hl=en.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Francesco Paparella]

Thanks a lot for the answer.

Il giorno mercoledì 6 febbraio 2013 03:09:38 UTC+1, Aaron Meurer ha scritto:

On Tue, Feb 5, 2013 at 6:46 AM, Francesco Paparella
<francesco...@unile.it> wrote:
> Is there any way to teach sympy to yield a (useful and correct) result when
> convolving two piecewise functions?
> Here is an example (SymPy 0.7.2-git)
>
> --------------------------------------------------------------------------------------------------------------------------------------
> In [1]: f = Piecewise((0, x<=-1), (Rational(1,2), x<1), (0, True))
> In [2]: integrate(f.subs(x, x-y, simultaneous=True)*f(x), y)
> Out[2]:
> ⎧        0           for x ≤ -1
> ⎪
> ⎪⎧0  for x - y ≤ -1  for x < 1
> ⎪⎪
> ⎪⎪y
> ⎨⎨─  for x - y < 1
> ⎪⎪4
> ⎪⎪
> ⎪⎩0    otherwise
> ⎪
> ⎩        0           otherwise
> ---------------------------------------------------------------------------------------------------------------------------------------
> Note that I have used f.subs(x, x-y, simultaneous=True) because
> piecewise-defined functions seem to
> have memory of the name of the dummy variable used to define them.

Why do you have to do this? It seems to be the same without simultaneous=True.

I'm unfamiliar with sympy. I think to understand that 'simultaneous=True' is a true change of variable because it delays any evaluation/simplification to after all substitutions have been done. In any case it is irrelevant for this example.

Probably this would require something like
http://code.google.com/p/sympy/issues/detail?id=2128.  Also, we don't
currently have any kinds of algorithms to simplify the piecewise
conditions.

One idea is that in that issue it is noted that integration of
Heaviside functions works.  So maybe you could try that.
Unfortunately. Piecewise.rewrite(Heaviside) is not implemented yet, so
you'll have to do it manually.  For example, your function is
Heaviside(x + 1)/2 - Heaviside(x - 1)/2 (except at the points -1 and
1, because of the way Heaviside is defined at 0).  Integrating the
convolution of that gives a complex piecewise expression, which is
based on whether |y| is less than or greater than |x + 1| and/or |x -
1|. It's actually more complex than it looks, because it gives some
(probably) singular G-functions in the cases where |y| = |x + 1| or |x
- 1|.

Using Heaviside functions might help, if this yielded clear, useful results. Unfortunately sympy yields something fishy:

In [1]: integrate( Heaviside(x + 1)/2 - Heaviside(x - 1)/2, (x, -10, -5))

Out[1]: -5/2

...I expected to see zero as a result.

 

In the end I think I'll have to write some convolve_piecewise(f,g) function by myself. How do I extract the pieces of a piecewise function from its instance? I tried the obvious:

In [2]: q=Piecewise((0,x<0),(x, x<1),(0,True))

In [3]: for i in q:

    pprint(i)

expecting to see the list of the (expression, condition) pairs, but it doesn't work: a piecewise function is not iterable.

[Francesco Paparella]

[Francesco Paparella]

I think I'll have to write some code to do the convolution piece by piece. Therefore any help is welcome. Thanks.

[Aaron Meurer]

On Thu, Feb 7, 2013 at 2:15 AM, Francesco Paparella

simultaneous=True doesn't even make sense if you are only making one
substitution. It's for when you want to do some kind of "swap". For
example:

>>> (x*y).subs([(x, y), (y, z)])
z**2
>>> (x*y).subs([(x, y), (y, z)], simultaneous=True)
y*z

Note that if I had used a dictionary in the first one, the result
wouldn't even be well-defined, because it would depend on the order of
iteration from the dictionary.

>>
>>
>> Probably this would require something like
>> http://code.google.com/p/sympy/issues/detail?id=2128. Also, we don't
>> currently have any kinds of algorithms to simplify the piecewise
>> conditions.
>>
>> One idea is that in that issue it is noted that integration of
>> Heaviside functions works. So maybe you could try that.
>> Unfortunately. Piecewise.rewrite(Heaviside) is not implemented yet, so
>> you'll have to do it manually. For example, your function is
>> Heaviside(x + 1)/2 - Heaviside(x - 1)/2 (except at the points -1 and
>> 1, because of the way Heaviside is defined at 0). Integrating the
>> convolution of that gives a complex piecewise expression, which is
>> based on whether |y| is less than or greater than |x + 1| and/or |x -
>> 1|. It's actually more complex than it looks, because it gives some
>> (probably) singular G-functions in the cases where |y| = |x + 1| or |x
>> - 1|.
>>
> Using Heaviside functions might help, if this yielded clear, useful results.
> Unfortunately sympy yields something fishy:
> In [1]: integrate( Heaviside(x + 1)/2 - Heaviside(x - 1)/2, (x, -10, -5))
> Out[1]: -5/2
> ...I expected to see zero as a result.

That's a bug. Can you report it at http://code.google.com/p/sympy/issues/list?

>
> In the end I think I'll have to write some convolve_piecewise(f,g) function
> by myself. How do I extract the pieces of a piecewise function from its
> instance? I tried the obvious:
> In [2]: q=Piecewise((0,x<0),(x, x<1),(0,True))
> In [3]: for i in q:
> pprint(i)

To extract the parts of any SymPy expression, use .args. For example

>>> Piecewise((0,x<0),(x, x<1),(0,True)).args
((0, x < 0), (x, x < 1), (0, True))

Aaron Meurer

>
> expecting to see the list of the (expression, condition) pairs, but it
> doesn't work: a piecewise function is not iterable.
>