引子^[1]

單擺，這個(gè)在中學(xué)物理都學(xué)過(guò)的東西，應(yīng)該是非常熟悉了。

圖片來(lái)源-維基百科

小角度簡(jiǎn)單擺

若最高處(

)的繩子和最低處(速度最大值)的繩子的角度為 ，則可使用下列公式算出它的振動(dòng)周期。

公式證明

擺球受力分析

繩與對(duì)稱線夾角為

，繩長(zhǎng)為 ，繩距離對(duì)稱線的水平距離為 。對(duì)于處在某一點(diǎn)的小球進(jìn)行受力分析，小球受到的力只有重力和拉力。

，這個(gè)力指向平衡位置，稱為回復(fù)力。

由于

很小，所以有近似 ，所以有下列成立

，令

下面推導(dǎo)彈簧的周期公式(如果知道直接跳過(guò))

彈簧簡(jiǎn)諧振動(dòng)

設(shè)彈簧的彈性系數(shù)為

，彈簧距離平衡點(diǎn)的距離為 ，對(duì)于簡(jiǎn)諧振動(dòng)的彈簧，有

，將其與上式聯(lián)立，有

，根據(jù)三角函數(shù)的周期公式，有

回到單擺的話題上來(lái)

剛剛推導(dǎo)的回復(fù)力

，代入公式彈簧的周期公式，得到

這是中學(xué)時(shí)候的單擺公式的推導(dǎo)，看得出，主要是由于條件

這個(gè)條件限制產(chǎn)生的 所導(dǎo)致的，下面來(lái)看看一般情況下的單擺，也就是任意條件的 。

一般情況下的單擺

設(shè)與對(duì)稱線的夾角為

，繩子的長(zhǎng)度為 ，球的角加速度為 ，忽略空氣阻力由受力分析，根據(jù)牛頓第二運(yùn)動(dòng)定理，有 。根據(jù)角加速度和(線)加速度之間的關(guān)系 ^[2]，有下列微分方程成立

（右邊取負(fù)號(hào)是因?yàn)榛貜?fù)力的方向）

圖片來(lái)源-百度百科

如果

,根據(jù)Taylor展開式， ，，所以對(duì)于充分小的 ，有 ，使用Taylor展開來(lái)理解和使用極限來(lái)理解是一樣的。

求解特殊情況，微分方程

，這是一個(gè)二階常系數(shù)齊次微分方程，

1.猜測(cè)解法

移項(xiàng)

，從這個(gè)等式可以看出，函數(shù) 經(jīng)過(guò)兩次微分之后，除了除了系數(shù)及其前面的±號(hào)改變了，本身并沒(méi)有發(fā)生其他改變，對(duì)于微分之后本身還不變的函數(shù)，很容易想到余弦函數(shù)，令 ，其中 為常數(shù)。

對(duì)

，求二階導(dǎo)，代入微分方程， ，即

周期為

，跟前面一致。

2.公式解法

微分方程

，形如 ，存在通用解法。

特征方程為

,這是一對(duì)共軛復(fù)根：對(duì)復(fù)數(shù) ， 。根為 , 通解為 ，即 ，

令

， ，由于 ，那么肯定有 ，可以化為

根據(jù)上面的結(jié)果，也可輕松知道周期為

求解一般情況，微分方程

能量守恒視角看待問(wèn)題

這個(gè)方程求解比較麻煩，下面換個(gè)思路，通過(guò)能量守恒的方式來(lái)看這個(gè)問(wèn)題。該系統(tǒng)的自由度為

，使用廣義坐標(biāo) 來(lái)描述^[3]。設(shè)系統(tǒng)動(dòng)能為 ，系統(tǒng)勢(shì)能為 ，小球質(zhì)量為 ，繩長(zhǎng) ，不計(jì)繩質(zhì)量。取擺點(diǎn)最低處為零勢(shì)能點(diǎn)。

從上式計(jì)算角速度

從這個(gè)式子里面可以看到，

都是常數(shù)，是因變量角速度 與自變量角度 的函數(shù)關(guān)系。下面來(lái)討論一下這個(gè)式子，令

得到式子為

角度和角速度之間的關(guān)系,圖標(biāo)對(duì)應(yīng)的是不同的e的值

上圖研究了角度與角速度之間的關(guān)系，這樣的圖叫做相平面圖。從上面圖可以看出

由于能量不足，擺僅做平衡位置附近的周期運(yùn)動(dòng)。

是一種臨界狀態(tài)，相當(dāng)于擺錘擺到最高位置的過(guò)程。

能量過(guò)大，而使擺角的絕對(duì)值隨時(shí)間之增加而無(wú)限增加，對(duì)應(yīng)于擺繞支點(diǎn)無(wú)窮次旋轉(zhuǎn)的運(yùn)動(dòng)過(guò)程。

微分方程數(shù)值解法

• 圖1
• 圖2
• 圖3
• 圖4

微分方程計(jì)算的相平面圖

求解的時(shí)域圖

在求解上面的微分方程的時(shí)候，使用的方法是

