Articles

Status Update—Week 5

In Uncategorized on June 24, 2011 by nessgrh

This was one busy week. Here is a list of things I did:

  • Start inverse mellin transform.
  • Implement Laplace transform, start inverse Laplace transform.
  • Implement expansion of G functions of finite confluence.
  • Implement the recursive mellin transform method.
  • Re-organise the integration code, Performance improvements.

Let’s look at everything in turn.

Inverse Mellin Transform

The inverse mellin transform is a cornerstone of the recursive mellin transform method. Since the meijer G function is essentially defined as the inverse mellin transform of a fairly general function, this is “only” about recognising a well-defined set of building blocks (rational functions, exponentials, trigonometric functions, gamma function), rewriting them using gamma functions, and turning this information into a G function. As usual the devil is in the details, but the code works quite well by now. Here are a few examples. Note that inverse mellin transforms are only defined relative to a fundamental strip (i.e. a strip parallel to the imaginary axis, free of singularities): they are defined as a contour integral along a line in the strip. This is what is being passed as the last parameter. But notice that the fundamental strip can be recovered unambiguously from only one of its edges. Namely if we say know the left edge, then any contour sufficiently close to the right of it will give the right answer—indeed any contour to the left of all singularities (since we can only work with a small set of functions anyway, finding the singularities is easy). This can be handy sometimes, and is triggered by passing None for one of the edges (see last example).

Currently unimplemented are inverse mellin transforms of trigonometric functions, and some transformations on gamma functions (Basically for the meijer G function, all gamma functions must have argument \pm s + c. If e.g. it is an integer multiple of s, we can use the gamma function multiplication theorem.)

Laplace Transforms

The laplace transform of a function is defined as F(s) = \int_0^\infty e^{-st} f(t) \mathrm{d}t. This is evidently an integral of the form we can handle. Notice also that the laplace transform is related by a change of variables to the mellin transform. But computing laplace transforms via mellin transforms is unlikely to work since the functions become much more complicated (multiply by a step function, and introduce a logarithm in the innermost part). Hence the laplace_transform function directly tries to compute the integral. Here are some examples (the very first example used to be quite slow—see the section below). Similarly to mellin transforms, laplace transforms of all “sensible” functions converge for Re(s) > a, i.e. in a half plane. Results are returned as (transform, a, cond) triples where `a` designates the half plane of convergence, and cond are auxiliary convergence conditions (these are always True here, because of the assumptions we set up on a, b, c, d).

The inverse laplace transform consists of a Mellin-type contour integral. Again it can be expressed via inverse mellin transforms, but this time chances that this works are much higher: the change of variables effectively happens after the integration, so the integrand is not made more complicated. This is what I currently implemented, and it works for many interesting examples. Notice that the inverse laplace transform is also defined (in principle) relative to a half plane. But here we can use the capabilities of inverse_mellin_transform to pass the left edge of the “strip” (a, \infty) as None. This is why the plane does not have to be specified.

However there are common examples that are unlikely to ever work this way (e.g. inverse laplace transform of laplace transform bessel functions, or even inverse laplace transform of s^a for symbolic/non-integer a). But it turns out that the laplace inversion integral of a G function is again expressible as a G function. This yields another route to computation, but I have not yet implemented it.

G functions of finite confluence

Slater’s theorem can be used to expand a G function in terms of hypergeometric functions if all poles of the integrand in the definition are simple. If there are multiple poles, this does not work; the resulting series are not hypergeometric. We hit such a case when computing the inverse mellin (or laplace) transform of e.g. 1/s^2.

Even though it is not expressible (directly) using hypergeometric functions only, if there are only finitely many multiple poles, then we can write the G function as a finite sum of residues, plus some hypergeometric series (I didn’t find this worked out in the literature, but its a fairly simple generalisation of Slater’s theorem. Just messier. I wrote it down in detail in a documentation file.) Thus we need to compute residues of gamma products, for which we need laurent expansions of gamma functions. I added these as well. This all works now. The last answer comes out incredibly messy without the simpilfy—a bug.

Recursive Mellin Transforms; Optimisations

I described the rescursive mellin transform method in my last blog post. This is now implemented, allowing us to compute e.g. this.

However, this recursive code can be fairly slow. I sped it up a little, using a simple heuristic and caching. For the heuristic, recall that we try to write the integrand as a product of two G functions. The heuristic is to try to make only one power of x occur in every term that is supposed to be a G function. E.g. if we want to rewrite e^{x^2} e^x e^{3x^2} e^{2x}, this code will first try to rewrite this as \left(e^{x^2} e^{3x^2} \right) \left( e^x e^{2x} \right).

Also I made the use of dummies in my code more caching-friendly (not creating too many new ones), and declared a key function @cacheit. The algorithm implementation could be changed so as to avoid redundant calls, but this way is much easier. Computing above integral now takes about 500ms, mostly due to the heuristic. The above-mentioned laplace transform now takes about 1.7 seconds, mostly due to caching. These times are not great, but I think they are ok for now. I don’t want to hit the “premature optimisation trap”.

Conclusion

That’s it for this week. Next week I will add the laplace inversion integral, and Fourier transforms. I also have a growing list of assorted bugs and todos. With some luck, by the end of next week I will be able to start extending the lookup tables.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: