|
From: Prahas D. N. <pra...@gm...> - 2015-03-10 15:30:29
|
Friends, I thought you'd like to see the solution. Many thanks to Jake Vanderplas for his code and teachings: https://jakevdp.github.io/blog/2013/02/16/animating-the-lorentz-system-in-3d/ If you start a new IP Notebook session, run as your first entry: %pylab and then copy and paste the text below and run it, you should be good to go (on a Mac, at least). There are several parameters I've changed from his original, and I've commented as I've changed. The original code is at the link above. There is one error in his code -- I've documented it below. Again, thanks to the community, Jake, and Ben Root. --Prahas ****************** import numpy as np from scipy import integrate from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import Axes3D from matplotlib.colors import cnames from matplotlib import animation # orig value of N_traj was 20 -- very cool this way. N_trajectories = 1 def lorentz_deriv((x, y, z), t0, sigma=10., beta=8./3, rho=28.0): """Compute the time-derivative of a Lorentz system.""" return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z] # Choose random starting points, uniformly distributed from -15 to 15 np.random.seed(1) # changing from -15,30 to 10,5 below starts the drawing in the middle, # rather than getting the long line from below.... # if using N_Traj > 1, return to orig values. # x0 = -15 + 30 * np.random.random((N_trajectories, 3)) x0 = 10 + 5 * np.random.random((N_trajectories, 3)) # Solve for the trajectories # orig values: 0,4,1000 # 3rd value -- lower it, it gets choppier. # 2nd value -- increase it -- more points, but speedier. # change middle num from 4 to 15 -- this adds points!!!!!!!! t = np.linspace(0, 40, 3000) x_t = np.asarray([integrate.odeint(lorentz_deriv, x0i, t) for x0i in x0]) # Set up figure & 3D axis for animation fig = plt.figure() ax = fig.add_axes([0, 0, 1, 1], projection='3d') # changing off to on below adds axises. slows it down but you # can fix that with interval value in the animation call ax.axis('on') # choose a different color for each trajectory colors = plt.cm.jet(np.linspace(0, 1, N_trajectories)) # set up lines and points -- this is a correction from # the orig jake code. the next four lines... lines = [ax.plot([], [], [], '-', c=c)[0] for c in colors] pts = [ax.plot([], [], [], 'o', c=c)[0] for c in colors] # prepare the axes limits ax.set_xlim((-25, 25)) ax.set_ylim((-35, 35)) ax.set_zlim((5, 55)) # set point-of-view: specified by (altitude degrees, azimuth degrees) ax.view_init(30, 0) # initialization function: plot the background of each frame def init(): for line, pt in zip(lines, pts): line.set_data([], []) line.set_3d_properties([]) pt.set_data([], []) pt.set_3d_properties([]) return lines + pts # animation function. This will be called sequentially with the frame number def animate(i): # we'll step two time-steps per frame. This leads to nice results. i = (2 * i) % x_t.shape[1] for line, pt, xi in zip(lines, pts, x_t): x, y, z = xi[:i].T line.set_data(x, y) line.set_3d_properties(z) pt.set_data(x[-1:], y[-1:]) pt.set_3d_properties(z[-1:]) # changed 0.3 to 0.05 below -- this slows the rotation of the view. # changed 30 to 20 below # changing 20 to (20 + (.1 * i)) rotates on the Z axis. trippy. ax.view_init(10, 0.1 * i) # ax.view_init(10, 100) fig.canvas.draw() return lines + pts # instantiate the animator. I've deleted the blit switch (for Mac) # enlarging frames=500 works now -- it failed before because I didn't give it # enough data -- by changing the t=np.linspace line above I generate more points. # interval larger slows it down # changed inteval from 30 to 200, frames from 500 to 3000 anim = animation.FuncAnimation(fig, animate, init_func=init, frames=3000, interval=200) # Save as mp4. This requires mplayer or ffmpeg to be installed. COMPLEX!!!!! # Instead, use a screen record program: Quicktime on the Mac; MS Expression Encoder on PC. # anim.save('PDNlorentz_attractor.mp4', fps=15, extra_args=['-vcodec', 'libx264']) plt.show() ******************************** |