maths with Python 2: Rössler system

After installing python and the packages to code some math, it’s time to use them.

Open the GUI and File>New Window. This time instead of writting the comands directly on the comand window, we are going to make a script with all the operations.

We are going to solve numerically the Rössler equation system and learn some tricks in the journey. The system is composed of this set of eaquations which describe the evolution of variables x\quad y \quad z.

\left\{ \begin{matrix} \frac{\mathrm{d}x}{\mathrm{d}t}=-y-z \\ \frac{\mathrm{d}y}{\mathrm{d}t}=x+ay \\  \frac{\mathrm{d}z}{\mathrm{d}t}=b+z(x-c)) \end{matrix}\right.

For parameters

a=0.13 \quad b=0.2 \quad c=6.5

and initial conditions

x(0)=0 \quad y(0)=0 \quad z(0)=0

using Maple with it’s differential plotting tools we can see how the solutions looks like:

roslerr2The reason to use these values it’s just fun, no especial interest on them. Only on parameter a, because for 0<a<2 the origin is an unstable point.

Ok, so now turning the equations into numerical equations, how to do it?

In this case we are going to use the most easiest and simple way, the forward Euler method. This is the easiest and more simple method because you don’t need to solve any equation and the solution is built in a recursive manner.

The Euler method consist in using these definition of derivative as the approximation of the derivative

We have

\lim_{h\rightarrow 0} { \frac{ f(t+h) - f(t) }{ h} } = \frac{\mathrm{d} f(t)}{ \mathrm{d}t}

and we approximate:

{\frac{ f(t+h) - f(t) }{ h}} \simeq \frac{\mathrm{d} f(t)}{ \mathrm{d}t}

and the approximation is better the smaller h is.

So in our case we want to find

x(t) \quad y(t) \quad z(t)

which satisfy the Rössler system. We start by splitting the time into equally spaced steps t_n such

t_{n+1}-t_n=h

and we name

x_n=x(t_n)

Using the Euler approximation we can rewrite the Rössler system as the next recursion relationship:

\left\{ \begin{matrix} \frac{x_{n+1}-x_n}{h}=-y_n-z_n \\ \frac{y_{n+1}-y_n}{h}=x_n+ay_n \\  \frac{z_{n+1}-z_n}{h}=b+z_n(x_n-c)) \end{matrix}\right.

We can now turn this equations into it’s explicit form

\left\{ \begin{matrix} x_{n+1}=x_n+h(-y_n-z_n) \\ y_{n+1}=y_n+h(x_n+ay_n) \\  z_{n+1}=z_n+h(b+z_n(x_n-c))) \end{matrix}\right.

and since we know the nitial conditions

x_0 \quad y_0 \quad z_0

we can find easily a numerical solution to the system.

Let’s now see how to do that with Python. Get back to the open scripting window. This is the code to do what we need:

##This is a simple code to solve Rössler equation system in python and plot solutions.
##Made by Héctor Corte-León leo_corte@yahoo.es on 23/02/2013
##The numerical method used is first order forward Euler method.

#############################################

##First we load some packages.
from numpy import *
from matplotlib import *
from scipy import *
from pylab import figure, show, setp
from mpl_toolkits.mplot3d import Axes3D

#We define a function which is going to be the recursive function.
def num_rossler(x_n,y_n,z_n,h,a,b,c):
    x_n1=x_n+h*(-y_n-z_n)
    y_n1=y_n+h*(x_n+a*y_n)
    z_n1=z_n+h*(b+z_n*(x_n-c))   
    return x_n1,y_n1,z_n1

#Now we prepare some variables
#First the parameters
a=0.13
b=0.2
c=6.5

#Them the time interval and the step size
t_ini=0
t_fin=32*pi
h=0.0001
numsteps=int((t_fin-t_ini)/h)

#using this parameters we build the time.
t=linspace(t_ini,t_fin,numsteps)
#And the vectors for the solutions
x=zeros(numsteps)
y=zeros(numsteps)
z=zeros(numsteps)

#We set the initial conditions
x[0]=0
y[0]=0
z[0]=0