公式 。具體求解過(guò)程可以參考MATLAB的幫助文檔和參考書^[4]。

帶阻尼的單擺

下面來(lái)看看帶阻尼的單擺是怎么使用的其實(shí)在求解上面的代碼里面就已經(jīng)加入空氣阻力一項(xiàng)了，只不過(guò)將其值設(shè)為了0，下面來(lái)看看改變這個(gè)值會(huì)發(fā)生什么。一般來(lái)說(shuō)，空氣阻力公式為

這里設(shè)

， 為阻尼因子。

那么上面的微分方程變?yōu)?

跟上面一樣將，使用

求解。

時(shí)域圖

時(shí)間T=50s，阻力系數(shù)μ=0.1

時(shí)間T=50s，阻力系數(shù)μ=0.5

時(shí)間T=50s，阻力系數(shù)μ=0.8

相平面圖

時(shí)間T=500s，阻力系數(shù)μ=0.1

時(shí)間T=50s，阻力系數(shù)μ=0.1

時(shí)間T=50s，阻力系數(shù)μ=0.5

時(shí)間T=50s，阻力系數(shù)μ=0.8

從上面這些圖來(lái)看，加入空氣阻力之后確實(shí)是一種帶阻尼的震動(dòng)圖像的樣子。并且阻尼越大，能量耗散的也越快。

---

最后，其實(shí)要說(shuō)的是，強(qiáng)烈推薦這個(gè)視頻。其實(shí)前面所有討論的東西都在下面的幾張圖里面了，可以回味回味。

代碼

隱函數(shù)畫圖代碼

for i=0:0.2:2 e = 1+i; f = @(theta,theta_bar) theta_bar^2-cos(theta)-e+1; fimplicit(f);legend(['e=',num2str(e)]);hold on; end xlabel('角度theta');ylabel('角速度omega');title('相平面圖'); labels=num2cell(1:0.2:3); labels = cellfun(@num2str,labels,'UniformOutput',false); legend(labels);plotset; print('推導(dǎo)函數(shù)相平面圖','-dpng');

微分方程數(shù)值解代碼

subplot_er(2,2,1) [t,y] = ode45(@solve,[0,50],[pi/4,0]); plot(y(:,1),y(:,2)); xlabel('theta');ylabel('omega');plotset;axis equal; legend('theta_0=pi/4,omega_0=0'); subplot_er(2,2,2) [t,y] = ode45(@solve,[0,50],[pi/2,0]); plot(y(:,1),y(:,2)); xlabel('theta');ylabel('omega');plotset;axis equal; legend('theta_0=pi/2,omega_0=0'); subplot_er(2,2,3) [t,y] = ode45(@solve,[0,50],[pi,1]); plot(y(:,1),y(:,2)); xlabel('theta');ylabel('omega');plotset;axis equal; legend('theta_0=pi/2,omega_0=1'); subplot_er(2,2,4) [t,y] = ode45(@solve,[0,50],[2*pi,2]); plot(y(:,1),y(:,2)); xlabel('theta');ylabel('omega');plotset;axis equal; legend('theta_0=2pi,omega_0=2'); suptitle('時(shí)間T=50s'); function dydt=solve(t,y) g = 9.8;l=2*g;mu=0; dydt = [y(2);-mu*y(2)-g/l*sin(y(1))]; end

求解的時(shí)域圖

subplot_er(2,2,1) [t,y] = ode45(@solve,[0,50],[pi/4,0]); plot(t,y(:,1),'-o',t,y(:,2),'-o'); title('單擺數(shù)值解法');xlabel('Time t');ylabel('Solution y'); legend('theta=pi/4','omega=0'); subplot_er(2,2,2) [t,y] = ode45(@solve,[0,50],[pi/2,0]); plot(t,y(:,1),'-o',t,y(:,2),'-o'); title('單擺數(shù)值解法');xlabel('Time t');ylabel('Solution y'); legend('theta=pi/2','omega=0'); subplot_er(2,2,3) [t,y] = ode45(@solve,[0,50],[pi/2,1]); plot(t,y(:,1),'-o',t,y(:,2),'-o'); title('單擺數(shù)值解法');xlabel('Time t');ylabel('Solution y'); legend('theta=pi/2','omega=1'); subplot_er(2,2,4) [t,y] = ode45(@solve,[0,50],[pi/2,2]); plot(t,y(:,1),'-o',t,y(:,2),'-o'); title('單擺數(shù)值解法');xlabel('Time t');ylabel('Solution y'); legend('theta=pi/2','omega=2'); function dydt=solve(t,y) g = 9.8;l=2*g;mu=0; dydt = [y(2);-mu*y(2)-g/l*sin(y(1))]; end

參考

^單擺-維基百科?https://zh.wikipedia.org/wiki/%E6%93%BA

^理論力學(xué)-哈工大?https://book.douban.com/subject/3866935/

^力學(xué)-朗道?https://book.douban.com/subject/2059252/

^數(shù)值方法?https://book.douban.com/subject/4780614/

總結(jié)

以上是生活随笔為你收集整理的偏微分方程数值解法pdf_单摆-微分方程浅谈的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。