Python程序设计期末作业

球-弹簧物理模型模拟程序一、编写目的1.模拟一个非线性动力系统，体会系统参数的混沌现象。2.熟悉VPython的使用。二、物理模型描述模型空间的四周及顶面为刚性壁，底为一有质量的刚性平板，平板由一根轻弹簧支撑。模型空间内有数个刚性球作自由运动。球的初始位置在空间内均匀分布，初始速度满足正态分布。所有碰撞均为完全弹性碰撞。并假设平板只能作Z方向一维运动。三、设计思想由物理学中的能量守恒和动量守恒推导碰撞前后物体的速度变化。每次迭代先使物体状态按时间小量dt发展，再判断是否有碰撞发生，若有则由前面推导的公式改变相关物体的速度。以适当的速率重复上述步骤，再用visual库绘制出对应图形，则得到这个物理模型的模拟演示动画。四、使用方法直接运行程序，则开始模拟并实时绘制。也可修改程序开始处的几个全局变量，改变模型的参数。五、参考资料VPython官方文档六、说明在Python2.7下调试通过。本程序中没有一行代码直接来自网络、书籍等处，参考资料仅为VPython官方文档。附截图及源代码。程序运行效果截图fromvisualimport*fromrandomimportrandomfrommathimportsin,cos,sqrt,log,pi#物理模型的参数以全局变量形式定义g=9.81#重力加速度k=10#弹簧胡克系数m1=1#球质量m2=10#平板质量N=10#球数R=1#球半径#物体基类，提供设置、返回位置、速度信息，及迭代一个时间小量dt的方法classobject:defsetVel(self,v):self.v=vdefsetPos(self,p):self.p=pdefgetPos(self):returnself.pdefgetVel(self):returnself.vdefgo(self,dt):self.v+=self.a*dtself.p+=self.v*dt#球，继承于物体类classC_Ball(object):def__init__(self,p,v):self.p=pself.v=vself.a=vector(0,0,-g)#平板，继承于物体类，有自己的迭代方法，即加速度与位置有关classC_Board(object):def__init__(self):self.p=vector(0,0,0)self.v=vector(0,0,0)self.a=vector(0,0,0)defgo(self,dt):self.a=self.p*-kself.v+=self.a*dtself.p+=self.v*dt#生成随机位置向量（模型空间内均匀分布）defrandPos():x=(random()-0.5)*2*(15-R-0.5)y=(random()-0.5)*2*(15-R-0.5)z=random()*(30-R*2-0.5)+(R+0.5)returnvector(x,y,z)#生成三维标准正态分布向量，使用Box-Muller公式defrandVel():u1,u2=random(),random()x=sqrt(-2*log(u1))*cos(2*pi*u2)y=sqrt(-2*log(u1))*sin(2*pi*u2)u1,u2=random(),random()z=sqrt(-2*log(u1))*cos(2*pi*u2)returnvector(x,y,z)#物理模型，含有一个平板和N个球classC_System:#放置球，并检查冲突def__init__(self):self.ball=[]foriinrange(N):self.ball.append(C_Ball(randPos(),randVel()*3))j=0whileji:#若有两个球干涉则重新放置n=self.ball[i].getPos()-self.ball[j].getPos()ifmag(n)=2*R:self.ball[-1]=C_Ball(randPos(),randVel()*3)j=-1j+=1self.board=C_Board()#迭代一个时间小量dtdefgo(self,dt):foriinrange(N):self.ball[i].go(dt)self.board.go(dt)foriinrange(N):#检查与平板的碰撞ifself.ball[i].getPos().z-R=self.board.getPos().z+0.5:self.ball[i].setPos(self.ball[i].getPos()+(self.board.getPos().z-self.ball[i].getPos().z+1.5)*vector(0,0,1))v1=self.ball[i].getVel()v2=self.board.getVel()v1z=v1.zv2z=v2.zself.ball[i].setVel(v1-vector(0,0,v1z)+vector(0,0,1)*(-v1z*m2+m1*v1z+2*m2*v2z)/(m2+m1))self.board.setVel(vector(0,0,(m2*v2z-v2z*m1+2*m1*v1z)/(m2+m1)))#检查与上壁的碰撞ifself.ball[i].getPos().z+R=30-0.5:self.ball[i].go(-dt)v=self.ball[i].getVel()n=vector(0,0,1)l=dot(n,v)*2*nself.ball[i].setVel(v-l)#检查与前后壁的碰撞ifabs(self.ball[i].getPos().x)+R=15:self.ball[i].go(-dt)v=self.ball[i].getVel()n=vector(1,0,0)l=dot(n,v)*2*nself.ball[i].setVel(v-l)#检查与左右壁的碰撞ifabs(self.ball[i].getPos().y)+R=15:self.ball[i].go(-dt)v=self.ball[i].getVel()n=vector(0,1,0)l=dot(n,v)*2*nself.ball[i].setVel(v-l)#检查两球之间的碰撞forjinrange(i):n=self.ball[i].getPos()-self.ball[j].getPos()ifmag(n)=2*R:self.ball[i].go(-dt)self.ball[j].go(-dt)n=norm(n)v1=self.ball[i].getVel()v2=self.ball[j].getVel()l1=dot(n,v1)*2*nl2=dot(n,v2)*2*nself.ball[i].setVel(v1-l1)self.ball[j].setVel(v2-l2)#返回第i个球的位置向量defgetBallPos(self,i):returnself.ball[i].getPos()#返回平板位置向量defgetBoardPos(self):returnself.board.getPos()mySys=C_System()scene=display(title='Animation',x=300,y=100,width=480,height=640,center=(0,0,0),autocenter=False,range=(40,40,60),up=(0,0,1),forward=(-8,-1,-2))B=[]#球数组foriinrange(N):B.append(sphere(pos=mySys.getBallPos(i),radius=R,material=materials.plastic))P=box(pos=mySys.getBoardPos(),length=28.5,height=28.5,width=1,material=materials.wood)#平板H=helix(pos=(0,0,-30),axis=(0,0,1),coils=10,radius=7.5,length=30,material=materials.chrome)#弹簧#生成三面壁和上壁WallL=box(pos=(0,-15,0),length=30,height=1,width=59,material=materials.rough)WallR=box(pos=(0,15,0),length=30,height=1,width=59,material=materials.rough)WallF=box(pos=(-15,0,0),length=1,height=30,width=59,material=materials.rough)WallU=box(pos=(0,0,30),length=31,height=31,width=1,material=materials.rough)whileTrue:foriinrange(100):mySys.go(4e-4)#迭代100次，每次0.4毫秒，则系统状态更新40毫秒#根据物理模型重新绘制物体foriinrange(N):B[i].pos=mySys.getBallPos(i)P.pos=mySys.getBoardPos()H.length=P.pos.z+30rate(50)