摄动理论求解非线性单摆振动问题

概念

扰动理论也称摄动理论，求解微分方程的一种方法，精确解=具有精确解的简单部分+不具有精确解的微扰部分，扰动理论的精确解通常会表示为一个微小参数的幂级数。
扰动分析分为：regular perturbation与singular perturbation，扰动将一个不可解问题作为可解问题的扰动，奇异摄动问题中不可解问题与可解问题的解有一个本质的区别。当$ϵ\to 0$，解$u(x,ϵ)$一致收敛converge uniformly，则为regular perturbation，反之为singular perturbation。

regular perturbation，$u=\Sigma {ϵ}^m{u}_m(t)$
singular perturbation，$u\approx \Sigma {ϵ}^m{u}_m(t)$

扰动理论 
http://www.math.unm.edu/~vageli/courses/Ma570/text.pdf
https://www.sjsu.edu/faculty/watkins/perturbregsing.htm
物理学咬文嚼字之五十七 简并

非线性单摆振动问题

方程：${\theta }^{\prime \prime }+w_0^{2}sin(\theta )=0,w_0^2=g/l$
条件：$\theta (0)=\alpha ,{\theta }^{\prime }(0)=0$
将问题尺度化，令$u=\theta /\alpha$
方程：$u^{\prime \prime }+w_0^2/\alpha sin(u\alpha )=0$
条件：$u(0)=1,u^{\prime }(0)=0$

[1]赵炳林. 应用叠代法求单摆周期公式[J]. 大学物理, 1984, 1(5).
[2]许玉兴. 奇异摄动理论导数展开法在常微分方程中的应用[J]. 大学数学, 1987(04):33-35+12.
[3]程荣龙,宫昊,许永红.非线性单摆周期的摄动解法[J].蚌埠学院学报,2017,6(06):39-41.
[4]叶慧群.单摆周期近似解法综述[J].浙江师范大学学报(自然科学版),2004(03):34-38.
[5]黄秀兰.单摆周期近似解的讨论[J].云南师范大学学报(自然科学版),1995(02):38-45.
[6]钱伟长. 奇异摄动理论及其在力学中的应用[M]. 科学出版社, 1981.
[7] Nayfeh A H . Introduction to Perturbation Techniques[J]. Journal of Mathematical Physics, 1962, 3(5):936-945.

迭代法[1]

$u^{\prime \prime }+w_0^2/\alpha sin(u\alpha )+w_0^{2}u-w_0^{2}u=0$
$u^{\prime \prime }+w_0^{2}u=w_0^2(u-1/\alpha sin(u\alpha ))$

将上式右端泰勒展开
$u-1/\alpha sin(u\alpha )=u-1/\alpha (u\alpha -(u\alpha {)}^3/3!+(u\alpha {)}^5/5!-(u\alpha {)}^7/7!+...)=u^3{\alpha }^2/6-u^5{\alpha }^4/60+u^7{\alpha }^6/5040-...$
即
$u-1/\alpha sin(u\alpha )=\Sigma [(-1{)}^{n+1}{\alpha }^{2n}{u}^{2n+1}/(2n+1)!]$
假设$\alpha$很小，$ϵ={\alpha }^2$
取$u-1/\alpha sin(u\alpha )=u^3ϵ/6$
则$u^{\prime \prime }+w_0^{2}u=w_0^2{u}^3ϵ/6,式1$

设解为$u=\Sigma {ϵ}^m{u}_m(t)$

非齐次方程的通解 = 齐次的通解 + 非齐次的特解

第0次迭代，将$u=u_0$带入式1
$u_0^{\prime \prime }+w_0^2{u}_0=0$
求齐次方程的通解，得$u_0$，将其带入式1的右端，得
$u^{\prime \prime }+w_0^{2}u=w_0^2{u}_0^3ϵ/6,式2$

第1次迭代，将$u=u_0+u_1$带入式2
$u_1^{\prime \prime }+w_0^2{u}_1=w_0^2{u}_0^3/6$
求非齐次方程的特解，得$u_1$，将$u_0+u_1$带入式1的右端，得
$u^{\prime \prime }+w_0^{2}u=w_0^2(u_0+u_1{)}^3ϵ/6,式3$

