尝试编译以下程序时收到错误消息。目标是分析非线性摆的行为和一个特定的庞加莱截面。我试图将钟摆行为的数据打印到一个文件中,并将对应于 Poincare 部分的输出打印到另一个文件中。
python shell中列出的错误如下:
Traceback (most recent call last):
File "C:\Python32\Lib\site-packages\visual\examples\JMB312.py", line 39
thetaDotDot = -(g / L) * sin(theta) - q * thetaDot + Fd * sin(Ohm * t)
TypeError: unsupported operand type(s) for /: '_io.TextIOWrapper' and 'float'
我真的不明白为什么 io.TextIOWrapper 会出现在这个错误中,或者如何解决这个问题。任何帮助表示赞赏!
代码:
from visual import *
from visual.graph import *
scene.title = 'Pendulum'
scene.autoscale = False
scene.height = scene.width = 500
pivot = array([0,0,0]) # pivot position of pendulum
scene.center = pivot # center scene on pendulum pivot
framesPerSecond = 400
myWindow1 = gdisplay(x=500, title = "Phase Space", xtitle="Theta", ytitle="ThetaDot", height=500, width=500)
f1 = gdots(gdisplay=myWindow1, color = color.cyan, size = 1)
myWindow2 = gdisplay(title="Poincare Section", xtitle="Theta", ytitle="ThetaDot",y=400, x=500, height=500, width=500)
f2 = gdots(gdisplay=myWindow2, color = color.yellow)
#initial conditions for pendulum 1
g = 9.8 # gravity
m = 1 # mass at end of rod
L = 9.8 # length of rod
theta = 0.2 # initial angle (from vertical)
thetaDot = 0 # start pendulum at rest
thetaDotDot = 0 # acceletation
dt = 0.04 # time step
t = 0 # time
Fd = 1.2 # drive force coefficient
Ohm = 2/3 # oscillatory coefficient for drive force
q = 0.5 # coefficient of friction
#pendulum 1 rod
rod = cylinder(pos=pivot, axis = (L*sin(theta),-L*cos(theta),0), radius = 0.1, color=color.cyan)
f=open('JMB312 data 1.txt','w')
g=open('JMB312 data 2.txt','w')
while t < 100000: #a lot of data
thetaDotDot = -(g / L) * sin(theta) - q * thetaDot + Fd * sin(Ohm * t)
thetaDot = thetaDot + thetaDotDot * dt
theta = theta + thetaDot * dt
time = str(t)
angle = str(theta)
angular_velocity = str(thetaDot)
if t < 750:
f.write(time)
f.write(",")
f.write(angle)
f.write(",")
f.write(angular_velocity)
f.write("\n")
f1.plot(pos=(theta,thetaDot))
if -dt**2/4 < Fd * sin(Ohm * t) < dt**2/4: #print poincare section data to separate file
g.write(time)
g.write(",")
g.write(angle)
g.write(",")
g.write(angular_velocity)
g.write("\n")
f2.plot(pos=(theta,thetaDot))
if theta > 2 * pi:
theta = theta - 2 * pi
if theta < -2 * pi:
theta = theta + 2 * pi
rod.axis = (L*sin(theta), -L*cos(theta), 0)
t = t+dt
g.close()
f.close()