Tricky branch cuts

Articles

Tricky branch cuts

In Uncategorized on July 7, 2011 by nessgrh

This post will elaborate on the recent problems I’m having with “multivalued” functions. My usual status update will come tomorrow. I shall try to be almost pedantically clear in my description, both in order to clear my own thoughts and to increase the chance of someone suggesting a way out. In what follows, $\mathcal{S}$ will denote the riemann surface of the logarithm (this only means infinitely many copies of the complex line, glued together approprietely so that $\log{z}$ becomes an entire, single-valued function on $\mathcal{S}$ ). Zero is not an element of $\mathcal{S}$ .

The problem we will study is how to compute the inverse fourier transform $f(t) = \int_{-\infty}^\infty e^{i t x} \frac{\mathrm{d}x}{i x + 1}$ . First of all it is easy to see that the integral converges only for real, non-zero $t$ . Also using contour integration it is easily established that $f(t) = 2\pi H(t) e^{-t}$ , where $H(t)$ denotes the Heaviside step function.

How does the Meijer G-function algorithm compute this integral? First of all it is split up into two parts, $I_1 = \int_0^\infty e^{i t x} \frac{\mathrm{d}x}{i x + 1}$ and $I_2 = \int_{-\infty}^0 e^{i t x} \frac{\mathrm{d}x}{i x + 1}$ . Let’s concentrate on $I_1$ for now. The integral exists if the exponential has a decaying component, or if it is oscillatory. That is to say it converges for $Im(t) \ge 0, t \ne 0$ . But let’s be a bit more specific: all meijer G functions are defined on $\mathcal{S}$ , and may or may not descend to $\mathbb{C}$ . Thus define $I_1: \mathcal{S} \supset D \rightarrow \mathbb{C}$ to be our integral (where $D$ denotes the region of convergence). Of course $I_1(z) = I_1(e^{2\pi i} z)$ when the integral converges (note that this equation is not vacuous: $e^{2\pi i} \ne 1$ on $\mathcal{S}$ ). We thus find:

The integral $I_1(t)$ converges if $(2n-1) \pi \le Arg(t) \le 2n \pi$ , for some integer $n$ .

Having established this, let’s see what the algorithm computes. I should explain the output. Note that debugging output is enabled. The first piecewise result is what would ordinarily be presented to the user: there is a closed expression, holding under various conditions. And there is the “otherwise” part of the piecewise, which just says that if the conditions do not hold, we don’t know the answer. Consider the conditions. unbranched_argument(t) means the argument of $t \in \mathcal{S}$ (whereas arg(t) would mean some value between $-\pi$ and $\pi$ , and that function is not continuous either). Thus the conditions are just a funny way of saying $0 \le Arg(t) \le \pi$ . Now consider the claimed result. It looks a bit daunting, but this is because the functions in which it is written are defined on $\mathbb{C}$ , not $\mathcal{S}$ . The debugging output line “Result before branch substitutions …” shows what the result would be on $\mathcal{S}$ . I have computed this in output [8] (notice that what looks like $\pi$ here is a dummy, to stop evaluation of $e^{\pi i}$ etc). Of course $\exp(z)$ is an unbranched function on $\mathbb{C}$ , so the first two factors are uninteresting. But the last factor is not: we see that the result is related to the upper incomplete gamma function. This is a branched function, i.e. it is properly defined on $\mathcal{S}$ and it does not descend to $\mathbb{C}$ . If $t$ is a positive real number, then the incomplete gamma function is to be evaluated at argument $-\pi$ . This is where the extra terms (like argument_period()) come from: on $\mathbb{C}$ argument $-\pi$ means negative reals. But with the branch cut conventions as implemented (in mpmath, effectively, but also enforced throughout sympy), negative reals mean argument $+\pi$ . So the answer cannot just be $\Gamma(0, -t)$ , because this would end up on the wrong branch. Thus we conclude (or rather the algorithm did for us):

Define $g: \mathcal{S} \rightarrow \mathbb{C}$ by $g(z) = -i e^{-t} \Gamma(0, e^{-i \pi} t)$ . Then for $0 \le Arg(t) \le \pi, I_1(t) = g(t)$ .

Notice in particular that our original integral has two natural extensions to $\mathcal{S}$ : the naive definition $I_1(t)$ , and the analytic continuation $g(t)$ . The algorithm has to return some subset of the first definition, and not the second. In a sense this is why the convergence condition is necessary (although technically it is there because the integral representation of the G-functions used to establish this formula diverges outside the region specified by the convergence conditions).

