Menu
▾
▴
Re: [Matplotlib-users] Lorenz - solution
|
From: Prahas D. N. <pra...@gm...> - 2015-03-10 15:54:01
|
I think you need to ask Jake Vanderplas -- the code is all his! See the link in the email to get to his blog... Thanks again! On Tue, Mar 10, 2015 at 8:49 AM, Benjamin Root <ben...@ou...> wrote: > +1000!! > > Great job! Would you mind if I clean it up a bit and add it to the > mplot3d/animation gallery? Full credit, of course. > > > On Tue, Mar 10, 2015 at 11:30 AM, Prahas David Nafissian > <pra...@gm...> wrote: >> >> 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() >> >> ******************************** >> >> >> ------------------------------------------------------------------------------ >> Dive into the World of Parallel Programming The Go Parallel Website, >> sponsored >> by Intel and developed in partnership with Slashdot Media, is your hub for >> all >> things parallel software development, from weekly thought leadership blogs >> to >> news, videos, case studies, tutorials and more. Take a look and join the >> conversation now. http://goparallel.sourceforge.net/ >> _______________________________________________ >> Matplotlib-users mailing list >> Mat...@li... >> https://lists.sourceforge.net/lists/listinfo/matplotlib-users > > |
