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
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 Get Luis Morales’s stories in your inbox Join Medium for free to get updates from this writer.
Enter your email Subscribe #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() Press enter or click to view image in full size
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.
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()
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!