我有以下的python代码:
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()
我希望把它做得更紧凑，做一个彩色的情节，同时提高它的效率.在这方面的任何帮助都将是真正有益的.
我原以为代码会运行得很快，但它需要15分钟以上才能运行.
