[ 常微分方程 ] 04 高阶微分方程实例

高阶微分方程一般用于一些具体的物理情景中，下面以质点振动和宇宙速度的推导为例。参考书：王高雄《常微分方程（第四版）》

文章目录

• 一、质点振动
• • 01 无阻尼自由振动
  • • （1）物理推导
    • （2）微分方程推导
  • 02 有阻尼自由摆动
  • 03 无阻尼强迫振动
  • 04 有阻尼强迫振动
  • 05 质点振动小结

一、质点振动

01 无阻尼自由振动

（1）物理推导

当我们在研究一个数学摆的简谐运动的时候，分析其受力如下图所示，我们将重力分为切向和法向两个方向。

设绳子的拉力为 T T T，有在法向和切向的受力分析有如下表达式：
{ T − m g cos ⁡ φ = F n m g sin ⁡ φ = F t = − m a t ⏟ 牛顿第二定律 = − m d 2 l d t 2 = − m L d 2 φ d t 2 \begin{cases} T-mg\cos \varphi=F_n\\ \boxed{mg\sin\varphi}=\underbrace{F_t=-ma_t}_{牛顿第二定律}=-m\displaystyle\frac{d^2l}{dt^2}=\boxed{-mL\frac{d^2\varphi}{dt^2}} \end{cases} ⎩ ⎨ ⎧​T−mgcosφ=Fn​mgsinφ​=牛顿第二定律 Ft​=−mat​​​=−mdt2d2l​=−mLdt2d2φ​​​
对上面画框的部分进行整理可以得到
d 2 φ d t 2 = − g L sin ⁡ φ \frac{d^2\varphi}{dt^2}=-\frac{g}{L}\sin\varphi dt2d2φ​=−Lg​sinφ
怎么和书上的不一样？事实上，物理人为了简化模型，由于单摆的角度极其微小，此时认为 sin ⁡ φ ≈ φ \sin\varphi\approx \varphi sinφ≈φ，并令 ω 2 = g L \omega^2=\displaystyle\frac{g}{L} ω2=Lg​，得到书上的表达式
d 2 φ d t 2 + ω 2 φ = 0 (1) \frac{d^2\varphi}{dt^2}+\omega^2\varphi=0 \tag{1} dt2d2φ​+ω2φ=0(1)

（2）微分方程推导

这是一个二阶常系数微分方程，其特征方程是
λ 2 + ω 2 = 0 \lambda^2+\omega^2=0 λ2+ω2=0
解出共轭复根 λ 1 , 2 = ± ω i \lambda_{1,2}=\pm\omega i λ1,2​=±ωi，其通解为
φ = c 1 cos ⁡ ω t + c 2 sin ⁡ ω t = c 2 2 + c 2 2 ( c 1 c 1 2 + c 2 2 cos ⁡ ω t + c 2 c 1 2 + c 2 2 sin ⁡ ω t ) ⏟ 辅助角公式 = A ( sin ⁡ θ cos ⁡ ω t + cos ⁡ θ sin ⁡ ω t ) = A sin ⁡ ( ω t + θ ) (2) \begin{align*} \varphi&=c_1\cos\omega t+c_2\sin\omega t\\ &=\underbrace{\sqrt{c_2^2+c_2^2}\left(\frac{c_1}{\sqrt{c_1^2+c_2^2}}\cos\omega t+\frac{c_2}{\sqrt{c_1^2+c_2^2}}\sin\omega t\right)}_{辅助角公式}\\ &=A\left( \sin\theta\cos\omega t+\cos\theta\sin \omega t\right)\\ &=A\sin(\omega t+\theta) \end{align*}\tag{2} φ​=c1​cosωt+c2​sinωt=辅助角公式 c22​+c22​ ​(c12​+c22​ ​c1​​cosωt+c12​+c22​ ​c2​​sinωt)​​=A(sinθcosωt+cosθsinωt)=Asin(ωt+θ)​(2)
可以看出单摆的运动轨迹服从正弦函数，我们利用常微分方程的初值问题定义，给不同的初值赋予不同的物理条件，就可以得到不同的物理意义。
{ 周期 : T = 2 π ω 角频率： ω = 2 π ν 频率 : ν = 1 T = ω 2 π \begin{cases} \displaystyle周期:T=\frac{2\pi}{\omega}\\ 角频率：\omega=2\pi \nu\\ \displaystyle频率:\nu=\frac{1}{T}=\frac{\omega}{2\pi} \end{cases} ⎩ ⎨ ⎧​周期:T=ω2π​角频率：ω=2πν频率:ν=T1​=2πω​​
由于这是二阶微分方程，需要两个常数作为初值(主要就是指位置和一阶导速度），赋予 t = 0 t=0 t=0时的初值条件为
{ φ ∣ t = 0 = φ 0 ( 初始角度为 φ 0 ) d φ d t ∣ t = 0 = 0 ( 初始速度为 0 ) \begin{cases} \varphi\big|_{t=0}=\varphi_0\quad(初始角度为\varphi_0)\\\\ \displaystyle\frac{d\varphi}{dt}\bigg|_{t=0}=0\quad(初始速度为0) \end{cases} ⎩ ⎨ ⎧​φ ​t=0​=φ0​(初始角度为φ0​)dtdφ​ ​t=0​=0(初始速度为0)​
此时代入式 ( 2 ) (2) (2)中所得到的通解为
φ = φ 0 sin ⁡ ( ω t + π 2 ) \varphi=\varphi_0\sin\left(\omega t+\frac{\pi}{2}\right) φ=φ0​sin(ωt+2π​)
实际就是一个物理情景：我们将该单摆拉到 φ 0 \varphi_0 φ0​的位置并由静止状态松开，让其自由摆动。

02 有阻尼自由摆动

通常我们认为空气阻力和一个物体运动的速度成正比，因此我们在式 ( 1 ) (1) (1)中加入阻尼因子得到
d 2 φ d t 2 + μ m d φ d t ⏟ 表示阻尼 + g L φ = 0 (3) \frac{d^2\varphi}{dt^2}+\underbrace{\frac{\mu}{m}\frac{d\varphi}{dt}}_{表示阻尼}+\frac{g}{L}\varphi=0 \tag{3} dt2d2φ​+表示阻尼 mμ​dtdφ​​​+Lg​φ=0(3)
为了方程更加简洁，物理人会给常数一个更简单的符号，（其中 n n n指示阻尼的大小），替换得到
d 2 φ d t 2 + 2 n d φ d t + ω 2 φ = 0 \frac{d^2\varphi}{dt^2}+2n\frac{d\varphi}{dt}+\omega^2\varphi=0 dt2d2φ​+2ndtdφ​+ω2φ=0
其特征方程为
λ 2 + 2 n λ + ω 2 = 0 ⟹ λ 1 , 2 = − n ± n 2 − ω 2 \lambda^2+2n\lambda+\omega^2=0\\ \Longrightarrow \lambda_{1,2}=-n\pm\boxed{\sqrt{n^2-\omega^2}} λ2+2nλ+ω2=0⟹λ1,2​=−n±n2−ω2 ​​
可以看到这里面有根号，如果说被开根的数是负数的话，就要涉及复数了，这里我们进行分类讨论:

• 小阻尼： 当 n < ω n<\omega n<ω，被开根数小于0， λ 1 , 2 = − n + ω 2 − n 2 ⏟ ω 1 i \lambda_{1,2}=-n+\underbrace{\sqrt{\omega^2-n^2}}_{\displaystyle\omega_1}i λ1,2​=−n+ω1​ ω2−n2 ​​​i，此时通解为
  φ = e − n t ( c 1 cos ⁡ ω 1 t + c 2 sin ⁡ ω 1 t ) = A e − n t sin ⁡ ( ω 1 t + θ ) \begin{align*} \varphi&=e^{-nt}(c_1\cos\omega_1 t+c_2\sin\omega_1 t)\\ &=\boxed{Ae^{-nt}}\sin(\omega_1t+\theta) \end{align*} φ​=e−nt(c1​cosω1​t+c2​sinω1​t)=Ae−nt​sin(ω1​t+θ)​
  可以看到上面振幅和 t t t有关，并且是负相关，意味着有阻尼时的振动强度是逐渐减弱的，虽然已经没有了标准的周期性（指振幅不一致了），但是摆动的最大偏离的间隔时间仍然相同，也就是仍具有振动性质。
  
  
  

• 大阻尼： 当 n > ω n>\omega n>ω，此时 λ 2 < λ 1 < 0 \lambda_2<\lambda_1<0 λ2​<λ1​<0，其通解为
  φ = c 1 e λ 1 t + c 2 e λ 2 t \varphi=c_1 e^{\lambda_1t}+c_2e^{\lambda_2t} φ=c1​eλ1​t+c2​eλ2​t
  研究解的图像可以发现，由
  { φ = c 1 e λ 1 t + c 2 e λ 2 t = 0 ( 解曲线最多只有一个根 ) d φ d t = e λ 1 t [ c 1 λ 1 + c 2 λ 2 e ( λ 2 − λ 1 ) t ] ( t 很大时符号与 c 1 相反 ) \begin{cases} \varphi=c_1e^{\lambda_1t}+c_2e^{\lambda_2t}=0\quad(解曲线最多只有一个根)\\ \displaystyle\frac{d\varphi}{dt}=e^{\lambda_1t}[c_1\lambda_1 +c_2\lambda_2e^{(\lambda_2-\lambda_1)t}]\quad (t 很大时符号与c_1相反) \end{cases} ⎩ ⎨ ⎧​φ=c1​eλ1​t+c2​eλ2​t=0(解曲线最多只有一个根)dtdφ​=eλ1​t[c1​λ1​+c2​λ2​e(λ2​−λ1​)t](t很大时符号与c1​相反)​
  由 e e e指数可以看出，在 t t t很大时，其一阶导绝对值逐渐减小，也就是曲线逐渐趋于平稳，此时摆并不具备振动性质了。
  

• 临界阻尼： 也就是 n = ω n=\omega n=ω的临界状态（临界阻尼值），此时特征方程有重根 λ 1 = λ 2 = − n \lambda_1=\lambda_2=-n λ1​=λ2​=−n，此时阻尼恰好可以抑制振动，而呈现出如上图所示的效果，而当 n n n向下变化一个很小的范围时，摆就又具有振动性质了，这就又回到了小阻尼的图像。

03 无阻尼强迫振动

加入外力的因素，且不考虑阻尼，即方程 ( 1 ) (1) (1)变为
d 2 φ d t 2 + g L φ = 1 m L F ( t ) (4) \frac{d^2\varphi}{dt^2}+\frac{g}{L}\varphi=\frac{1}{mL}F(t) \tag{4} dt2d2φ​+Lg​φ=mL1​F(t)(4)
当我们施加的力以正弦函数为周期进行变化时，设 F ( t ) m L = H sin ⁡ p t \displaystyle\frac{F(t)}{mL}=H\sin pt mLF(t)​=Hsinpt，其对应齐次线性微分方程的通解为
φ = A sin ⁡ ( ω t + θ ) \varphi=A\sin(\omega t+\theta) φ=Asin(ωt+θ)

• 若 ω ≠ p \omega\neq p ω=p，式 ( 4 ) (4) (4)有形如下式的特解，并对照系数之后解得：
  φ ~ = M ⏟ 0 cos ⁡ p t + N ⏟ H ω 2 − p 2 sin ⁡ p t \widetilde{\varphi}=\underbrace{M}_{0}\cos pt +\underbrace{N}_{\displaystyle\frac{H}{\omega^2-p^2}}\sin pt φ ​=0 M​​cospt+ω2−p2H​ N​​sinpt
  则式 ( 4 ) (4) (4)的通解为
  φ = A sin ⁡ ( ω t + θ ) ⏟ 固有振动 + H ω 2 − p 2 sin ⁡ p t ⏟ 强迫振动 \varphi=\underbrace{A\sin(\omega t+\theta)}_{固有振动}+\underbrace{\frac{H}{\omega^2-p^2}\sin pt}_{强迫振动} φ=固有振动 Asin(ωt+θ)​​+强迫振动 ω2−p2H​sinpt​​
  可以发现，若外力的角频率 p p p越接近固有振动的角频率 ω \omega ω，则强迫振动的振幅就越大，那如果二者相等呢？
• 若 ω = p \omega=p ω=p，特别的，此时振幅一定是最大的。式 ( 4 ) (4) (4)有如下形式的特解：
  φ ~ = t ( M ⏟ − H 2 ω cos ⁡ ω t + N ⏟ 0 sin ⁡ ω t ) \widetilde{\varphi}=t(\underbrace{M}_{\displaystyle-\frac{H}{2\omega}}\cos\omega t+\underbrace{N}_{0}\sin\omega t) φ ​=t(−2ωH​ M​​cosωt+0 N​​sinωt)
  即方程的通解为
  φ = A sin ⁡ ( ω t + θ ) − H 2 ω t cos ⁡ ω t \varphi=A\sin(\omega t+\theta)-\frac{H}{2\omega}t\cos \omega t φ=Asin(ωt+θ)−2ωH​tcosωt
  该方程描述了，随着时间的增大，摆的偏离程度将无限增加，这称为共振现象，但这是不符合物理直觉和实际生活的，别忘了，我们在推导单摆运动方程的时候， sin ⁡ φ ≈ φ \sin\varphi\approx \varphi sinφ≈φ的成立条件是在角度极小的时候，也就是说当摆的偏离到达一定程度时，已经不能用该方程来描述单摆运动了。

04 有阻尼强迫振动

对于
d 2 φ d t 2 + 2 n d φ d t + ω 2 φ = H sin ⁡ p t (5) \frac{d^2\varphi}{dt^2}+2n\frac{d\varphi}{dt}+\omega^2\varphi=H\sin pt\tag{5} dt2d2φ​+2ndtdφ​+ω2φ=Hsinpt(5)
这里我们只讨论小阻尼情形，即 n < ω n<\omega n<ω，对应齐次线性微分方程通解及式 ( 5 ) (5) (5)的特解为
φ = A e − n t sin ⁡ ( ω 1 t + θ ) \varphi=Ae^{-nt}\sin(\omega_1 t+\theta) φ=Ae−ntsin(ω1​t+θ)
φ ~ = M ⏟ − 2 n p H ( ω 2 − p 2 ) 2 + 4 n 2 p 2 cos ⁡ p t + N ⏟ ( ω 2 − p 2 ) H ( ω 2 − p 2 ) 2 + 4 n 2 p 2 sin ⁡ p t \widetilde{\varphi}=\underbrace{M}_{\displaystyle\frac{-2npH}{(\omega^2-p^2)^2+4n^2p^2}}\cos pt+\underbrace{N}_{\displaystyle\frac{(\omega^2-p^2)H}{(\omega^2-p^2)^2+4n^2p^2}}\sin pt φ ​=(ω2−p2)2+4n2p2−2npH​ M​​cospt+(ω2−p2)2+4n2p2(ω2−p2)H​ N​​sinpt
由于观察得到
M 2 + N 2 = H 2 ( ω 2 − p 2 ) 2 + 4 n 2 p 2 M^2+N^2=\frac{H^2}{(\omega^2-p^2)^2+4n^2p^2} M2+N2=(ω2−p2)2+4n2p2H2​
则我们可以将其简化，得到
{ M = H ( ω 2 − p 2 ) 2 + 4 n 2 p 2 cos ⁡ θ ∗ = H ∗ cos ⁡ θ ∗ N = H ( ω 2 − p 2 ) 2 + 4 n 2 p 2 sin ⁡ θ ∗ = H ∗ sin ⁡ θ ∗ \begin{cases} M=\displaystyle\boxed{\frac{H}{\sqrt{(\omega^2-p^2)^2+4n^2p^2}}}\cos\theta^*=H^*\cos\theta^*\\ N=\displaystyle\boxed{\frac{H}{\sqrt{(\omega^2-p^2)^2+4n^2p^2}}}\sin\theta^*=H^*\sin\theta^* \end{cases} ⎩ ⎨ ⎧​M=(ω2−p2)2+4n2p2 ​H​​cosθ∗=H∗cosθ∗N=(ω2−p2)2+4n2p2 ​H​​sinθ∗=H∗sinθ∗​
则式 ( 5 ) (5) (5)的通解为
φ = A e − n t sin ⁡ ( ω 1 t + θ ) + H ∗ sin ⁡ ( p t + θ ∗ ) \varphi=Ae^{-nt}\sin(\omega_1 t+\theta)+H^*\sin (pt+\theta^*) φ=Ae−ntsin(ω1​t+θ)+H∗sin(pt+θ∗)
可以观察发现，振动仍由两部分组成，但第二部分虽然和外力的频率一样，可是振幅已经不同了。此时我们仍然想知道当我们施加的外力符合什么条件时，振幅最大，此时需要 H ∗ H^* H∗最大，即需要 ( ω 2 − p 2 ) 2 + 4 n 2 p 2 (\omega^2-p^2)^2+4n^2p^2 (ω2−p2)2+4n2p2最小。记
Φ ( p ) = ( ω 2 − p 2 ) 2 + 4 n 2 p 2 \varPhi(p)=(\omega^2-p^2)^2+4n^2p^2 Φ(p)=(ω2−p2)2+4n2p2
求一次导并令其为零得：
Φ ′ ( p ) = − 4 p ( ω 2 − p 2 ) + 8 n 2 p = 0 ⟹ p ∗ = ω 2 − 2 n 2 ( 需要保证被开根数大于零，即为小阻尼情形 ) \varPhi'(p)=-4p(\omega^2-p^2)+8n^2p=0\\ \Longrightarrow p^*=\sqrt{\omega^2-2n^2}\quad(需要保证被开根数大于零，即为小阻尼情形) Φ′(p)=−4p(ω2−p2)+8n2p=0⟹p∗=ω2−2n2 ​(需要保证被开根数大于零，即为小阻尼情形)
此时 Φ ′ ′ ( p ∗ ) = 8 p 2 > 0 , H max ⁡ ∗ = H 2 n ω 2 − n 2 \varPhi''(p^*)=8p^2>0,\displaystyle H^*_{\max}=\frac{H}{2n\sqrt{\omega^2-n^2}} Φ′′(p∗)=8p2>0,Hmax∗​=2nω2−n2 ​H​， p p p达到了共振频率，单摆产生共振现象。

05 质点振动小结