4.彗星轨迹图
彗星轨迹图是最简单的动画,可以展现质点运动的过程,其使用格式为: comet(x,y) comet(x,y,p) comet3(x,y,z) comet3(x,y,z,p) 二维平面彗星轨迹,x,y的含义与plot命令中x,y含义相同,绘出的轨迹慧长为p*length(x),0
?n?1?设空气阻力与抛体运动速度的n次方成正比,F??bvv,b为正值常量.其运动微分方程为
????b(x?2?y?2)(n?1)2x?,my????b(x?2?y?2)(n?1)2y? mx令
?,y(3)?y,y(4)?y?,则化为 y(1)?x,y(2)?xby(2)[y(2)2?y(4)2](n?1)2 mb?(4)??9.8?y(4)[y(2)2?y(4)2](n?1)2 ?(3)?y(4)yym?(2)???(1)?y(2),yy编写微分方程函数文件znxpfun.m如图39;再编写解微分方程的主程序znxp.m如图41;运行结果如图40.
图39 图40
图41
图42
注意例子中第4,5,6行中的命令.comet语句在例6图26中就已经用过了,从现在的图上看不出动态效果,
13
读者自己编程运行后就可以看到了.
把例9图35中第7行plot3改为comet3(x,y,z,p)就是三维彗星轨迹图了,如图42.程序中解微分方程的步长由0.1改为0.002,是为了使彗星轨迹图运行的速度减慢.由于不同计算机计算速度不同,可调整步长使图形运行速度适当.
§3 一维非线性振动
一、一维振动
一维振动一般微分方程形式为
???FR(x?)?F(x)?F0cos?t (3.1) mx二、一维线性振动
在恢复力和阻力满足线性关系条件下,一维线性振动的微分方程为
????x??kx?F0cos?t mx引入2???m,?0?km,f?F0m,则
2???2?x???02x?fcos?t (3.2) x下面用数值计算的方法对一维线性振动做一点讨论.主要是通过x?t图和x?运动情况.
依据(3.2)式编写微分方程函数文件ywzdfun.m如图43;再编写解微分方程的主程序ywdz.m如图44.在主程序中改变参量,⑴令b=0,w0=1,f=0,w=2;运行,得简谐振动的x?t图和相图,如图45和fl.46所示.⑵令b=0.1,w0=1,f=0,w=2;运行,得弱阻尼振动的x?t图和相图,如图47和fl.48所示.⑶令b=0.1,w0=1,f=0.5,w=2;运行,得受迫振动的x?t图和相图,如图49和fl.50所示;再令w=0.99;运行,得共振状态的x?t图和相图,如图51和fl.52所示.
?图(即相图)x,去了解质点的
图43 图44
图45 图46 图47 图48
图49 图50 图51 图52
读者还可以继续改变参量,比如对阻尼振动情况,改变b的值就可以得到临界阻尼和过阻尼相应的图形.
14
三、用数值计算和相图研究大幅度单摆的运动
本节研究最大摆角不受限制的单摆的运动, 单摆由质点和刚性轻杆构成,它的运动微分方程为
????2sin??0 (3.3) ?0这个非线性方程是可积的,它的解可用第一类椭圆积分表示;现在,用数值计算和相图方法研究它.
依据(3.3)式编写微分方程函数文件djdbfun.m如图53;再编写解微分方程并画x?t图的主程序djdb1.m如图54;运行,计算结果如图55.
图53 图54
图55 图56
图57 图58
编写解微分方程并画相图的主程序djdb2.m如图56;运行,计算结果如图57.
最后,编写解微分方程,作快速傅里叶变换并画功率谱的主程序djdb3.m如图58;运行,计算结果如图59. 由图55可见,振幅(最大摆角)越大,曲线与简谐振动正弦曲线的差别越显著.同时,周期随振幅增加而变大,当振幅接近?时,周期迅速增大(实际上,振幅??时,周期??);图60为由数值计算画出的,振幅不太大时的“周期-振幅”曲线.
由图57可见,相轨迹分为两类,它们之间存在相应于总能量E?mgl的
分界线(以悬点为重力势能零点);分界线之内的相轨迹是围绕中心的闭合曲线,相应于周期性的往复摆动,总能量满足0?E?mgl;分界线之外的相轨迹是
不闭合的,它代表另一种周期运动,即向一个方向的转动,相应总能量满足
E?mgl.
对应于振幅较大(即总能量较大)的往复摆动,谐频的影响亦不容忽视,它
的相轨迹虽然闭合,但已不是椭圆,与简谐振动的相图是不同的. 图59
15
图59为对数值解进行快速傅里叶变换(fft)得到的功率谱图(纵轴为对数坐标),清晰可见相应?的基频谱线,相应3?和5?的谐频谱线;隐约可见相应7?的谐频谱线.由于谐频振动的振幅随其级数的增加而迅速减少,所以更高倍的谐频谱线无法显示.
四、自激振动 范·德·波尔方程
????(x02?x2)x???02x?0 (3.4) 图60 x依据(3.4)式编写微分方程函数文件zjzdfun.m如图61;再编写解微分方程并画相图的主程序zjzd.m如图62,参量取值为??0.3,x0?1,?0?1;运行,计算结果如图63?65.
图61 图62
图63 图64 图65
§4 能导致混沌的倒摆受迫振动
一、运动微分方程的建立
实验装置如图66,图67为示意原理图图67中,R为可绕O轴转动的均匀铝制圆盘,m为与盘固定的配重,并用m处指针标志盘的位置,总体上看是一个m在O轴上方的摆,故称为倒摆.螺旋形弹簧提供弹性回复力;D为电磁铁,提供线性阻尼;可调速电机E带动轮B旋转,轮B又经连杆L驱动弹簧A端,提供驱动力.
图66 图67
如图67,建立直角坐标系
Oxyz,以
?描述系统位置.系统所受对
z轴的力矩有: 重力矩
?,螺旋形弹簧线性恢复力矩?c?rmgsin??asin?,阻力矩?b?矩M,由于卷弹簧
A端被驱动而受周期性驱动力
cos?t.无驱动时系统有三个平衡位置,??0为不稳定平衡位置,其两侧还有O1和O2两个稳定平衡位置.
系统对z轴的转动惯量为I,由对z轴的角动量定理,得
???asin??b???c??Mcos?t I?
16