Increasing the efficiency of a python code by using arrays

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: