[Matplotlib-users] Lorenz - solution

[Prahas D. N. <pra...@gm...>]

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

********************************