Issues with integrate and piecewise functions

[bonc...@udel.edu]

I'm trying to write a generic function to convolve two functions (in one dimension, for now) and came across a few issues I don't know how to solve.  Rather than give up, maybe I can start a discussion.  Some Sympy follows:

a , b = symbols('a b', positive=True)
f = exp(-a*x)
g = exp(-b*x)
integrate(f.subs(x,tau)*g.subs(x,t-tau), (tau,0,t))

1. I had to compute the range of integration myself (tau, 0, t).  This isn't a fatal problem, but piecewise doesn't have enough information to compute it simply.  (My long term solution is to propose a replacement for the piecewise function, but I'm not there yet.)

The answer is below:

Piecewise((t*exp(-b*t), Eq(a, b)), (1/(a*exp(b*t) - b*exp(b*t)) - 1/(a*exp(a*t) - b*exp(a*t)), True))

The answer is correct, but there are problems trying to incorporate it into a larger function.

2. It is not obvious in advance the answer will depend on whether a = b.  Is there some way to assume in advance that a=b or a!=b?

3. The second condition, True, isn't helpful.  To understand what True means, one has to keep track of the previous conditions.  Is there some way to replace the True with a != b.

4. This is a simpler question, but I can't figure it out.  How do I manipulate the second answer to my preferred form?

(exp(-b*t)-exp(-a*t))/(a-b)

Is there a simple way to do it programmatically?

Thanks

[Aaron Meurer]

On Wed, Jun 7, 2017 at 11:35 AM, <bonc...@udel.edu> wrote:
> I'm trying to write a generic function to convolve two functions (in one
> dimension, for now) and came across a few issues I don't know how to solve.
> Rather than give up, maybe I can start a discussion. Some Sympy follows:
>
> a , b = symbols('a b', positive=True)
> f = exp(-a*x)
> g = exp(-b*x)
> integrate(f.subs(x,tau)*g.subs(x,t-tau), (tau,0,t))
>
>
> 1. I had to compute the range of integration myself (tau, 0, t). This isn't
> a fatal problem, but piecewise doesn't have enough information to compute it
> simply. (My long term solution is to propose a replacement for the
> piecewise function, but I'm not there yet.)
>
> The answer is below:
>
> Piecewise((t*exp(-b*t), Eq(a, b)), (1/(a*exp(b*t) - b*exp(b*t)) -
> 1/(a*exp(a*t) - b*exp(a*t)), True))
>
>
> The answer is correct, but there are problems trying to incorporate it into
> a larger function.
>
> 2. It is not obvious in advance the answer will depend on whether a = b. Is
> there some way to assume in advance that a=b or a!=b?

The assumptions don't work with this yet, but the simplest way is to just do

expr = integrate(f.subs(x,tau)*g.subs(x,t-tau), (tau,0,t))
expr.subs(Eq(a, b), True)

>
> 3. The second condition, True, isn't helpful. To understand what True
> means, one has to keep track of the previous conditions. Is there some way
> to replace the True with a != b.

I don't think there's a simple way to do it, but it can be represented
that way. It might be useful to have a Piecewise method that replaces
the True condition with the negation of the other conditions.

>
> 4. This is a simpler question, but I can't figure it out. How do I
> manipulate the second answer to my preferred form?

This is harder, because it's hard to get SymPy to prefer exp(-a) over
exp(a). Using simplify gives

(exp(a*t) - exp(b*t))*exp(-t*(a + b))/(a - b)

There is a related bug here
https://github.com/sympy/sympy/issues/11506. If it were fixed you
could use cancel(expr, exp(-a*t), exp(-b*t)).

The only way I found to do it is

simplify(expr.subs({a: -a, b: -b})).subs({a: -a, b: -b})

Aaron Meurer

>
> (exp(-b*t)-exp(-a*t))/(a-b)
>
> Is there a simple way to do it programmatically?
>
> Thanks
>

> --
> 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 https://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/9d71a081-e038-4712-93de-aca17f951e15%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[bonc...@udel.edu]

BTW, when I'm demonstrating sympy to people, I show them we can test our preferred form for equality.  (Sympy is handy for verifying hand calculations or standard formulas found in books.)

h = integrate(f.subs(x,tau)*g.subs(x,t-tau), (tau,0,t))
simplify(h.args[1][0] - (exp(-b*t)-exp(-a*t))/(a-b))

0

Your last simplify trick is impressive.  It works, but isn't obvious.

My thoughts on Piecewise is that it's trying to do too many things.  There's a difference between representing a function whose form changes depending on the value of its argument (in my case, t) and an answer that depends on parameter values (a and b).  The latter distinction is best described with a cases statement.  I.e., t is different from a and b.

In signal processing (my area) we often have signals whose form changes as t changes.  Convolving them must be done piece by piece and the pieces assembled into a whole.  (The assembling difficulty is that each convolved piece has a bigger domain than either of the input pieces--but this can be managed.)  My roadblock is when the integration depends on parameter values (as above).

Charlie

[Aaron Meurer]

On Wed, Jun 7, 2017 at 3:12 PM, <bonc...@udel.edu> wrote:
> BTW, when I'm demonstrating sympy to people, I show them we can test our
> preferred form for equality. (Sympy is handy for verifying hand
> calculations or standard formulas found in books.)
>
> h = integrate(f.subs(x,tau)*g.subs(x,t-tau), (tau,0,t))
> simplify(h.args[1][0] - (exp(-b*t)-exp(-a*t))/(a-b))
>
> 0
>
> Your last simplify trick is impressive. It works, but isn't obvious.

It works because SymPy wants to treat exp(-a*t) as 1/exp(a*t), and
simplify doesn't want to have nested fractions. So I force it into a
form that isn't a fraction, call simplify(), then put it back.

>
> My thoughts on Piecewise is that it's trying to do too many things. There's
> a difference between representing a function whose form changes depending on
> the value of its argument (in my case, t) and an answer that depends on
> parameter values (a and b). The latter distinction is best described with a
> cases statement. I.e., t is different from a and b.

I'm not sure I follow here. Are you suggesting to use a separate
object. What would it look like?

>
> In signal processing (my area) we often have signals whose form changes as t
> changes. Convolving them must be done piece by piece and the pieces
> assembled into a whole. (The assembling difficulty is that each convolved
> piece has a bigger domain than either of the input pieces--but this can be
> managed.) My roadblock is when the integration depends on parameter values
> (as above).

I believe this old issue about integrating piecewise expressions where
the variable is in the conditions is still relevant
https://github.com/sympy/sympy/issues/5227#issuecomment-36998678.

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/e0836c94-efe2-429d-ace6-c5385b181222%40googlegroups.com.

[bonc...@udel.edu]

On Wednesday, June 7, 2017 at 5:49:34 PM UTC-4, Aaron Meurer wrote:

>
> My thoughts on Piecewise is that it's trying to do too many things.  There's
> a difference between representing a function whose form changes depending on
> the value of its argument (in my case, t) and an answer that depends on
> parameter values (a and b).  The latter distinction is best described with a
> cases statement.  I.e., t is different from a and b.

I'm not sure I follow here. Are you suggesting to use a separate
object. What would it look like?

I'm thinking of a "range" object (I wish the name "Piecewise" was still available). Something like this:

f = range_function([(f1,t1,t2),(f2,t3,t4),...],t)

The first argument is a list of functions with the range of the arguments for that function. The final argument is the variable.  The notation mimics the integrate function. 

I've written a partially completed convolution function for these range functions.  Convolution requires knowing the range of each function.  But I got stuck when I ran the example before whose answer depends on parameter values.  That the answer depends on parameter values is the job for a "cases" function (like the Latex cases function).

The last step (which I haven't finished, but it's straightforward) is to combine the outputs of the various convolutions into a range function.

I haven't given much thought how to extend the range function to multiple variables.  It's probably straightforward when the ranges are rectangles, but maybe less for other shapes (e.g., circles).

 

[Chris Smith]

In comment 1 you said you had to compute the range yourself. What do you mean and why did you have to do so?

Regarding a later comment about handing convolutions where different functions have different ranges, that *should* be handled with Piecewise. Perhaps you could give my `lts` branch in PR #12587 a try. A major overhaul of Piecewise is presented there.

[Chris Smith]

Here is an example of integrating a function with conditions independent of the integration variable:

>>> A, B, a, b x = var('A B a b x')
>>> Piecewise((A, And(x<0, a<1)), (B, Or(x<1, a>2))).integrate(x)
Piecewise(
    (B*x, a > 2),
    (Piecewise((A*x, x <= 0), (B*x, x <= 1), (nan, True)), a < 1),
    (Piecewise((B*x, x <= 1), (nan, True)), True))

[Chris Smith]

Here is an example from the `lts` branch computing the product of two functions

that have different ranges in each:

>>> p = Piecewise((a,abs(x-1)<1),(b,abs(x-2)<2),(c,True))*Piecewise((d,x>1),(e,True))
>>> piecewise_fold(p)
Piecewise(
    (a*d, (x > 1) & (x < 2)),
    (b*d, (x > 1) & (x < 4)),
    (c*d, x > 1),
    (a*e, Abs(x - 1) < 1),
    (b*e, Abs(x - 2) < 2),
    (c*e, True))
>>> _.integrate(x).diff(x)
Piecewise(
    (c*e, x <= 0),
    (a*e, x <= 1),
    (a*d, x <= 2),
    (b*d, x <= 4),
    (c*d, True))

[Note that value `b*e` does not appear in the final result because it has been determined to not appear because higher priority expression are defined on the range in which `b*e` is defined.]

I like the idea of range_function. Here are two ways it might be done:

>>> def range_function(x, *args):
...  return Piecewise(*[
...   (a[0], And(x>=a[1],x<=a[2])) for a in args])
...
>>> range_function(x, (a,1,2),(b,3,4),(c,0,6))
Piecewise((a, (x >= 1) & (x <= 2)), (b, (x >= 3) & (x <= 4)), (c, (x >= 0) & (x <= 6)))


>>> def interval_function(x, *args):
...  return Piecewise(*[
...   (a[0], a[1].contains(x)) for a in args])
...
>>> interval_function(x, (a,Interval(1,2)),(b,Interval.Lopen(0,4)))
Piecewise((a, (x >= 1) & (x <= 2)), (b, (x <= 4) & (x > 0)))