大角度摆线图未按预期显示

问题描述

当小角度近似分解失败时，我试图绘制一个无阻尼且无驱动的摆的周期与幅度之间的关系，但是，我的代码没有达到我的期望...

我想我应该期待视频中显示的严格增加的图表：https://www.youtube.com/watch?v=34zcw_nNFGU

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

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

enter image description here

解决方法

您可以将solve_ivp的事件机制用于此类任务，它是针对此类“简单”情况而设计的

def halfperiod(t,y): return y[1]
halfperiod.terminal=True  #  stop when root found
halfperiod.direction=1    #  find sign changes from negative to positive

for i in range (1,len(theta0)): # amp==0 gives no useful result
    sol = solve_ivp(derivatives,(0.0,100.0),(theta0[i],0.0),method = 'RK45',events =(halfperiod,) )
    if sol.success and len(sol.t_events[-1])>0:
        period.append(2*sol.t_events[-1][0]) # the full period is twice the event time
        amp.append(theta0[i])

这导致了情节