Python and Physics: Lorenz and Rossler Systems

Building upon my previous article on the Runge-Kutta method, we are going to explore some applications of the Runge-Kutta method to graph the famous Lorenz and Rossler systems.

Lorenz

The Lorenz system are a group of equations that show chaotic behavior for certain parameters, which most notably produce an interesting graph. The equations that we are going to work with are the following

Lorenz Equations

These will be the equations that we will plugging in to our derivative function so that we can apply the RK4 (Runge-Kutta Four-Step) method. First, we start off by importing the necessary packages.

import numpy as np #For arrays
import matplotlib.pyplot as plt #For plotting

Then, we start with the variables we will be using for these equations.

sigma = 10.0 #Variable for dx/dt
rho = 28.0 #Variable for dy/dt
beta = 8/3 #Variable for dz/dt
t = 0 #Starting time
tf = 40 #Ending time
h = 0.01 #Step size for RK4
#These variables were used since they can be easily found for comparison, for example, in Wikipedia under the Lorenz system article

With our variables set, we now move on to our derivative function.

#Derivative function to work with RK4 loop
def derivative(r,t):
    x = r[0]
    y = r[1]
    z = r[2]
    return np.array([sigma * (y - x), x * (rho - z) - y, (x * y) - (beta * z)])

A quick refresher for the derivative function. The r variable is for the initial conditions set inside an array. We assign each variable with the values of the of the array, three in this case, and then the function will return the differential equations that were pictured at the beginning, but now with initial conditions. Each time loop goes through an iteration, the function will continue updating with new values, which will be stored in the following arrays.

time = np.array([]) #Empty time array to fill for the x-axis
x = np.array([]) #Empty array for x values
y = np.array([]) #Empty array for y values
z = np.array([]) #Empty array for z values
r = np.array([1.0, 1.0, 1.0]) #Initial conditions array

Now, with our main pieces set, we just go ahead and apply everything inside the RK4 loop and we end up with

while (t <= tf ):
        #Appending values to graph
        time = np.append(time, t)
        z = np.append(z, r[2])
        y = np.append(y, r[1])
        x = np.append(x, r[0])
        #RK4 Step method
        k1 = h*derivative(r,t)
        k2 = h*derivative(r+k1/2,t+h/2)
        k3 = h*derivative(r+k2/2,t+h/2)
        k4 = h*derivative(r+k3,t+h)
        r += (k1+2*k2+2*k3+k4)/6
        #Updating time value with step size
        t = t + h

After running the loop, our arrays should be filled with an equal amount of values, so now we can just go ahead and plot these values so we can see how they produce the graphs below, using the following code.

#Multiple graph plotting
fig, (ax1,ax2,ax3) = plt.subplots(1,3, figsize = (15, 5))
ax1.plot(x, y)
ax1.set_title("X & Y")
ax2.plot(x, z)
ax2.set_title("X & Z")
ax3.plot(y, z)
ax3.set_title("Y & Z")
plt.show()

Lorenz System graphs, plotted against each variable

Now those are some interesting graphs! Playing around with different values for our Greek variables, we can see changes in the behavior for our system, with some values producing complex graphs.

Rossler

Much like the Lorenz system, we apply our differential equations in our derivative function and run it through the RK4 loop to produce a visual representation of the Rossler System, by using the following equations.

Rossler Equations

With these equations, we just change the return line in our derivative function equation to these and we are all set. Below, you will find the complete code for this one, as well as the graphs produced.

import numpy as np
import matplotlib.pyplot as plt
a = 0.2
b = 0.2
c = 5.7
t = 0
tf = 100
h = 0.01
def derivative(r,t):
    x = r[0]
    y = r[1]
    z = r[2]
    return np.array([- y - z, x + a * y, b + z * (x - c)])
time = np.array([])
x = np.array([])
y = np.array([])
z = np.array([])
r = np.array([0.1, 0.1, 0.1])
while (t <= tf ):
    
        time = np.append(time, t)
        z = np.append(z, r[2])
        y = np.append(y, r[1])
        x = np.append(x, r[0])
        
        k1 = h*derivative(r,t)
        k2 = h*derivative(r+k1/2,t+h/2)
        k3 = h*derivative(r+k2/2,t+h/2)
        k4 = h*derivative(r+k3,t+h)
        r += (k1+2*k2+2*k3+k4)/6
        
        t = t + h
fig, (ax1,ax2,ax3) = plt.subplots(1,3, figsize = (15, 5))
ax1.plot(x, y)
ax1.set_title("X & Y")
ax2.plot(x, z)
ax2.set_title("X & Z")
ax3.plot(y, z)
ax3.set_title("Y & Z")
plt.show()

Rossler System graphs, plotted against each variable

Much like the Lorenz system, the Rossler system also exhibits complex behavior depending on the variables used for the equation.

Overview

We can see how useful the Runge-Kutta Method is for programming equations that use ODE differential equations. While these are more for exploring physics concepts, this method can also be used to help with other more practical simulations. Futhermore, I highly suggest checking out the Wikipedia pages on these two systems since they actually show you how program these equations in Python but in a 3D space!