## Wednesday, February 15, 2017

### Plotting 3D graphs using Python and Tellurium

Originally Posted on  by

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:

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 teimport roadrunnerfrom mpl_toolkits.mplot3d import Axes3Dimport matplotlib as mplimport matplotlib.pyplot as pltimport 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; # Chaosr.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()