引子^[1]

單擺，這個在中學物理都學過的東西，應該是非常熟悉了。

圖片來源-維基百科

小角度簡單擺

若最高處(

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

公式證明

擺球受力分析

繩與對稱線夾角為

，繩長為 ，繩距離對稱線的水平距離為 。對于處在某一點的小球進行受力分析，小球受到的力只有重力和拉力。

，這個力指向平衡位置，稱為回復力。

由于

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

，令

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

彈簧簡諧振動

設彈簧的彈性系數為

，彈簧距離平衡點的距離為 ，對于簡諧振動的彈簧，有

，將其與上式聯立，有

，根據三角函數的周期公式，有

回到單擺的話題上來

剛剛推導的回復力

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

這是中學時候的單擺公式的推導，看得出，主要是由于條件

這個條件限制產生的 所導致的，下面來看看一般情況下的單擺，也就是任意條件的 。

一般情況下的單擺

設與對稱線的夾角為

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

（右邊取負號是因為回復力的方向）

圖片來源-百度百科

如果

,根據Taylor展開式， ，，所以對于充分小的 ，有 ，使用Taylor展開來理解和使用極限來理解是一樣的。

求解特殊情況，微分方程

，這是一個二階常系數齊次微分方程，

1.猜測解法

移項

，從這個等式可以看出，函數 經過兩次微分之后，除了除了系數及其前面的±號改變了，本身并沒有發生其他改變，對于微分之后本身還不變的函數，很容易想到余弦函數，令 ，其中 為常數。

對

，求二階導，代入微分方程， ，即

周期為

，跟前面一致。

2.公式解法

微分方程

，形如 ，存在通用解法。

特征方程為

,這是一對共軛復根：對復數 ， 。根為 , 通解為 ，即 ，

令

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

根據上面的結果，也可輕松知道周期為

求解一般情況，微分方程

能量守恒視角看待問題

這個方程求解比較麻煩，下面換個思路，通過能量守恒的方式來看這個問題。該系統的自由度為

，使用廣義坐標 來描述^[3]。設系統動能為 ，系統勢能為 ，小球質量為 ，繩長 ，不計繩質量。取擺點最低處為零勢能點。

從上式計算角速度

從這個式子里面可以看到，

都是常數，是因變量角速度 與自變量角度 的函數關系。下面來討論一下這個式子，令

得到式子為

角度和角速度之間的關系,圖標對應的是不同的e的值

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

由于能量不足，擺僅做平衡位置附近的周期運動。

是一種臨界狀態，相當于擺錘擺到最高位置的過程。

能量過大，而使擺角的絕對值隨時間之增加而無限增加，對應于擺繞支點無窮次旋轉的運動過程。

微分方程數值解法

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

微分方程計算的相平面圖

求解的時域圖

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

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

帶阻尼的單擺

下面來看看帶阻尼的單擺是怎么使用的其實在求解上面的代碼里面就已經加入空氣阻力一項了，只不過將其值設為了0，下面來看看改變這個值會發生什么。一般來說，空氣阻力公式為

這里設

， 為阻尼因子。

那么上面的微分方程變為

跟上面一樣將，使用

求解。

時域圖

時間T=50s，阻力系數μ=0.1

時間T=50s，阻力系數μ=0.5

時間T=50s，阻力系數μ=0.8

相平面圖

時間T=500s，阻力系數μ=0.1

時間T=50s，阻力系數μ=0.1

時間T=50s，阻力系數μ=0.5

時間T=50s，阻力系數μ=0.8

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

---

最后，其實要說的是，強烈推薦這個視頻。其實前面所有討論的東西都在下面的幾張圖里面了，可以回味回味。

代碼

隱函數畫圖代碼

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('推導函數相平面圖','-dpng');

微分方程數值解代碼

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('時間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

求解的時域圖

subplot_er(2,2,1) [t,y] = ode45(@solve,[0,50],[pi/4,0]); plot(t,y(:,1),'-o',t,y(:,2),'-o'); title('單擺數值解法');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('單擺數值解法');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('單擺數值解法');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('單擺數值解法');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

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

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

^數值方法?https://book.douban.com/subject/4780614/

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。