# Rcjp's Weblog

## July 28, 2006

### Bifurcation

Filed under: python — rcjp @ 1:14 pm

I really should have made a better job of this with numpy instead of translating C++ version… todo. (Axis are just showing off the tex features…)

 #from pyx import * from pylab import * from matplotlib import rc from Numeric import * StabalisingIterations = 1000 RecordingIterations = 30 def rmmay(r, x):     return r*x*(1-x) def stabalise_logistic_mapping(mapping, r, x=0.1):     """let the mapping settle down before recording the values by     apply mapping function nsettle times starting at 0.1, returning the last nrecord of these"""     for n in xrange(StabalisingIterations):         x = mapping(r, x)     res = []     for n in xrange(RecordingIterations):         x = mapping(r, x)         res.append(x)     return res def lyapunov_exponent(mapping, r):     data = stabalise_logistic_mapping(mapping, r)     data = [log(abs(r*(1-2*x))) for x in data]     return sum(data)/RecordingIterations def lyapunov_data (mapping, rmin, rmax, rstep):     f = open('/tmp/ldata', 'w')     xdata = []     ydata = []     for r in arange(rmin, rmax, rstep):         ly = lyapunov_exponent(mapping,r)         xdata.append(r)         ydata.append(ly)         f.write('%f %f\n' %(r, ly))     f.close()     # Plot graph #     g = graph.graphxy(width=8) #     g.plot(graph.data.file('c:/tmp/ldata', x=1, y=2)) #     g.writePDFfile("c:/tmp/lyapunov") #     rc('text', usetex=True)     figure(1)     ax = axes([0.15, 0.1, 0.8, 0.7])     plot(xdata, ydata, linewidth = 1.0)     xlabel(r'\bf{time (s)}')     ylabel(r'\it{voltage (mV)}',fontsize=16)     title(r"\TeX\ is Number $\displaystyle\sum_{n=1}^\infty\frac{-e^{i\pi}}{2^n}$!",           fontsize=16, color='r')     savefig('/tmp/lyapunov.eps') #    grid(True)     show() if __name__ == '__main__':     lyapunov_data(rmmay, 2.95, 4.0, 0.005) 

actually… found some *much* nicer code from Cornell Uni’s interesting course

 from scipy import * from pylab import * f = lambda x,r : r * x * ( 1 -x ) rlist = linspace( 2, 4, 800) X = [ 0.5 * ones_like(rlist), ] for i in arange(0,10000):   X += [ f(X[ -1], rlist) , ] X = hsplit( vstack(X[-2000:]), rlist.size) from scipy import stats H = map( lambda Z : stats.histogram( Z, defaultlimits=(0,1), numbins=300 )[0],X) H = map( lambda Z : 1-Z/Z.max(), H ) H = vstack(H) figure() imshow( rot90(H), aspect = 'auto' , extent = [2, 4, 0, 1]) bone() xlabel('r') ylabel(r'$X_{n \rightarrow \infty}$') show()