#!/usr/bin/python -Wall # ================================================================ # John Kerl # kerl at math dot arizona dot edu # 2009-10-08 # Please see http://math.arizona.edu/~kerl/doc/python-talk.pdf for # more information. # ================================================================ import sys, pylab # ---------------------------------------------------------------- def project_3D_to_2D(x, y, z, Pxhat, Pyhat, Pzhat): u = Pxhat[0] * x + Pyhat[0] * y + Pzhat[0] * z v = Pxhat[1] * x + Pyhat[1] * y + Pzhat[1] * z return [u, v] # ---------------------------------------------------------------- # Simulation parameters tmax = 10.0 h = 0.002 N = int(tmax/h) # Or: # h = 0.002 # N = 5000 # tmax = h * N # ---------------------------------------------------------------- # Initial conditions x0 = 10.0 y0 = 0.0 z0 = 10.0 # ---------------------------------------------------------------- # Functional parameters s = 10.0 r = 28.0 b = 2.666667 # ---------------------------------------------------------------- # Output parameters if len(sys.argv) == 2: h = float(sys.argv[1]) N = int(tmax/h) # ---------------------------------------------------------------- # Integration # For projecting 3D (x,y,z) down to 2D (u,v) for plotting. Pxhat = [-0.25, -0.25]; Pyhat = [1, 0]; Pzhat = [0, 1] # Apply initiail conditions. x = x0; y = y0; z = z0 # Build up lists (starting from empty lists) of (u,v) pairs. us = [] vs = [] t = 0 while t <= tmax: # Cheesy Euler integrator. xp = s * (y - x) yp = r * x - y - x * z zp = x * y - b * z x += h * xp y += h * yp z += h * zp [u, v] = project_3D_to_2D(x, y, z, Pxhat, Pyhat, Pzhat) us.append(u) vs.append(v) t += h # ---------------------------------------------------------------- # Plot pylab.figure() pylab.plot(us,vs) pylab.show() # Or: #pylab.savefig('myfile.eps', dpi=150) # To disk