## Status Update—Week 4

In Uncategorized on June 17, 2011 by nessgrh

## Introduction

This week I finally started writing code to compute integrals. Let’s quickly recap the general strategy employed: given an integral $\int_I f(x) \mathrm{d}x$, rewrite it in one of the following forms:

1. $\int x^s G(x) \mathrm{d}x$
2. $\int_0^\infty x^s G(ax) \mathrm{d}x$
3. $\int_0^\infty x^s G(ax) G(bx) \mathrm{d}x$

Here $G(x)$ represents any G function. If this is possible, we have a chance of evaluating the integral in terms of more G functions, and the code from the previous weeks takes over to expand these in terms of known functions (hopefully). Note that in cases (2) and (3) there are of course very many conditions under which the integrals may or may not converge. These have been tabulated extensively, and this way we can, given an integrand with parameters, actually determine convergence regions! My current code contains more than 30 conditions for (3) which seem to take care of almost all cases.

There are two parts to the re-writing: (a) converting the limits into the form $0, \infty$, and (b) rewriting the integrand using G-functions. Both seem to be more of an art than a science. Much of my efforts during the rest of GSOC will probably involve improving the code for (a) and (b) in some way. For the curious, (b) is essentially equivalent to computing the Mellin transform of a function. This of course again involves integration, and opens the door for a recursive algorithm. I have not implemented that, yet.

## Examples

Let’s take a look at some examples. Another GSOC student needs to compute integrals of Gaussian functions. This works to some extent. There are several interesting things to notice. First of all, declaring $\sigma, \mu$ to be positive is crucial, otherwise the code would return a piecewise expression with conditions of convergence. Second, we see in (2, 3) that general gaussian integrals can be computed. What happens is that the code expands the integrand and then ends up with a product of two exponentials, both of which are easily rewritten as G-functions. On the other hand, (4) does not work: we end up with three exponentials, and we have no formula for that. Of course the two with linear  exponent could be combined, but the algorithm does not recognise that. Here the recursive mellin transform technique would have been helpful. (Also note the new parameter meijerg for integrate(): it can be either True, False, or None (default), which mean, respectively, to always use meijer G function methods, never use meijer G function methods, and to use all methods available.)

That brings us to another very important application: integral transforms. It is clear that the code should be able to compute any common integral transform (Fourier, Laplace, Mellin, etc) of any G function, which means integral transforms of *lots* of functions. I will add code that wraps integrate() to apply various conventions of integral transforms, e.g. to provide an analytic continuation of the laplace transform, and the critical strip of a mellin transform, etc. For now there is only mellin_transform. Given any function $f(x)$, it computes the Mellin Transform $F(s)$, the critical strip $(a, b)$ (such that the defining integral converges for $a < Re(s) < b$), and possibly auxiliary convergence conditions. Here are a few examples. Again there are a few noteworthy things. First of all we see that the convergence conditions returned by the code are strong enough to compute the fundamental strips—those shown indeed agree with what can be found in the literature. Secondly, we see that sometimes the result is expressed using factorials, where gamma functions would be more appropriate. This comes from the function combsimp(), which I will have to improve. Thirdly, we see that all mellin transforms shown are essentially products exponentials and gamma functions involving $s$. This is no surprise, since this is exactly the condition for being expressible as a G function. In particular, example (3) demonstrates that using recursive mellin transforms, the non-working example (4) from the previous link could be made to work. Finally, there is a catch: the algorithm could not determine any convergence conditions for positive real t. This probably (hopefully) means that I mis-typed some of the convergence conditions—I will have to investigate this.

## Conclusion

So much for this week. Next week I will concentrate on improving all the problems pointed out above (that’s of course quite open-ended…), and on extending the tables used to rewrite functions as G functions. This should hopefully yield many non-trivial integrals.

### 3 Responses to “Status Update—Week 4”

1. Don’t forget that there are at least three other transforms from Mellin, Fourier, and Laplace, namely the inverse of those three. Will your code be able to handle those? I’m particularly wondering about the inverse Mellin, whose integration limits are c – I*oo and c + I*oo.

There’s also a bunch of other integral transforms that you might as well implement if they work.

2. Inverse fourier should be fine (of course). Other transforms (e.g. Hankel transform, fourier sine transform) should be possible, I just didn”t mention them because I don’t really know anything about them.

Inverse Mellin and Laplace transforms are more worrying, and I don’t immediately see what to do about them. Note that the “c” actually depends on the function to be inverted (although when computing the direct mellin or laplace transform, we find what the c has to be – it is just not usually specified).

Computing the inverse mellin transform of a quotient of gamma functions and exponentials is of course easy, it is just a G function. But the general case seems very hard.

3. That reminds me that there are different conventions for the Fourier transforms relating to the coefficient (some put $\frac{1}{2\pi}$ in front of one and nothing in front of the inverse, while some put $\frac{1}{\sqrt{2\pi}}$ in front of both). I suppose we should just pick a convention and stick with it.