python - 大角摆图未按预期显示
问题描述
我正在尝试绘制一个无阻尼和无驱动摆的周期和振幅之间的关系,以便在小角度近似失效时,但是,我的代码没有达到我的预期......
我想我应该期待一个严格增加的图表,如本视频所示: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()
解决方案
您可以将事件机制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])
这导致情节
推荐阅读
- java - 在数据集中转换 Flink CEP 表时出错
目的
- python - 绘制一个显示直到下次单击的形状
- javascript - jQuery Ajax 发布 JSON - 响应
- angular - 离子图标启动时不充电
- javascript - jquery中的对象未在提交时设置
- elasticsearch - 当我尝试编辑 elasticsearch.yml 文件时,Elasticsearch 未运行
- angular - 如何使用路由制作活动按钮?
- android - Android R8 - 缺少类
- django - 单击两次链接时出现Django URL页面错误
- ckeditor - 如何将html添加到图像插件ckeditor 4的标签上传