#This is the main loop where we use the recursive system to obtain the solution
for k in range(x.size-1):
    #We use the previous point to generate the new point using the recursion
    [x[k+1],y[k+1],z[k+1]]=num_rossler(x[k],y[k],z[k],t[k+1]-t[k],a,b,c)
   
#Now that we have the solution in vectors t,x,y,z is time to plot them.

#We create a figure and 4 axes on it. 3 of the axes are going to be 2D and the fourth one is a 3D plot.

fig = figure()
ax1 = fig.add_axes([0.1, 0.7, 0.4, 0.2])
ax2 = fig.add_axes([0.1, 0.4, 0.4, 0.2])
ax3 = fig.add_axes([0.1, 0.1, 0.4, 0.2])
ax4 = fig.add_axes([0.55, 0.25, 0.35, 0.5],projection='3d')

#And we add vectors to each plot
ax1.plot(t, x,color='red',lw=1,label='x(t)')
ax1.set_xlabel('t')
ax1.set_ylabel('x(t)')
ax1.legend()
ax1.axis((t_ini,t_fin,min(x),max(x)))

ax2.plot(t, y,color='green',lw=1,label='y(t)')
ax2.set_xlabel('t')
ax2.set_ylabel('y(t)')
ax2.legend()
ax2.axis((t_ini,t_fin,min(y),max(y)))

ax3.plot(t, z,color='blue',lw=1,label='z(t)')
ax3.set_xlabel('t')
ax3.set_ylabel('z(t)')
ax3.legend()
ax3.axis((t_ini,t_fin,min(z),max(z)))

ax4.plot(x, y,z,color='black',lw=1,label='Evolution(t)')
ax4.set_xlabel('x(t)')
ax4.set_ylabel('y(t)')
ax4.set_zlabel('z(t)')
ax4.set_title('Evolution')
#When finished we show the figure with all the plots.
show()

Once you have the code written and understand what it is doing, press F5 to execute the script. After a short time a figure like this one will appear. And voilá!

pythonross

On next tutorial we are going to construct the bifurcation diagram for x when we sweep parameter b. Enjoy this meanwhile!
beauty

Advertisements

9 thoughts on “maths with Python 2: Rössler system”

    1. Hi Kimi. Sorry for not giving you an answer before. The most tricky thing is to decide when it has “landed”. I don’t know of any condition it has to full fill, so I cannot tell you how to check when it has gone into a periodic orbit in the attractor, but I can tell you how to plot it. If you found for instance that after 60 time units the solution is in the attractor, you can put this code instead the one on the post. Notice the change in the plot command.

      #And we add vectors to each plot
      tini=round(60/h); #60 suppose to be for isntance the moment it “landed”
      ax1.plot(t[tini:size(t)], x[tini:size(t)],color=’red’,lw=1,label=’x(t)’)
      ax1.set_xlabel(‘t’)
      ax1.set_ylabel(‘x(t)’)
      ax1.legend()
      ax1.axis((t_ini,t_fin,min(x),max(x)))

      ax2.plot(t[tini:size(t)], y[tini:size(t)],color=’green’,lw=1,label=’y(t)’)
      ax2.set_xlabel(‘t’)
      ax2.set_ylabel(‘y(t)’)
      ax2.legend()
      ax2.axis((t_ini,t_fin,min(y),max(y)))

      ax3.plot(t[tini:size(t)], z[tini:size(t)],color=’blue’,lw=1,label=’z(t)’)
      ax3.set_xlabel(‘t’)
      ax3.set_ylabel(‘z(t)’)
      ax3.legend()
      ax3.axis((t_ini,t_fin,min(z),max(z)))

      ax4.plot(x[tini:size(t)], y[tini:size(t)],z[tini:size(t)],color=’black’,lw=1,label=’Evolution(t)’)
      ax4.set_xlabel(‘x(t)’)
      ax4.set_ylabel(‘y(t)’)
      ax4.set_zlabel(‘z(t)’)
      ax4.set_title(‘Evolution’)
      #When finished we show the figure with all the plots.
      show()

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s