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_acousticThe 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. timeSo 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 pltimport vortex_acoustic as vafrom scipy import optimizefrom constants import *T = 5000dt = 0.05tn_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, seedy = np.load('simulation.npy')with open('params.p', 'r') as f:    params = pickle.load(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?
Calculating fixed point inside limit cycle
roots;fixed point