from scipy.integrate import odeint from math import * import numpy as np import pylab as pl s = 10.0 b = 8.0/3.0 x0 = 0.0 y0 = 0.0 z0 = 0.0 # Lorenz System def F(X, t, r): x, y, z = X return [-s * x + s * y, r * x - y - x * z, x * y - b * z] # Fixpunkte C+ und C- def Cpm(pm, r): xi = sqrt(b*(r-1.0)) return np.array((pm * xi, pm * xi, r -1.0)) t = np.linspace(0, 10, 20001) def calc_orbit(r): # unstable eigenvector vec = np.array([-((-1 + s - sqrt(1 - 2 * s + 4 * r * s + s**2))/(2 * s)), 1, 0]) delta = -1E-8 erg, info = odeint(F, delta*vec, t, args=(r,), full_output=True) x, y, z = erg.transpose() return x, y, z def plot_orbit(r): x, y, z =calc_orbit(r) pl.plot(x, y) Cplus = Cpm(+1, r) Cminus = Cpm(-1, r) C0 = (0, 0, 0) FPs = np.array([Cplus, Cminus, C0]).T pl.plot(FPs[0], FPs[1], 'ko') #ko black circle pl.xlabel(r'$x$') pl.ylabel(r'$y$') #ra und rb Grenzen des Intervallschachtlungsverfahrens ra = 1 rb = 20 r0 = 0.5*(rb + ra) #plot_orbit(13) #pl.show() #import sys #sys.exit(0) while rb-ra > 1E-10: x, y, z = calc_orbit(r0) last_point = np.array([x[-1], y[-1], z[-1]]) Cplus = Cpm(+1, r0) Cminus = Cpm(-1, r0) if (abs(last_point-Cplus)).sum() < (abs(last_point-Cminus)).sum(): rb = r0 else: ra = r0 r0 = 0.5 * (ra + rb) print r0 t = np.linspace(0, 4, 20001) plot_orbit(r0) pl.show() #plot_orbit(r0) #### 3D plot # #from mpl_toolkits.mplot3d import Axes3D #fig = pl.figurel() #ax = Axes3D(fig) #def plot_orbit3D(r): # x, y, z = calc_orbit(r) # pl.plot(x, y) # Cplus = Cpm(+1, r) # Cminus = Cpm(-1, r) # C0 = (0, 0, 0) # ax.plot(x, y, z) # ax.plot(FPs, 'ko') # ax.set_xlabel('$x$') # ax.set_ylabel('$y$') # ax.set_zlabel('$z$') #plot_orbit(r0) #pl.show() #pl.savefig('period-doubling_roessler.pdf')