大角度摆图没有按预期显示

Question

我正在尝试绘制无阻尼和无驱动摆的周期和振幅之间的关系，以应对小角度近似失效时的情况，但是，我的代码没有达到我的预期......

这是我的代码，我使用过零法来计算周期：

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp

from 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)/dt

def 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, omega0

theta0 = np.linspace(0.0,np.pi,100)

q = 0.0       #alpha / (mass*g), resistive term

F = 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()

Answer 1

您可以将 的事件机制用于solve_ivp此类任务，它专为此类“简单”情况而设计def halfperiod(t,y): return y[1]halfperiod.terminal=True&nbsp; #&nbsp; stop when root foundhalfperiod.direction=1&nbsp; &nbsp; #&nbsp; find sign changes from negative to positivefor i in range (1,len(theta0)): # amp==0 gives no useful result&nbsp; &nbsp; sol = solve_ivp(derivatives, (0.0,100.0), (theta0[i], 0.0), method = 'RK45', events =(halfperiod,) )&nbsp; &nbsp; if sol.success and len(sol.t_events[-1])>0:&nbsp; &nbsp; &nbsp; &nbsp; period.append(2*sol.t_events[-1][0]) # the full period is twice the event time&nbsp; &nbsp; &nbsp; &nbsp; amp.append(theta0[i])这导致情节