Exploring the Lorenz System of Differential Equations

In this Notebook we explore the Lorenz system of differential equations:

\[\begin{split}\begin{aligned} \dot{x} & = \sigma(y-x) \\ \dot{y} & = \rho x - y - xz \\ \dot{z} & = -\beta z + xy \end{aligned}\end{split}\]

This is one of the classic systems in non-linear differential equations. It exhibits a range of different behaviors as the parameters ((\sigma\), \(\beta\), \(\rho\)) are varied.

Imports

First, we import the needed things from IPython, NumPy, Matplotlib and SciPy.

[1]:
%matplotlib inline
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-9e3324102725> in <module>
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')

~/checkouts/readthedocs.org/user_builds/ipywidgets/conda/latest/lib/python3.6/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line, _stack_depth)
   2312                 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
   2313             with self.builtin_trap:
-> 2314                 result = fn(*args, **kwargs)
   2315             return result
   2316

</home/docs/checkouts/readthedocs.org/user_builds/ipywidgets/conda/latest/lib/python3.6/site-packages/decorator.py:decorator-gen-108> in matplotlib(self, line)

~/checkouts/readthedocs.org/user_builds/ipywidgets/conda/latest/lib/python3.6/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
    185     # but it's overkill for just that one bit of state.
    186     def magic_deco(arg):
--> 187         call = lambda f, *a, **k: f(*a, **k)
    188
    189         if callable(arg):

~/checkouts/readthedocs.org/user_builds/ipywidgets/conda/latest/lib/python3.6/site-packages/IPython/core/magics/pylab.py in matplotlib(self, line)
     97             print("Available matplotlib backends: %s" % backends_list)
     98         else:
---> 99             gui, backend = self.shell.enable_matplotlib(args.gui.lower() if isinstance(args.gui, str) else args.gui)
    100             self._show_matplotlib_backend(args.gui, backend)
    101

~/checkouts/readthedocs.org/user_builds/ipywidgets/conda/latest/lib/python3.6/site-packages/IPython/core/interactiveshell.py in enable_matplotlib(self, gui)
   3400         """
   3401         from IPython.core import pylabtools as pt
-> 3402         gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)
   3403
   3404         if gui != 'inline':

~/checkouts/readthedocs.org/user_builds/ipywidgets/conda/latest/lib/python3.6/site-packages/IPython/core/pylabtools.py in find_gui_and_backend(gui, gui_select)
    274     """
    275
--> 276     import matplotlib
    277
    278     if gui and gui != 'auto':

ModuleNotFoundError: No module named 'matplotlib'
[2]:
from ipywidgets import interact, interactive
from IPython.display import clear_output, display, HTML
[3]:
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
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-3-a8584716cd6d> in <module>
      1 import numpy as np
----> 2 from scipy import integrate
      3
      4 from matplotlib import pyplot as plt
      5 from mpl_toolkits.mplot3d import Axes3D

ModuleNotFoundError: No module named 'scipy'

Computing the trajectories and plotting the result

We define a function that can integrate the differential equations numerically and then plot the solutions. This function has arguments that control the parameters of the differential equation ((\sigma\), \(\beta\), \(\rho\)), the numerical integration (N, max_time) and the visualization (angle).

[4]:
def solve_lorenz(N=10, angle=0.0, max_time=4.0, sigma=10.0, beta=8./3, rho=28.0):

    fig = plt.figure()
    ax = fig.add_axes([0, 0, 1, 1], projection='3d')
    ax.axis('off')

    # prepare the axes limits
    ax.set_xlim((-25, 25))
    ax.set_ylim((-35, 35))
    ax.set_zlim((5, 55))

    def lorenz_deriv(x_y_z, t0, sigma=sigma, beta=beta, rho=rho):
        """Compute the time-derivative of a Lorenz system."""
        x, y, z = x_y_z
        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, 3))

    # Solve for the trajectories
    t = np.linspace(0, max_time, int(250*max_time))
    x_t = np.asarray([integrate.odeint(lorenz_deriv, x0i, t)
                      for x0i in x0])

    # choose a different color for each trajectory
    colors = plt.cm.viridis(np.linspace(0, 1, N))

    for i in range(N):
        x, y, z = x_t[i,:,:].T
        lines = ax.plot(x, y, z, '-', c=colors[i])
        plt.setp(lines, linewidth=2)

    ax.view_init(30, angle)
    plt.show()

    return t, x_t

