7

I wrote this function NFourierTransform, which takes a function $f(k)$ and numerically calculates the fourier transform integral for discrete values of $k \in [k_{\text{min}},k_{\text{max}}]$, finally returning an InterpolatingFunction object.

NFourierTransform[f_Function, {kmin_, kmax_}] := 
 Interpolation@
  Table[{k, Chop@NIntegrate[f@x E^(-I k x), {x, -Infinity, Infinity}]},
   {k, kmin, kmax, (kmax - kmin)/100}]

In my application (calculating the time propagation of wave functions) I need to evaluate NFourierTransform for a function $\psi(k,t)$, where $t$ assumes discrete values in some interval $[t_{\text{min}},t_{\text{max}}]$. So effectively I create a table of NFourierTransform.

TimePropagate[f_Function, kl : {kmin_, kmax_}, {tl__}] :=
 Quiet@Table[
   NFourierTransform[f[#] Exp[-(#^2/2) t] &, kl], {t, tl}]

Calculating a very simple example with only 2 time values, e.g. TimePropagate[Exp[-Abs@#] &, {-3, 3}, {0, 0.1, 0.1}] takes about 20 seconds to evaluate.

Is there any way to use Compile to speed up the calculations? As far as I know that's only possible for numeric function arguments. What are, in your experience, suitable setings for NIntegrate options such as MaxRecursion or AccuracyGoal, and how do they effect evaluation time?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
einbandi
  • 4,024
  • 1
  • 23
  • 39
  • I never integrate over t. Each value for t is inserted, and the integration is over some variable x, where k ranges from kmin to kmax. – einbandi Nov 23 '12 at 14:28
  • I looked into that. My problem is that I don't really understand the mathematics behind discrete fourier transform (Fourier). I'm much more comfortable with continuous FT. – einbandi Nov 23 '12 at 14:34
  • It would be helpful to know what can we assume of your function. If you can be sure about some bandwidth and support then using Fourier is surely a fast way. Meaning, a maximum frequency in which the function has a not too small fourier transform, and a maximum time in which the function is not too small. – Rojo Nov 23 '12 at 14:43
  • If that's the case, I remembered a previous question in which my poor answer got the least votes and was barely seen. I don't even remember if it works but you could try it and let me know :P – Rojo Nov 23 '12 at 14:44
  • That estimate function is a function of the frequency in Hz. To make that fit your fourier transform conventions you should evaluate it in k/(2 Pi) I think. Also, as is it assumes that your function is near 0 for negative arguments – Rojo Nov 23 '12 at 14:46
  • Thanks for the tip, I'll look into your function. Might take me some time to fully understand what it does, though. – einbandi Nov 23 '12 at 14:59
  • If your function f is slow to evaluate and you are evaluating it in the same places for each t in tl, you can implement memoization which means to cache the previous results. Look for that word in this site and you'll find lots of examples – Rojo Nov 23 '12 at 15:14
  • One thing that comes to mind: you might want to explicitly set the Method option of NIntegrate[] to something that's equipped to handle oscillatory integrals, like "DoubleExponentialOscillatory" or "ExtrapolatingOscillatory". – J. M.'s missing motivation Nov 23 '12 at 15:46
  • 1
    What Rojo said. The responses here seem to cover this topic. My answer got but one vote more than Rojo's. I think that fish might do what you want though (at least once suitably scaled). – Daniel Lichtblau Nov 26 '12 at 07:20
  • Thanks. I still have some troubles understanding the mathematics behind deiscrete FT, but I'll definitely look into that. – einbandi Nov 26 '12 at 14:32

1 Answers1

10

I had a play with various Compile options and didn't get anywhere (I managed to make it slower though!). However, you can get a nice little speed boost using ParallelTable. Your original on my machine:

NFourierTransform[f_Function, {kmin_, kmax_}] := 
 Interpolation@
  Table[{k, 
    Chop@NIntegrate[f@x E^(-I k x), {x, -Infinity, Infinity}]}, {k, 
    kmin, kmax, (kmax - kmin)/100}]
TimePropagate[f_Function, kl : {kmin_, kmax_}, {tl__}] := 
 Quiet@Table[NFourierTransform[f[#] Exp[-(#^2/2) t] &, kl], {t, tl}]

TimePropagate[Exp[-Abs@#] &, {-3, 3}, {0, 0.1, 0.1}] // AbsoluteTiming runs in 5.65 seconds. I launch some kernels

LaunchKernels[]

and throw in a ParallelTable

NFourierTransform[f_Function, {kmin_, kmax_}] := 
 Interpolation@
  ParallelTable[{k, 
    Chop@NIntegrate[f@x E^(-I k x), {x, -Infinity, Infinity}]}, {k, 
    kmin, kmax, (kmax - kmin)/100}]
TimePropagate[f_Function, kl : {kmin_, kmax_}, {tl__}] := 
 Quiet@Table[NFourierTransform[f[#] Exp[-(#^2/2) t] &, kl], {t, tl}]

to get an execution time of 2.04 seconds for the same function call: TimePropagate[Exp[-Abs@#] &, {-3, 3}, {0, 0.1, 0.1}] // AbsoluteTiming

A speed-up of almost 3 times.

rm -rf
  • 88,781
  • 21
  • 293
  • 472
WalkingRandomly
  • 4,107
  • 1
  • 20
  • 36
  • Hi Mike, welcome to Mathematica at Stack Exchange! I've been reading your blog with great interest for a while now. It's great that you joined, and I hope you will like the community here. Gives me a pleasure to give you the first upvote on this site,+1. – Leonid Shifrin Nov 26 '12 at 12:51
  • Thanks. I have been keeping an eye on this site for a while (I promoted the original proposal - http://www.walkingrandomly.com/?p=4014 - hope it helped) and have learned many useful tricks here. I rarely answer anything because the site regulars here are too efficient for their own good :) – WalkingRandomly Nov 26 '12 at 13:43
  • Finally a useful tip that uses my original functions! Thanks a lot! – einbandi Nov 26 '12 at 14:30
  • Welcome to the site, Mike. I've edited your post to format it a little better. A tip: while backticks are useful for inline code, you can typeset blocks of code by indenting it by 4 spaces (or selecting it and pressing the {} button in the edit window). This way, it's also pretty printed with a custom Mathematica syntax highlighter (developed here!) – rm -rf Nov 26 '12 at 15:25
  • Ok, thanks for that. I did a very quick scan around the site to see how to do code formatting but didn't find anything. Probably my bad as it was a VERY quck scan :) – WalkingRandomly Nov 26 '12 at 16:03
  • Yeah, it can be a bit hard to find at times (they're linked to from the edit box, but I don't remember where), but for reference, here's the full list: http://mathematica.stackexchange.com/editing-help :) – rm -rf Nov 26 '12 at 18:58