I have a function which I have obtained from numerical integration of a differential equation, and I would like to take its Fourier transform. What are good practices for doing this?
To make things more precise, suppose I have a function like the following:
solution = u /. First@NDSolve[{
u'[t] == -E^(-(t^2/2)) t Cos[5 t] - 5 E^(-(t^2/2)) Sin[5 t] t,
u[0] == 1
}, u, {t, 0, 50}]
(Note that, while this example has an analytical solution, my full problem does not.) I am using NDSolve
since it appears to be the best way to obtain functions $F$ of the form $F(t)=\int_{t_0}^t f(\tau) \,\text d\tau$ where $f$ is known. I am interested in the Fourier transform of this solution.
I know of course that I can sample this function over a discrete set of points and then use the discrete Fourier transform Fourier
to obtain a fair approximation:
ListPlot[
Abs[
Fourier[solution /@ Range[0, 50, 0.01]]
]
, PlotRange -> {{0, 150}, {0, 1}}]
However, this completely ignores the fact that the InterpolatingFunction
itself is already a discrete representation of the desired function. As it stands, the code above includes interpolation from a discrete data set into a new discrete data set, and this can only
- lose information, in the case that the sample grid is too sparse, or alternatively
- give the impression that it contains more information than was already in the original data, if the sample grid is too closely spaced.
Either way, this extra sampling step has no communication with the internals of the InterpolatingFunction
and as such it cannot know how much information it contains or how good it is. I am looking for a method which uses this information as optimally as possible. The accuracy of the ODE solution, and therefore of its transform, are best controlled from the options of NDSolve
, and the subsequent transformations of the data should take heed of this.
I am aware of NFourierTransform
, which, as I understand, uses NIntegrate
internally; I'm also aware that NIntegrate
natively supports InterpolatingFunction
objects, which is presumably still the case if they are multiplied by a further factor of $e^{i\omega t}$. (On the other hand, the behaviour of this factor can rather radically change the overall behaviour of the whole integrand from 'slow' into highly oscillatory, so that some attention to Method
would be warranted.)
It could be, though, that the Fast nature of Fourier
might be enough to outstrip such advantages, and still be faster than NFourierTransform
even when over-blasted on the accuracy. Speed is definitely an issue, as I will be performing automated scans over many such spectra, so I want them to be accurate enough but not waste time overcalculating.
So: what are good ways of performing this transform?
Fourier
. $\endgroup$ – Daniel Lichtblau May 13 '14 at 14:20Compile
and I should never abandon discrete representations of my data. If I had this need and performance in mind, I would roll a discrete solver for this differential equation. $\endgroup$ – LLlAMnYP Jun 29 '17 at 12:08