Let’s call the function once to view the solutions. For this set of parameters, we see the trajectories swirling around two points, called attractors.

[5]:
t, x_t = solve_lorenz(angle=0, N=10)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-dc018b50a68d> in <module>
----> 1 t, x_t = solve_lorenz(angle=0, N=10)

<ipython-input-4-864093adb516> in solve_lorenz(N, angle, max_time, sigma, beta, rho)
      1 def solve_lorenz(N=10, angle=0.0, max_time=4.0, sigma=10.0, beta=8./3, rho=28.0):
      2
----> 3     fig = plt.figure()
      4     ax = fig.add_axes([0, 0, 1, 1], projection='3d')
      5     ax.axis('off')

NameError: name 'plt' is not defined

Using IPython’s interactive function, we can explore how the trajectories behave as we change the various parameters.

[6]:
w = interactive(solve_lorenz, angle=(0.,360.), max_time=(0.1, 4.0),
                N=(0,50), sigma=(0.0,50.0), rho=(0.0,50.0))
display(w)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
~/checkouts/readthedocs.org/user_builds/ipywidgets/conda/latest/lib/python3.6/site-packages/ipywidgets/widgets/interaction.py in update(self, *args)
    254                     value = widget.get_interact_value()
    255                     self.kwargs[widget._kwarg] = value
--> 256                 self.result = self.f(**self.kwargs)
    257                 show_inline_matplotlib_plots()
    258                 if self.auto_display and self.result is not None:

<ipython-input-4-864093adb516> in solve_lorenz(N, angle, max_time, sigma, beta, rho)
      1 def solve_lorenz(N=10, angle=0.0, max_time=4.0, sigma=10.0, beta=8./3, rho=28.0):
      2
----> 3     fig = plt.figure()
      4     ax = fig.add_axes([0, 0, 1, 1], projection='3d')
      5     ax.axis('off')

NameError: name 'plt' is not defined

The object returned by interactive is a Widget object and it has attributes that contain the current result and arguments:

[7]:
t, x_t = w.result
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-7-1964f38448fa> in <module>
----> 1 t, x_t = w.result

TypeError: 'NoneType' object is not iterable
[8]:
w.kwargs
[8]:
{'N': 10,
 'angle': 0.0,
 'max_time': 4.0,
 'sigma': 10.0,
 'beta': 2.6666666666666665,
 'rho': 28.0}

After interacting with the system, we can take the result and perform further computations. In this case, we compute the average positions in \(x), \(y) and \(z).

[9]:
xyz_avg = x_t.mean(axis=1)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-9-ac815f30bff3> in <module>
----> 1 xyz_avg = x_t.mean(axis=1)

NameError: name 'x_t' is not defined
[10]:
xyz_avg.shape
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-da1da28fd348> in <module>
----> 1 xyz_avg.shape

NameError: name 'xyz_avg' is not defined

Creating histograms of the average positions (across different trajectories) show that on average the trajectories swirl about the attractors.

[11]:
plt.hist(xyz_avg[:,0])
plt.title('Average $x(t)$');
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-11-64e431047b89> in <module>
----> 1 plt.hist(xyz_avg[:,0])
      2 plt.title('Average $x(t)$');

NameError: name 'plt' is not defined
[12]:
plt.hist(xyz_avg[:,1])
plt.title('Average $y(t)$');
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-12-f178f527acab> in <module>
----> 1 plt.hist(xyz_avg[:,1])
      2 plt.title('Average $y(t)$');

NameError: name 'plt' is not defined
[ ]: