# Lorenz - solution

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....

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

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

+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 <prahas.music@…287…> 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)

# x0 = -15 + 30 * np.random.random((N_trajectories, 3))

x0 = 10 + 5 * np.random.random((N_trajectories, 3))

# 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’)

ax.axis(‘on’)

# choose a different color for each trajectory

colors = plt.cm.jet(np.linspace(0, 1, N_trajectories))

# 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
``````

more points.

# changed inteval from 30 to 200, frames from 500 to 3000

anim = animation.FuncAnimation(fig, animate, init_func=init,

``````                           frames=3000, interval=200)
``````

# 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

Matplotlib-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/matplotlib-users

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.root@...1304...> 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 > <prahas.music@...287...> 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....

# 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,
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
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

That’s pretty swag!

···

On Tue, Mar 10, 2015 at 11:49 AM, Benjamin Root <ben.root@…1304…> 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.

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

Matplotlib-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/matplotlib-users

On Tue, Mar 10, 2015 at 11:30 AM, Prahas David Nafissian <prahas.music@…287…> 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)

# x0 = -15 + 30 * np.random.random((N_trajectories, 3))

x0 = 10 + 5 * np.random.random((N_trajectories, 3))

# 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’)

ax.axis(‘on’)

# choose a different color for each trajectory

colors = plt.cm.jet(np.linspace(0, 1, N_trajectories))

# 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
``````

more points.

# changed inteval from 30 to 200, frames from 500 to 3000

anim = animation.FuncAnimation(fig, animate, init_func=init,

``````                           frames=3000, interval=200)
``````

# 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

Matplotlib-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/matplotlib-users

@ Everyone: given the Lorenz code shared, is there a way
to generate a log file of the x,y,z points generated?

--Prahas

···

On Tue, Mar 10, 2015 at 8:59 AM, Adam Hughes <hughesadam87@...287...> wrote:

That's pretty swag!

On Tue, Mar 10, 2015 at 11:49 AM, Benjamin Root <ben.root@...1304...> 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 >> <prahas.music@...287...> 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....

# 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,
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
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

------------------------------------------------------------------------------
Dive into the World of Parallel Programming The Go Parallel Website,