为了账号安全,请及时绑定邮箱和手机立即绑定

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

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

慕丝7291255 2023-05-23 16:22:24
我正在尝试绘制无阻尼和无驱动摆的周期和振幅之间的关系,以应对小角度近似失效时的情况,但是,我的代码没有达到我的预期......这是我的代码,我使用过零法来计算周期: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()
查看完整描述

1 回答

?
一只名叫tom的猫

TA贡献1906条经验 获得超2个赞

您可以将 的事件机制用于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])

这导致情节

//img1.sycdn.imooc.com//646c77f10001b7b303840289.jpg

查看完整回答
反对 回复 2023-05-23
  • 1 回答
  • 0 关注
  • 71 浏览
慕课专栏
更多

添加回答

举报

0/150
提交
取消
意见反馈 帮助中心 APP下载
官方微信