![](https://secure.gravatar.com/avatar/8ff70510c293426a72aca34f2a00b21d.jpg?s=120&d=mm&r=g)
Dear all, Suppose that I have a vector with the numerical solution of a differential equation -- more concretely, I am working with evolutionary game theory models, and the solutions are frequencies of types in a population that follows the replicator dynamics; but this is probably irrelevant. Sometimes these solutions are cyclical, yet I sample at points which do not correspond with the period of the cycle, so that np.allclose() cannot be directly applied. Is there any way to check for cycles in this situation? Thanks for any advice, Manolo
![](https://secure.gravatar.com/avatar/e93e18f71a5821a54c233690506bdbf7.jpg?s=120&d=mm&r=g)
Manolo,
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.
![](https://secure.gravatar.com/avatar/8ff70510c293426a72aca34f2a00b21d.jpg?s=120&d=mm&r=g)
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. Manolo
![](https://secure.gravatar.com/avatar/664d320baa05c827ff08ed361fe77769.jpg?s=120&d=mm&r=g)
On 3 December 2015 at 08:45, Manolo Martínez <manolo@austrohungaro.com> wrote:
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_Fin... 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
![](https://secure.gravatar.com/avatar/8ff70510c293426a72aca34f2a00b21d.jpg?s=120&d=mm&r=g)
Dear Oscar,
If what you have works out fine for you then feel free to ignore this but... [snip]
Talk about things I cannot see :) Thanks a lot for that very detailed explanation. I will certainly look into settting my problem up as a BVP. Incidentally, is there any modern textbook on numerical solving of ODEs that you could recommend? Thanks again, Manolo
![](https://secure.gravatar.com/avatar/664d320baa05c827ff08ed361fe77769.jpg?s=120&d=mm&r=g)
On 3 December 2015 at 11:58, Manolo Martínez <manolo@austrohungaro.com> wrote:
Not particularly. The shooting and mesh methods are described in chapter 17 of the Numerical Recipes in C book which is relatively accessible: http://apps.nrbook.com/c/index.html In terms of out of the box software I can recommend auto and xpp. Each is esoteric and comes with a clunky interface. XPP has a strange GUI and auto is controlled through Python bindings using IPython as frontend. -- Oscar
![](https://secure.gravatar.com/avatar/7e9e53dbe9781722d56e308c32387078.jpg?s=120&d=mm&r=g)
On 2015/12/02 10:45 PM, Manolo Martínez wrote:
1) this func sorts the absolute value of the amplitudes to find the two most important components, and this seems overkill for large vectors.
Try inds = np.argpartition(-np.abs(ft), 2)[:2] Now inds holds the indices of the two largest components. Eric
![](https://secure.gravatar.com/avatar/e93e18f71a5821a54c233690506bdbf7.jpg?s=120&d=mm&r=g)
Manolo,
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.
![](https://secure.gravatar.com/avatar/8ff70510c293426a72aca34f2a00b21d.jpg?s=120&d=mm&r=g)
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. Manolo
![](https://secure.gravatar.com/avatar/664d320baa05c827ff08ed361fe77769.jpg?s=120&d=mm&r=g)
On 3 December 2015 at 08:45, Manolo Martínez <manolo@austrohungaro.com> wrote:
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_Fin... 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
![](https://secure.gravatar.com/avatar/8ff70510c293426a72aca34f2a00b21d.jpg?s=120&d=mm&r=g)
Dear Oscar,
If what you have works out fine for you then feel free to ignore this but... [snip]
Talk about things I cannot see :) Thanks a lot for that very detailed explanation. I will certainly look into settting my problem up as a BVP. Incidentally, is there any modern textbook on numerical solving of ODEs that you could recommend? Thanks again, Manolo
![](https://secure.gravatar.com/avatar/664d320baa05c827ff08ed361fe77769.jpg?s=120&d=mm&r=g)
On 3 December 2015 at 11:58, Manolo Martínez <manolo@austrohungaro.com> wrote:
Not particularly. The shooting and mesh methods are described in chapter 17 of the Numerical Recipes in C book which is relatively accessible: http://apps.nrbook.com/c/index.html In terms of out of the box software I can recommend auto and xpp. Each is esoteric and comes with a clunky interface. XPP has a strange GUI and auto is controlled through Python bindings using IPython as frontend. -- Oscar
![](https://secure.gravatar.com/avatar/7e9e53dbe9781722d56e308c32387078.jpg?s=120&d=mm&r=g)
On 2015/12/02 10:45 PM, Manolo Martínez wrote:
1) this func sorts the absolute value of the amplitudes to find the two most important components, and this seems overkill for large vectors.
Try inds = np.argpartition(-np.abs(ft), 2)[:2] Now inds holds the indices of the two largest components. Eric
participants (5)
-
Daniel Sank
-
Elliot Hallmark
-
Eric Firing
-
Manolo Martínez
-
Oscar Benjamin