I have the following python code:
import matplotlib.pyplot as plt
import numpy as np
class bifurcation_diagram(object):
def __init__(self):
self.omega = []
self.theta = []
self.dt = (2 * np.pi / (2.0 / 3.0)) / 600
self.time = []
self.theta_in_bifurcation_diagram = []
self.F_D = np.arange(1.35,1.5,0.001)
self.theta_versus_F_D = []
def calculate(self):
l = 9.8
g = 9.8
q = 0.5
Omega_D = 2.0 / 3.0
for f_d in self.F_D:
self.omega.append([0])
self.theta.append([0.2])
self.time.append([0])
for i in range(600000):
k1_theta = self.dt * self.omega[-1][-1]
k1_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1]) - q * self.omega[-1][-1] + f_d * np.sin(Omega_D * self.time[-1][-1]))
k2_theta = self.dt * (self.omega[-1][-1] + 0.5 * k1_omega)
k2_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1] + 0.5 * k1_theta) - q * (self.omega[-1][-1] + 0.5 * k1_omega) + f_d * np.sin(Omega_D * (self.time[-1][-1] + 0.5 * self.dt)))
k3_theta = self.dt * (self.omega[-1][-1] + 0.5 * k2_omega)
k3_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1] + 0.5 * k2_theta) - q * (self.omega[-1][-1] + 0.5 * k2_omega) + f_d * np.sin(Omega_D * (self.time[-1][-1] + 0.5 * self.dt)))
k4_theta = self.dt * (self.omega[-1][-1] + k3_omega)
k4_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1] + k3_theta) - q * (self.omega[-1][-1] + k3_omega) + f_d * np.sin(Omega_D * (self.time[-1][-1] + self.dt)))
temp_theta = self.theta[-1][-1] + (1.0 / 6.0) * (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta)
temp_omega = self.omega[-1][-1] + (1.0 / 6.0) * (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega)
while temp_theta > np.pi:
temp_theta -= 2 * np.pi
while temp_theta < -np.pi:
temp_theta += 2 * np.pi
self.omega[-1].append(temp_omega)
self.theta[-1].append(temp_theta)
self.time[-1].append(self.dt * i)
for i in range(500,1000):
n = i * 600
self.theta_in_bifurcation_diagram.append(self.theta[-1][n])
self.theta_versus_F_D.append(f_d)
def show_results(self):
plt.plot(self.theta_versus_F_D,self.theta_in_bifurcation_diagram,'.')
plt.title('Bifurcation diagram' + '\n' + r'$\theta$ versus $F_D$')
plt.xlabel(r'$F_D$')
plt.ylabel(r'$\theta$ (radians)')
plt.xlim(1.35,1.5)
plt.ylim(1,3)
plt.show()
bifurcation = bifurcation_diagram()
bifurcation.calculate()
bifurcation.show_results()
I wish to make it more compact and make a coloured plot, along with increasing its efficiency. Any help in this regard would be truly beneficial.
I was expecting the code to run fast but it takes more than fifteen minutes to run.
Python is not made for such a code because it is generally interpreted using CPython (it is meant to be done for glue codes and scripts) and CPython performs (nearly) no optimizations of the code so repeated expressions are recomputed. The usual way to speed up such a code is to use Numpy functions operating on large arrays (not lists). THis is called vectorization. That being said, this code is pretty expensive and using Numpy may not be enough and it is also a non-negligible work to use Numpy here. An alternative solution is to use Numba so to compile this code (thanks to a JIT compiler). Note that Numba is meant to speed-up codes using mostly Numpy and basic numerical computations (not general-purpose codes like plotting).
The first thing to do is to get ride of lists and replace them by Numpy arrays. Arrays can be preallocated at the right size so to avoid the expensive calls to append
. Then, the computing function should be moved away from the class so for Numba to be easily used. Then, the two while
loops can be replaced by np.remainder
. Finally, you can compute each row of the arrays in parallel using multiple threads.
Here is the resulting code (barely tested):
import matplotlib.pyplot as plt
import numpy as np
import numba as nb
@nb.njit('(float64[::1], float64)', parallel=True)
def compute_numba(F_D, dt):
l = 9.8
g = 9.8
q = 0.5
Omega_D = 2.0 / 3.0
size = 600000
loop2_start, loop2_end = 500, 1000
loop2_stride = 600
omega = np.empty((F_D.size, size+1), dtype=np.float64)
theta = np.empty((F_D.size, size+1), dtype=np.float64)
time = np.empty((F_D.size, size+1), dtype=np.float64)
theta_in_bifurcation_diagram = np.empty((loop2_end-loop2_start) * F_D.size, dtype=np.float64)
theta_versus_F_D = np.empty((loop2_end-loop2_start) * F_D.size, dtype=np.float64)
for i in nb.prange(F_D.size):
f_d = F_D[i]
omega[i, 0] = 0.0
theta[i, 0] = 0.2
time[i, 0] = 0.0
for j in range(size):
cur_omega = omega[i,j]
cur_theta = theta[i,j]
cur_time = time[i,j]
tmp = np.sin(Omega_D * (cur_time + 0.5 * dt))
k1_theta = dt * cur_omega
k1_omega = dt * ((-g / l) * np.sin(cur_theta) - q * cur_omega + f_d * np.sin(Omega_D * cur_time))
k2_theta = dt * (cur_omega + 0.5 * k1_omega)
k2_omega = dt * ((-g / l) * np.sin(cur_theta + 0.5 * k1_theta) - q * (cur_omega + 0.5 * k1_omega) + f_d * tmp)
k3_theta = dt * (cur_omega + 0.5 * k2_omega)
k3_omega = dt * ((-g / l) * np.sin(cur_theta + 0.5 * k2_theta) - q * (cur_omega + 0.5 * k2_omega) + f_d * tmp)
k4_theta = dt * (cur_omega + k3_omega)
k4_omega = dt * ((-g / l) * np.sin(cur_theta + k3_theta) - q * (cur_omega + k3_omega) + f_d * np.sin(Omega_D * (cur_time + dt)))
temp_theta = cur_theta + (1.0 / 6.0) * (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta)
temp_omega = cur_omega + (1.0 / 6.0) * (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega)
temp_theta = np.remainder(temp_theta + np.pi, 2 * np.pi) - np.pi
omega[i,j+1] = temp_omega
theta[i,j+1] = temp_theta
time[i,j+1] = dt * j
for k, j in enumerate(range(loop2_start, loop2_end)):
n = j * loop2_stride
offset = (loop2_end - loop2_start) * i + k
theta_in_bifurcation_diagram[offset] = theta[i,n]
theta_versus_F_D[offset] = f_d
return (omega, theta, time, theta_in_bifurcation_diagram, theta_versus_F_D)
class bifurcation_diagram(object):
def __init__(self):
self.omega = np.array([])
self.theta = np.array([])
self.dt = (2 * np.pi / (2.0 / 3.0)) / 600
self.time = np.array([])
self.theta_in_bifurcation_diagram = np.array([])
self.F_D = np.arange(1.35,1.5,0.001)
self.theta_versus_F_D = np.array([])
def calculate(self):
self.omega, self.theta, self.time, self.theta_in_bifurcation_diagram, self.theta_versus_F_D = compute_numba(self.F_D, self.dt)
def show_results(self):
plt.plot(self.theta_versus_F_D,self.theta_in_bifurcation_diagram,'.')
plt.title('Bifurcation diagram' + '\n' + r'$\theta$ versus $F_D$')
plt.xlabel(r'$F_D$')
plt.ylabel(r'$\theta$ (radians)')
plt.xlim(1.35,1.5)
plt.ylim(1,3)
plt.show()
bifurcation = bifurcation_diagram()
bifurcation.calculate()
bifurcation.show_results()
On my i5-9600KF processor with 6 cores, this code takes 1.07 secondes rather than about 19 min 20 for the initial one! As a result, this new code is about 1150 times faster!
Besides, I also find the resulting code easier to read.
I advise you to check the results, although they look fine at first glance. Indeed, here is the resulting graph I got so far:
Collected from the Internet
Please contact [email protected] to delete if infringement.
Comments