Status Update—Week 6

In Uncategorized on July 2, 2011 by nessgrh

This week saw less spectacular addition of functionality, and more fixing of details. This is what I did:

  • Laplace Inversion integral for G functions
  • Fourier Transforms
  • Simplification of Gamma functions
  • Miscellaneous extensions and fixes.

Laplace Inversion integral for G functions

The inverse laplae tranform can be defined as f(t) = \int_{c - i\infty}^{c + i\infty} e^{zt} F(z) \mathrm{d}z. As discussed last week, this is related to the mellin inversion formula, and indeed mellin inversion can be used to compute many inverse laplace transforms. On the other hand, if F is expressible as a G function of argument az^b for b rational, it turns out that the answer itself is expressible as a G function. This opens another route to computation. However, I could only find a formula in the literature for b = \pm 1. For this reason I worked out in some detail the formula, and, more importantly, the accompanying convergence conditions in a documentation file. As expected, this allows some more computations.

Fourier Transforms

The last big transform to take care of is the Fourier transform. This first involved extending my integration code to handle integrals from -\infty to \infty; it does this by splitting up the integral in (hopefully) intelligent way into two parts. Many of the basic fourier transforms work already, but some do not. This will require some more work. Here are a few examples. As can be seen in particular inverse fourier transforms tend not to work (where I omitted them); this is because of some bugs in the way the integral is split, and also because some integrands are just not recognised (e.g. I should add code to recognise general rational functions). However, it is very pleasing to see that many non-trivial transforms are computed correctly essentially from the integral definition, almost without intervention.

Simplification of Gamma functions

Another fun thing I did this week was to extend the combsimp() function to handle two gamma function identities, namely the multiplication theorem and the reflection formula. These are used to reduce the number of gamma functions as far as possible. The algorithm works quite well now, bit it is getting somewhat slow. Here are a few examples. Indeed the last one is evidently a bug…

Miscellaneous extensions and fixes

Now that the release is out, Aaron has started reviewing my gsoc pull request in detail. This stirs up lots of little bugs/inconsistencies I have been working on. Also my work on the integration code introduced some test failures here and there, which I also have been fixing today. Furthermore, I sometimes hit bugs in sympy and have to fix them along the way, for example the implementations of is_irrational() in Add and Mul are broken (one would not expect this to be much trouble, but actually from this the assumptions engine can infer in a number of steps that certain things are e.g. positive which most definitely are not; for example abs(I*pi) used to return I*pi for me…). Finally I implemented a few simplifications of bessel functions of integral orders, using e.g. symmetry J_n(-x) = (-1)^n J_n(x) etc.


7 Responses to “Status Update—Week 6”

  1. Cool. Maybe inverse_laplace_transform() should call expand_func() automatically, so you only get Bessel functions when necessary.

    How is the last combsimp example a bug? Is it wrong?

  2. What about symbolic/unevaluated transforms?

    Just as summation returns a symbolic sum sign
    if it can not find a closed from, the Fourier/Laplace/*
    transformations could return a symbolic Transformation
    object. (From which we can extract the transformation
    integral if desired by an instance method.)

    Maybe this is better that raising an exception?

    Just some thoughts about this topic …

  3. Aaron: what do you mean “only return bessel when necessary”? The bessel function essentially comes from calling hyperexpand() on a G function, and will return cos/sin for order 1/2 (unless there is a bug).

    The last combsimp is right, but actually it can be simplified further (namely the way I set it up it should use gamm(x-1/4)*gamma(x)*gamma(x+1/4)*gamma(x+1/2) = gamma(4x-1)*C, but instead it applies the duplication case, leaving three gammas in the end instead of one.

    raoul: yes, I thought about this. For now I am concentrating on transforms we can actually do, so the exception might be seen as a temporary thing (or rather “I have no opinion on this, yet”).

  4. I played around with your latest branch and the results are impressive!
    Just wanted to say that 🙂

    The idea about symbolic transformations is of course only the finishing bit
    after the hard work is done. I assume writing a Fourier class is trivial at least
    compared to the other things here …

  5. > Aaron: what do you mean “only return bessel when necessary”?…

    Hmm. I though I saw you calling expand_func in one of your example pastes. Did you change them, or was I just seeing things?

    > The last combsimp is right, but actually it can be simplified further…

    Well, the order you enter a Mul has no baring on the final order in the .args, because they’re always canonically sorted. You’ll just have to improve your algorithm 🙂

    > What about symbolic/unevaluated transforms?

    This is a good idea. The transform methods should at the very least return unevaluated Integrals when it can’t compute rather than throwing an exception. But it would be useful to have these, because, for example, when you use Laplace transforms or Fourier transforms to solve a differential equation, you always have to keep track of the transform of the unknown function, which would be a little clunky to carry around as an unevaluated integral (especially when the transform has a coefficient). And like you said, you could always get the integral by doing .rewrite(Integral).

    You could also in theory use an unevaluated form to do some manual transformations on the transforms (like applying the convolution theorem for Fourier transforms). And you could pretty print them with script L, F, etc.

    But like you say, this should be a relatively trivial step.

  6. Yes, the algorithm is just wrong. It tries to combine any gammas there are, and thus eagerly combines some of those with arguments differing by 1/2. I’ll look at this again this week.

    The only expand functions in my examples I can see are regarding fourier transforms. Here I did simplify(expand_trig(expand_complex(result)))). This is because I want to get the result as a sinc-squared function (it comes out as a combination of complex exponentials). I think this hardly the general case … if it turns out it is, I might add it to the fourier_transform code (e.g. the inverse laplace transform code does expand_complex etc).

    @raoul: thanks. It is much more work than I anticipated but I feel like I’m finally getting to a sensible state.

  7. Oh, I see now that I was talking about the first pastebin. For some reason, I thought the last example was the same as the second to last, except the Bessel function was being converted to an erf. But I see now that this is wrong. So just ignore all that about expand_func().

Leave a Reply

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

You are commenting using your 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: