# Calculating fixed point inside limit cycle

So I'm working with a rather complicated dynamical system. Instead of writing it all out. It's probably easier if you just clone my git repository.

git clone https://gitlab.com/mdornfe1/vortex_acoustic


The dynamics of the system is encoded in the function flow inside vortex_acoustic.py. If I integrate the system from rest I find oscillations for certain values of the parameters. The system is 13 dimensional, but here is a plot of the first variable vs. time

So it seems like the system is approaching a limit cycle attractor. From my understanding of dynamical systems there should be an unstable fixed point inside that limit cycle. I want to try to calculate it so I can analyze its stability. My idea is that the fixed point should be close to the mean of the oscillating quantities, so I use that as a starting value for a numerical root solver. Here's the code for finding the fixed point, but you can also just run stack_exchange.py.

import pickle, numpy as np, matplotlib.pyplot as plt
import vortex_acoustic as va
from scipy import optimize
from constants import *

T = 5000
dt = 0.05
tn_transients = int(0.9 * T / dt)

def calc_fixed_point(seed, params):
return optimize.fsolve(va.flow_star, seed, params)

def calc_closest_fixed_point(y, params):
seed = np.mean(y[tn_transients:, :], 0)
seed[Nq+1:2*Nq+1] = 0
y_star = calc_fixed_point(seed, params)

return y_star, seed

with open('params.p', 'r') as f:

plt.plot(y[:,0])

y_star, seed = calc_closest_fixed_point(y, params)


The function calc_closest_fixed_point takes the output of the simulation and the params passed to the simulation. It calculates the mean of the fluctuating quantities after the transients has died out. It sets the derivatives to 0 and uses that quantity seed as a starting point for fsolve. The function calc_fixed_point then finds the fixed point of the system y_star. I thought this approach would show that seed and y_star are very close to each other. But for this example I'm seeing they're quite different. Furthermore the actual value isn't inside the limit cycle.

print(seed)
>>>array([ 0.18307736, -0.11207637,  0.00382286,  0.        ,  0.        ,
-0.02016129,  0.08881281, -0.07852504,  0.07439292,  0.0109214 ,
-0.00894938, -0.02896082,  0.00833759])
print(y_star)
>>>array([  1.06466138e+00,  -6.65605305e-01,   2.15920825e-02,
-5.18569988e-25,   1.44087666e-24,   7.55288134e-03,
5.17875853e-01,  -4.41152924e-01,   4.72080951e-01,
4.86813549e-02,  -6.04678149e-02,  -1.42533599e-01,
5.90028044e-02])


Also seed is pretty obviously not a fixed point. If you calculate the function flow at seed the derivatives are pretty far from 0.

va.flow(seed,0,*params)
>>>array([-0.00543102,  0.        ,  0.        , -0.00113827, -0.00042674,
-0.00973662, -0.00775962,  0.01449168,  0.01369892,  0.00331889,
0.00502761, -0.00026469, -0.00076228])


Anyone know what's going? I'm correct in assuming there is a fixed point inside the limit cycle right?

• Can you check how close seed is to the fixed point by evaluating y_star(seed), as well as y_star's Jacobian? It could also be worth trying the other scipy's root finding methods too to see if they do better. – Kirill Oct 18 '16 at 0:28
• @Kirill See my edit. seed is definitely not a fixed point, but the theory makes me thing it should be. So that's where I'm kind of confused. I'm not sure how to analyze the stability of the limit cycle without a fixed point inside of it. – mdornfe1 Oct 18 '16 at 0:34
• No, you don't have to have an unstable fixed point to have a limit cycle unless you are in 2D. Also, if you want to find unstable fixed points, you can quickly do so by running the dynamics backwards in time (stable FP becomes unstable and vice versa). – Memming Oct 19 '16 at 8:10
• @Memming, So how can you analyze the stability of a limit cycle with no fixed point inside? I don't think I've read about this before. – mdornfe1 Oct 19 '16 at 15:25
• In 3D, an example would be a limit cycle on a torus. No need for unstable fixed points... – Memming Oct 19 '16 at 16:41