[Numpy-discussion] Recognizing a cycle in a vector
Oscar Benjamin
oscar.j.benjamin at gmail.com
Thu Dec 3 06:23:26 EST 2015
On 3 December 2015 at 08:45, Manolo Martínez <manolo at austrohungaro.com> wrote:
>> > >> Is there any way to check for cycles in this situation?
>> >
>> > > Fast fourier transform (fft)?
>> >
>> > +1 For using a discrete Fourier transform, as implemented by numpy.fft.fft.
>> > You mentioned that you sample at points which do not correspond with the
>> > period of the signal; this introduces a slight complexity in how the
>> > Fourier transform reflects information about the original signal. I attach
>> > two documents to this email with details about those (and other)
>> > complexities. There is also much information on this topic online and in
>> > signal processing books.
>>
>
> So, I thought I'd report back on what I've ended up doing. Given that
> the cycles I find in my data are usually very close to sine waves, the
> following works well enough:
>
>
> def periodic_vector(vector):
> """
> Take the FFT of a vector, and eliminate all components but the
> two main ones (i.e., the static and biggest sine amplitude) and
> compare the reconstructed wave with the original. Return true if
> close enough
> """
> rfft = np.fft.rfft(vector)
> magnitudes = np.abs(np.real(rfft))
> choice = magnitudes > sorted(magnitudes)[-3]
> newrfft = np.choose(choice, (np.zeros_like(rfft), rfft))
> newvector = np.fft.irfft(newrfft)
> return np.allclose(vector, newvector, atol=1e-2)
>
>
> This is doing the job for me at the moment, but there are, that I can
> see, a couple of things that could be improved (and surely more that I
> cannot see):
>
> 1) this func sorts the absolute value of the amplitudes to find the two
> most important components, and this seems overkill for large vectors.
>
> 2) I'm running the inverse FFT, and comparing to the initial vector, but
> it should be possible to make a decision solely based on the size of
> terms in the FFT itself. I'm just not confident enough to design a test
> based on that.
>
> Anyway, thanks to those who pointed me in the right direction.
If what you have works out fine for you then feel free to ignore this but...
The more common way to find periodic orbits in ODEs is to pose the
question as a boundary-value problem (BVP) rather than seek orbital
patterns in the solution of an initial value problem (IVP). BVP
methods are more robust, can find unstable orbits, detect period
doubling etc. I would use the DFT method in the case that the ODEs are
of very high dimension and/or if the orbits in question are perhaps
quasi-periodic or if I can only really access the ODEs through the
output of an IVP solver. In the common case though the BVP method is
usually better. Something is written about finding periodic orbits
here:
http://www.scholarpedia.org/article/Periodic_orbit#Numerical_Methods_for_Finding_Periodic_Orbits
There are two basic ways to do this as a BVP problem: the shooting
method and the mesh. For your purposes the shooting method may suffice
so I'll briefly describe it.
Your ODEs must be of dimension at least 2 or you wouldn't have
periodic solutions. Consider x[n] to be the state vector at the nth
timestep. Suppose that x~ is part of an exact periodic orbit of the
ODE x' = f(x) with f(x) some vector field. Define P as the plane
containing x~ and normal to f(x~). The periodic orbit (if it exists)
must curve around and end up back at x~ pointing in the same
direction. For sinusoidal-ish orbits it will cross P twice, once at x~
and once somewhere else heading in something like the opposite
direction. If the orbit is more wiggly it may cross P more time but
always an even number of times before reaching x~.
Now suppose you have some guess x[0] which is close to a periodic
orbit. The true orbit should cross the plane P' generated by x[0],
f(x[0]) somewhere near x[0] pointing in approximately the same
direction. So use an ODE solver to iterate forward making x[1], x[2]
etc. until you cross the plane once and then twice coming back
somewhere near x[0]. Now you have x[n] and x[n+1] close-ish to x[0]
which lie on either side of the plane crossing in the same direction
as f(x[0]). You can now use the bisect method to take smaller and
larger timesteps from x[n] until your trajectory hits the plane
exactly. at this point your orbit is at some point x* which is on P'
near to x[0]. We now have an estimate of the period T of the orbit.
What I described in the last paragraph may be sufficient for your
purposes: if x* is sufficiently close to x[0] then you've found the
orbit and if not then it's not periodic. Usually though there is
another step: Define a function g(x, T) which takes a point x on the
plane P' and iterates the ODEs through a time T. You can put this into
a root-finder to solve g(x, T) = x for T and x. Since x is
N-dimensional we have N equations. However we have constrained x to
lie on a plane so we have N-1 degrees of freedom in choosing x but we
also want to solve for T which means we have N equations for N
unknowns.
Putting g into a root-finder as I described is called the shooting
method for BVPs. A more robust method uses a mesh and something like
the central difference method to solve for a set of points on the
orbit but this may not be necessary in your case.
Libraries for doing this (using more advanced methods than I have just
described) already exist so you may want to simply use them rather
than reinvent this particular wheel.
--
Oscar
More information about the NumPy-Discussion
mailing list