Analogmachine Blog: Plotting 3D graphs using Python and Tellurium

Originally Posted on February 15, 2017 by hsauro

As an example I wanted to show how one could plot a 3D phase plot. A great example to use for this is the Lorenz Attractor. This system is interesting because it displays chaotic behavior. The differential equations for the system are given by the following three:

${\begin{aligned}{\frac {\mathrm {d} x}{\mathrm {d} t}}&=\sigma (y-x),\\{\frac {\mathrm {d} y}{\mathrm {d} t}}&=x(\rho -z)-y,\\{\frac {\mathrm {d} z}{\mathrm {d} t}}&=xy-\beta z.\end{aligned}}$

Different values for the parameters, sigma, rho and beta, lead to different behaviors. The Python script below that uses Tellurium plots a 3D phase plot:

# The Lorenz Model
 
# The Lorenz system of equations is the classic set of
# equations which exhibit chaotic dynamics. The equations
# represent a highly simplified model of 2D thermal
# convection known as Rayleigh-Benard convection;
 
# The three variables, x, y and z represent the
# convective overturning, the horizontal temperature
# and vertical temperature variation respectively;
 
## The parameters a, b and c represent approximately
# the energy losses within the fluid due to viscosity,
# the temperature difference between the two plates
# which separate the fluid and the ratio of the vertical
# height of the fluid layer to the horizontal extent
# of the convective rolls respectively;
import tellurium as te
import roadrunner
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
 
r = te.loada ('''
// Define the model in terms of three differential equations;
$s -> x; -a*(x - y);
$s -> y; -x*z + b*x - y;
$s -> z; x*y - c*z;
 
a = 10; b = 28; c = 2.67;
''')
 
# Uncomment/comment the following lines to obtain behaviour you wish to observe;
 
# Non-chaotic oscillations;
#r.x = 5.43027; r.y = 26.74167; r.z = 46.0098;
#r.a = 4; r.b = 91.9; r.c = 2.67;
 
# Single oscillation;
#r.x = 32.82341; r.y = -44.84651; r.z = 258.02389;
#r.a = 4; r.b = 191.9; r.x = 2.67;
 
# Stable steady-state;
#r.x = 0.1; r.y = 0.1; r.z = 0.1;
#r.a = 10; r.b = 10; r.c = 2.67;
 
# Chaos
r.x = 0.96259; r.y = 2.07272; r.z = 18.65888;
 
m = r.simulate (0, 40, 3000, ['x', 'y', 'z']);
 
fig = plt.figure(figsize=(8,8))
ax = fig.gca(projection='3d')
 
#theta = np.linspace(-4 * np.pi, 4 * np.pi, 100)
ax.plot(m[:,0], m[:,1], m[:,2], label='Lorenz Attractor')
ax.legend()
 
plt.show()

Analogmachine Blog

Pages

Wednesday, February 15, 2017

Plotting 3D graphs using Python and Tellurium

No comments: