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