We could now do a similar analysis for $I_2$ . The story is essentially the same, with the upper half plane replaced by the lower half plane ( $-\pi \le Arg(t) \le 0$ ) throughout. So when can we compute the combined integral? We already noted that this should be the case exactly when $t$ is real, and working in $\mathbb{C}$ this is easy to see: the intersection of the (closed) upper and lower half planes is the real axis. But notice how this fails on $\mathcal{S}.$ The intersection of $0 \le Arg(t) \le \pi$ (where we can evaluate the first integral) and $-\pi \le Arg(t) \le(t)$ is $Arg(t) = 0$ , i.e. only positive $t$ .

To further amplify the confusion, computing said fourier transform with argument $-t$ , where $t$ is declared positive, works. This is because before “switching” to working on $\mathcal{S}$ , the algorithm carelessly cancels some minus signs … which is fine for entire functions, of course.

Here ends this rather lengthy explanation. It remains to figure out how we can teach the algorithm that the integrals are periodic in argument $2\pi$ , instead of knowing them only on one half-sheet. I have some vague ideas for this, but input is appreciated. In any case, I’ll think about this for a while. For now there is enough other work to do.

2 Comments

2 Responses to “Tricky branch cuts”

asmeurer July 8, 2011 at 00:35 | Reply

This may or may not help, but it might be instructive to temporarily disable automatic simplification of exp(n*I*pi). Just add return to line 58 of sympy/functions/elementary/exponential.py. I get:


In [8]: t = Symbol('t', positive=True)

In [9]: exp(-I*pi)
Out[9]: 
 -ⅈ⋅π
ℯ    

In [10]: inverse_fourier_transform(1/(1+2*pi*I*x), x, exp(-pi*I)*t, noconds=False)
Out[10]: 
                       ⎛     1             -ⅈ⋅π⎞
InverseFourierTransform⎜───────────, x, t⋅ℯ    ⎟
                       ⎝2⋅ⅈ⋅π⋅x + 1            ⎠

In [11]: inverse_fourier_transform(1/(1+2*pi*I*x), x, t, noconds=False)
Out[11]: 
⎛                                                                       -ⅈ⋅π      ⎞
⎜⎛    -ⅈ⋅π                          -ⅈ⋅π       ⅈ⋅π                   ⎞  ────      ⎟
⎜⎜ t⋅ℯ      ⎛      ⅈ⋅π⎞          t⋅ℯ        t⋅ℯ    + ⅈ⋅π  ⎛      ⅈ⋅π⎞⎟   2        ⎟
⎜⎝ℯ       ⋅Γ⎝0, t⋅ℯ   ⎠ + 2⋅ⅈ⋅π⋅ℯ        + ℯ            ⋅Γ⎝0, t⋅ℯ   ⎠⎠⋅ℯ          ⎟
⎜───────────────────────────────────────────────────────────────────────────, True⎟
⎝                                    2⋅π                                          ⎠

In [12]: inverse_fourier_transform(1/(1+2*pi*I*x), x, t, noconds=False)[0].subs(t, exp(-pi*I)*t)
Out[12]: 
                                                            -ⅈ⋅π
⎛    -2⋅ⅈ⋅π                     -2⋅ⅈ⋅π                   ⎞  ────
⎜ t⋅ℯ                        t⋅ℯ          t + ⅈ⋅π        ⎟   2  
⎝ℯ         ⋅Γ(0, t) + 2⋅ⅈ⋅π⋅ℯ          + ℯ       ⋅Γ(0, t)⎠⋅ℯ    
────────────────────────────────────────────────────────────────
                              2⋅π

I don’t even know if that’s right or not. If it’s wrong, it’s could just be because some part of your code relies on the automatic simplification for correctness.

nessgrh July 8, 2011 at 02:27 | Reply

This is exactly the expected behaviour (you can either disable the simplification, or simulate everything using dummies). The first transform not being evaluated is correct. [I did not explain this above, but if you use exp(-I*pi)*t, then the *first* integral cannot be done. The thing is that the first has do be done with argument pi, and the second with argument -pi—no matter which one we pass, one of them will not work.]

The second equation is again the analytic continuation and not what we want.

I think I know a simple sufficient condition for a G-function to be entire, which I believe is also fairly close to necessary. It does apply to all functions currently in the table (which is admittedly not many). Building this in should make these examples work.

But I will instead extend the tables tomorrow (finally), and sleep about this for a few more days.

nessgrh

About

Hello, World

Meta