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()

Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

%d bloggers like this: