我正在嘗試繪制無阻尼和無驅動擺的周期和振幅之間的關系，以應對小角度近似失效時的情況，但是，我的代碼沒有達到我的預期......這是我的代碼，我使用過零法來計算周期：import numpy as npimport matplotlib.pyplot as pltfrom scipy.integrate import solve_ivpfrom itertools import chain# Second order differential equation to be solved:# d^2 theta/dt^2 = - (g/l)*sin(theta) - q* (d theta/dt) + F*sin(omega*t)# set g = l and omega = 2/3 rad per second# Let y[0] = theta, y[1] = d(theta)/dtdef derivatives(t,y,q,F):? ? return [y[1], -np.sin(y[0])-q*y[1]+F*np.sin((2/3)*t)]t = np.linspace(0.0, 100, 10000)#initial conditions:theta0, omega0theta0 = np.linspace(0.0,np.pi,100)q = 0.0? ? ? ?#alpha / (mass*g), resistive termF = 0.0? ? ? ?#G*np.sin(2*t/3)value = []amp = []period = []for i in range (len(theta0)):? ? sol = solve_ivp(derivatives, (0.0,100.0), (theta0[i], 0.0), method = 'RK45', t_eval = t,args = (q,F))? ? velocity = sol.y[1]? ? time = sol.t? ? zero_cross = 0? ? for k in range (len(velocity)-1):? ? ? ? if (velocity[k+1]*velocity[k]) < 0:? ? ? ? ? ? zero_cross += 1? ? ? ? ? ? value.append(k)? ? ? ? else:? ? ? ? ? ? zero_cross += 0? ? if zero_cross != 0:? ? ? ? amp.append(theta0[i])? ? ? # period calculated using the time evolved between the first and last zero-crossing detected? ? ? ? period.append((2*(time[value[zero_cross - 1]] - time[value[0]])) / (zero_cross -1))plt.plot(amp,period)plt.title('Period of oscillation of an undamped, undriven pendulum \nwith varying initial angular displacemnet')plt.xlabel('Initial Displacement')plt.ylabel('Period/s')plt.show()