基于Matlab70的數(shù)學(xué)擺求解

1、基于Matlab7.0的數(shù)學(xué)擺求解NORTHWESTERN POLYTECHNICAL UNIVERSITY基于Matlab7.0的數(shù)學(xué)擺求解 一 數(shù)學(xué)擺分析數(shù)學(xué)擺是系于一根長度為的線上而質(zhì)量為的質(zhì)點(diǎn)M，在重力作用下，它在垂直于地面的平面上沿圓周運(yùn)動(dòng)，如下圖所示。我們來確定擺的運(yùn)動(dòng)方程。設(shè)取逆時(shí)針運(yùn)動(dòng)的方向作為計(jì)算擺與鉛垂線所成的角的正方向。質(zhì)點(diǎn)M沿圓周的切向速度可以表示為 。作用于質(zhì)點(diǎn)M的重力將白拉回平衡位置A 。把重力分解為兩個(gè)分量和，第一個(gè)分量沿半徑OM方向，與線的拉力相抵消，它不會(huì)引起質(zhì)點(diǎn)M的速度的數(shù)值的改變。第二個(gè)分量沿著圓周的切線方向，它引起質(zhì)點(diǎn)M的速度的數(shù)值的改變。因?yàn)榭偸鞘官|(zhì)點(diǎn)

2、M向著平衡位置A的方向運(yùn)動(dòng)，即當(dāng)角為正時(shí)，向減小的方向運(yùn)動(dòng)；當(dāng)角為負(fù)時(shí)，向增大的方向運(yùn)動(dòng)，所以的數(shù)值等于。因此，擺的運(yùn)動(dòng)方程是 (1-1)即 (1-2)如果只研究擺的微小振動(dòng)，即當(dāng)比較小時(shí)的情況，我們可以取的近似值代入方程(2)。這樣，就得到微小振動(dòng)時(shí)擺的運(yùn)動(dòng)方程 (1-3)如果我們假設(shè)擺是在一個(gè)粘性的介質(zhì)中擺動(dòng)，那么，沿著擺的運(yùn)動(dòng)方向就存在一個(gè)與速度成比例的阻力。如果阻力系數(shù)是，則擺的運(yùn)動(dòng)方程變?yōu)槿绻刂鴶[的運(yùn)動(dòng)方向恒有一個(gè)外力作用于它，這時(shí)擺的運(yùn)動(dòng)稱為強(qiáng)迫微小振動(dòng)，其方程為當(dāng)要確定擺的某一個(gè)特定的運(yùn)動(dòng)時(shí)，我們應(yīng)該給出擺的初始狀態(tài)：當(dāng)時(shí)， ， ，這里代表擺的初始位置，代表擺的初始角速度的大小

3、。二 無阻尼自由振動(dòng)2.1 解析法如果只研究擺的微小振動(dòng)，即當(dāng)比較小時(shí)的情況，我們可以取的近似值代入方程(1-2)。這樣，就得到微小振動(dòng)時(shí)擺的運(yùn)動(dòng)方程 (2-1)記 ，這里 是常數(shù)，(3)變?yōu)?(2-2)這是二階常系數(shù)齊次線性微分方程，它的特征方程為 (2-3)特征根為共軛復(fù)根 (2-4)因此，方程(3-1)的通解為 (2-5)其中，為常數(shù)。為了獲得明顯的物理意義，令 (2-6)因此，若取 (2-7)則(3-4)可以寫成 (2-8)即 (2-9)這里代替了作為通解中所含的兩個(gè)任意常數(shù)。2.2 軟件計(jì)算首先把二階常微分方程改寫為一階常微分方程組： (2-10)然后，把該一階常微分方程組用一個(gè)M文

