In [1]:

    

%pylab inline --no-import-all
import matplotlib
import matplotlib.pyplot as plt
from pylab import *
import pickle
import sympy


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

eigenvalues = pickle.load(open('eigenvals.pickle', 'rb'))


    

In [3]:

    

eigenvalues[1][1].as_real_imag()
is_complex = lambda x: x.as_real_imag()[1] > 1e-20
filtered_eigenvalues = [(row[0], max(filter(is_complex, row[1:]), key=lambda x:x.as_real_imag()[0]).as_real_imag()[0]) for row in eigenvalues[4:]]
a15s = [row[0] for row in filtered_eigenvalues]
lambdas = [row[1] for row in filtered_eigenvalues]


    

In [5]:

    

figure()
plot(a15s, lambdas)
xlim(0, 2)
xlabel("Parameter—delta_p = a15")
ylim(-0.5, 0.1)
ylabel("Re(lambda)")
axvline(0.5707, color='black')
axhline(color='black')
show()


    

In [45]:

    

matplotlib.rcParams.update({'font.size': 16})
matplotlib.rcParams.update({'figure.autolayout': True})
figure()
plot(a15s, lambdas)
xlabel("Parameter—delta_p = a15")
ylabel("Real Part of Eigenvalue")
axvline(0.5707, color='black')
axhline(color='black')
ylim(-0.5, 0.1)
xlim(-5, 20)
annotate("Hopf bifurcation", (0.5707, 0.05), xytext=(5.2,0.04), arrowprops=dict(facecolor='black',lw=0.01))
#show()
savefig("eigenvalue_plot.png")


    

In [23]:

    

import scipy.integrate
scipy.integrate.simps(lambdas[:51], a15s[:51])


    

    
Out[23]:





-0.0365763523726227

In [94]:

    

a15s[50]


    

    
Out[94]:





0.570751

In [95]:

    

a15s[139]


    

    
Out[95]:





1.45449

In [88]:

    

upper = 140
scipy.integrate.simps(lambdas[51:upper], a15s[51:upper])


    

    
Out[88]:





0.00997222710664623

In [87]:

    

a15s[140]


    

    
Out[87]:





1.46442

In [69]:

    

import numpy


    

In [83]:

    

fitted = numpy.polyfit(a15s, lambdas, 20)
p = numpy.poly1d(fitted)


    

    




/usr/lib/python3.3/site-packages/numpy/lib/polynomial.py:587: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

In [84]:

    

figure()
plot(a15s, p(a15s))
axvline(0.5707, color='black')
axhline(color='black')
show()


    

In [ ]: