Hi.
On Mon, Jul 4, 2011 at 5:49 PM, Tom Bachmann <[email protected]> wrote:
> Hi,
>
> I hit a bit of a set-back today with my integration code. Specifically, I am
> getting wrong results, because of improper treatment of branch cuts.
>
> As a first example consider
>
> sin(sqrt(z)) - sin(sqrt(w)),
>
> where we want to substitute z=exp(I*pi)=-1 and w=exp(-I*pi)=-1. If we do
> this naively, the answer comes out as zero, but what we (at least I, in my
> integration code) mean is 2*sin(I).
>
> I have for now been alleviating this problem (since I started the
> integration code) by explicitely tracking the branch we are on by
> (essentially) substituting early z=|z|*exp(I*arg(z)/pi*my_pi), where my_pi
> is a dummy that is not recognised by any of the simplifications code. When
> doing this, the above expression gets changed to
>
> sin(exp(I*my_pi/2)) - sin(exp(-I*my_pi/2))
>
> and in the end we can safely substitute pi for my_pi and get a correctly
> computed branch.
>
>
> However, this does not work for all functions. Generically, this will fail
> when we replace sin(sqrt(z)) above by a non-elementary function (i.e. one
> that is not recognised by the simplification machinery) branched at the
> origin. A particular example I have hit is uppergamma(0, z). As you may
> know, the function uppergamma(s, z) (for fixed s) is an entire function when
> s is a positive integer, satisfies uppergamma(s, z) = z**s * F(s, z) for an
> entire function F if s is not an integer, and satisfies uppergamma(s,
> exp(2*pi*I*n)*z) = uppergamma(s, z) + 2*pi*n/factorial(-s) for non-positive
> integral s.
Which simplification machinery are you referring to?
>
> My question thus is how to put this information into code. Notice that all
> functions in question come from meijer g-functions, but I don't think there
> is a chance of understanding the monodromy action in this generality (but I
> would love to be wrong!). Assuming this is true, the best I could come up
> with would be to add ._eval_rewrite_as_{un,}branched functions to special
> functions, perhaps even with a default implementation raising
> NotImplementedError. This would be somewhat akin to the
> ._eval_rewrite_as{in,}tractable used for gruntz' algorithm.
>
> [Since the answer is likely messy, any code should obviously try to avoid
> calling these functions at any cost, but that is something I can take care
> of locally.]
>
> In more detail, somewhere in functions/elementary/complexes.py I would add
> principal_branch(z) which represents z with its argument converted to range
> (-pi, pi], and argument_period(z) which then satisfies z =
> principal_branch(z)*exp(I*2*pi*argument_period(z)). Then e.g.
> uppergamma._rewrite_as_unbranched could return an expression involving these
> functions. The integration code would use a sequence of operations as
> follows to compute the branch of the result:
>
> # TODO try to avoid this when possible
> res = res.rewrite('unbranched')
> res = res.subs(my_pi, pi)
> res = res.rewrite('branched')
>
> My real questions are then as follows, I guess:
>
> 1. Does anyone see a better way around the issue?
> 2. Any suggestions for a better interface?
>
> Thanks,
> Tom
>
I don't see any problems with it, but complex analysis is not my
strength, as I haven't take the course yet (I will have it next
semester), so I only know what I've picked up. I think the rewrite
machinery should be rewritten to be more powerful, but this is likely
beyond the scope of your project.
Aaron Meurer
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