4、件形式的函數(shù)來表示，代碼設(shè)置如下：%ODE.m%常微分方程function d;最后，用函數(shù)ode(45)來解這個(gè)微分方程，并畫圖計(jì)算結(jié)果，代碼設(shè)置如下：%ode45.m%解微分方程;；由上述語句得到如下圖所示的示意圖。三 有阻尼自由振動(dòng)3.1 解析法如果我們假設(shè)擺是在一個(gè)粘性的介質(zhì)中擺動(dòng)，那么，沿著擺的運(yùn)動(dòng)方向就存在一個(gè)與速度成比例的阻力。如果阻力系數(shù)是，則擺的運(yùn)動(dòng)方程變?yōu)?(3-1)記，這里是正常數(shù)，(4)可以寫成 (3-2)它的特征方程為 (3-3)特征根為 (3-4)對(duì)于不同的阻尼值，微分方程有不同形式的解，它表示不同的運(yùn)動(dòng)形式，現(xiàn)分下面三種情況進(jìn)行討論： 小阻尼的情形：即的情形，這時(shí)

5、為一對(duì)共軛復(fù)根，記，則 。而方程(4-1)的通解為 (3-5)和前面無阻尼的情形一樣，可以把上述通解改寫成如下形式： (3-6)這里為任意常數(shù)。 大阻尼的情形：即的情形，這時(shí)，特征方程(4-2)有兩個(gè)不同的負(fù)實(shí)根。方程(4-1)的通解為 (3-7)這里為任意常數(shù)。 臨界阻尼的情形：即的情形，這時(shí)特征方程(4-2)有重根。方程(4-1)的通解為 (3-8)這里為任意常數(shù)。3.2 軟件計(jì)算首先把二階常微分方程改寫為一階常微分方程組： (3-9)然后，把該一階常微分方程組用一個(gè)M文件形式的函數(shù)來表示，代碼設(shè)置如下：%ODE.m%常微分方程function d;；;最后，用函數(shù)ode(45)來解這個(gè)微

6、分方程，并畫圖計(jì)算結(jié)果，代碼設(shè)置如下：%ode45.m%解微分方程;；由上述語句得到如下圖所示的示意圖。四 無阻尼強(qiáng)迫振動(dòng)4.1 解析法這時(shí)擺的運(yùn)動(dòng)稱為強(qiáng)迫微小振動(dòng)，其方程為 (4-1)令，設(shè)，為已知常數(shù)，為外力圓頻率。這時(shí)(5)變?yōu)?(4-2)方程(5-1)的對(duì)應(yīng)齊次線性微分方程的通解為 (4-3)這里是任意常數(shù)。現(xiàn)求(5-1)的一個(gè)特解。如果，則(5-1)有形如因此，方程(5-1)的通解為 (4-4)4.2軟件計(jì)算首先把二階常微分方程改寫為一階常微分方程組： (4-5)然后，把該一階常微分方程組用一個(gè)M文件形式的函數(shù)來表示，代碼設(shè)置如下：%ODE.m%常微分方程function d;；;最

7、后，用函數(shù)ode(45)來解這個(gè)微分方程，并畫圖計(jì)算結(jié)果，代碼設(shè)置如下：%ode45.m%解微分方程;；由上述語句得到如下圖所示的示意圖。五 有阻尼強(qiáng)迫振動(dòng)5.1 解析法這時(shí)擺的運(yùn)動(dòng)方程為 (5-1)根據(jù)實(shí)際的需要，我們只討論小阻尼的情形，即的情形。這時(shí)(6)對(duì)應(yīng)的齊次線性微分方程的通解為 (5-2)這里為任意常數(shù)，。現(xiàn)求(6)的一個(gè)特解，這時(shí)可以尋求形如 (5-3)的特解，這里是待定常數(shù)。將(6-2)代入(6)，比較同類項(xiàng)系數(shù)，得到為了獲得更明顯的物理意義，令即令 (5-4)及這時(shí)(6-2)可以寫成 (5-5)因此，(6)的通解為 (5-6)5.2 軟件計(jì)算首先把二階常微分方程改寫為一階常微分方程組： (5-7)然后，把該一階常微分方程組用一個(gè)M文件形式的函數(shù)來表