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