第2次迭代，将$u=u_0+u_1+u_2$带入式3
求非齐次方程的特解，得$u_2$

…

最终解为$u=\Sigma {ϵ}^m{u}_m(t)$，只取到第1次迭代，则$u=u_0+ϵu_1$
再将初始条件带入$u$得到最终解，但结果并不保证一定收敛。

下面是我的两次迭代，但sympy算出来的文章中的不一样

# -*- coding: utf-8 -*-
"""
@time: 2021-08-27 下午 09:41
@author: leslie lee
"""
import sympy as sp

# # -------------------------------ODE1
# # 自变量
# t = sp.symbols('t')
# # 常数
# w0 = sp.symbols('w0')
# # 因变量
# u0 = sp.Function('u0')
# # ODE
# eq = sp.Eq(u0(t).diff(t,2) + w0**2*u0(t), 0)
# # initial conditions
# ics = {u0(0): 1, u0(t).diff(t).subs(t,0): 0}
# # Solve
# print(sp.dsolve(eq, u0(t), ics = ics))
# # print(sp.dsolve(eq, u0(t)))

# -------------------------------ODE2
# 自变量
t = sp.symbols('t')
# 常数
w0 = sp.symbols('w0')
# 因变量
u1 = sp.Function('u1')
# ODE
eq = sp.Eq(u1(t).diff(t,2) + w0**2*u1(t), w0**2*sp.cos(w0*t)**3/6)
# initial conditions
ics = {u1(0): 0, u1(t).diff(t).subs(t,0): 0}
# Solve
print(sp.sympify(sp.dsolve(eq, u1(t), ics = ics)))

摄动法[2]

导数展开法[7]
方程的解$u(t)$，根据不同的时间尺度$T_0=t,T_1=ϵt,...,T_n={ϵ}^{n}t$可将解变为$\hat{u}(T_0,T_1,T_2,...,T_n)$，$ϵ$为小参数
$\frac{d}{dt}=\frac{\partial }{\partial T_0}+\frac{\partial }{\partial T_1}+\frac{\partial }{\partial T_2}+...$
$\frac{d^2}{d{t}^2}=\frac{{\partial }^2}{\partial T_0^2}+2ϵ\frac{{\partial }^2}{\partial T_0\partial T_1}+{ϵ}^2(2\frac{{\partial }^2}{\partial T_0{T}_2}+\frac{{\partial }^2}{\partial T_1^2})+...$

$u^{\prime \prime }+w_0^{2}u=w_0^2{u}^3ϵ/6$
将导数展开带入方程
$\frac{d^{2}u}{d{t}^2}=\frac{{\partial }^{2}u}{\partial T_0^2}+2ϵ\frac{{\partial }^{2}u}{\partial T_0\partial T_1}+{ϵ}^2(2\frac{{\partial }^{2}u}{\partial T_0{T}_2}+\frac{{\partial }^{2}u}{\partial T_1^2})$
$w_0^{2}u+\frac{{\partial }^{2}u}{\partial T_0^2}+2ϵ\frac{{\partial }^{2}u}{\partial T_0\partial T_1}+{ϵ}^2(2\frac{{\partial }^{2}u}{\partial T_0{T}_2}+\frac{{\partial }^{2}u}{\partial T_1^2})=w_0^2{u}^3ϵ/6,式1$
解的形式为$u=u_0+ϵu_1+{ϵ}^2{u}_2$
将解带回式1，并令$ϵ$同次幂的系数等于0，则可以得到三个常微分方程
${ϵ}^0$阶，关于$u_0$的常微分方程1
${ϵ}^1$阶，关于$u_0$与$u_1$的常微分方程2
${ϵ}^2$阶，关于$u_0$与$u_1$与$u_2$的常微分方程3
求解方程1的通解得到$u_0$，将$u_0$带入方程2求特解得到$u_1$，将$u_0$与$u_1$带入方程3求特解得到$u_2$
最终解为$u=u_0+ϵu_1+{ϵ}^2{u}_2$