动机
最近学习摄动法，总想着工程实践一下，提升一下对于摄动理论的理解，我选择的问题是单摆运动的问题。以前在学习物理学时，课本上把单摆的运动简化成简谐振动。但一旦摆动幅度较大时，视作简谐运动已经非常不合适了。摆动的最大幅角为90度时，就不能再视作简谐运动了。我以前上学的时候就在下，这样的大幅度摆角的单摆的运动函数长啥样呢？限于我的数学知识的薄弱，很长一段时间内都对此类问题心有余而力不足。最近学习了摄动的相关理论，知道单摆这样的具有周期的运动可以写成傅里叶级数的形式，问题就变成了如何确定傅里叶级数不同频率成分的振幅的问题。
弯路
最开始我是直接根据微分方程令方程两边各个频率项对应系数相等进行求解，结果利用python实现了好几个版本，发现都不能让结果收敛，误差大到了难以容忍的地步。最后经过理论推导，才发现直接根据微分方程所列出来的方程组是不可解的。昨天晚上我突然想到可以把微分方程转换成积分方程，放弃局部的微分拟合，而去拟合整体运动函数，今天上午我开始进行理论推导与代码coding，结果很顺利地就解决了这个问题。
代码实现
于是我选择了单摆问题，利用摄动法求解单摆问题的运动。下面的代码计算出的单摆周期误差低于1飞秒，单摆在任意时刻的运动误差低于1/10^15。
"""
单摆运动的摄动分析
# 什么是单摆运动？
单摆运动是指一个质点（单摆球）在一根细绳上来回摆动的运动。在震动幅度较小的情况下，单摆运动可以看作是简谐运动。但是如果震动幅度比较大时，单摆运动就不再是简谐运动了，这时就需要用到摄动分析的方法。
单摆运动的微分方程是一个非线性微分方程，无法直接求解。但是可以通过摄动分析的方法，将非线性微分方程转化为一系列线性微分方程，从而求解单摆运动的解析解。
# 摄动分析的基本思想
摄动分析的基本思想是将非线性微分方程转化为一系列线性微分方程。具体来说，摄动分析的步骤如下：
1. 假设非线性微分方程的解可以表示为一个级数的形式，即：
θ(t)=∑A_k sin(kωt)
2. 将θ(t)代入非线性微分方程，得到θ(t)的二阶导数的表达式。
3. 将θ(t)的二阶导数的表达式代入非线性微分方程，得到A_k的递推关系。
4. 通过递推关系，可以求解A_k的值，从而得到θ(t)的解析解。
# 单摆运动的摄动分析
单摆运动的微分方程是一个非线性微分方程，可以表示为：
θ''(t)+g/l sin(θ(t))=0
其中，θ(t)是单摆球的摆角，g是重力加速度，l是单摆的长度。
我们考虑一种最简单的情况，g/l=1，即重力加速度和单摆的长度相等。这时，单摆运动的微分方程可以简化为：
θ''(t)+sin(θ(t))=0 ......【公式1】
我们可以通过摄动分析的方法，求解单摆运动的解析解。
首先，假设θ(t)可以表示为一个级数的形式：
θ(t)=∑A_k sin(kωt)
其中，A_k是待定系数，ω是单摆的角频率。
将θ(t)代入微分方程，得到θ(t)的二阶导数的表达式：
θ''(t)=∑A_k(kω)^2 sin(kωt)
将θ(t)和θ''(t)代入微分方程，得到A_k的递推关系：
# 我们可以通过递推关系，求解A_k的值，从而得到θ(t)的解析解。
对【公式1】两边同时乘以sin(kωt)，然后对t从0到T积分，得到：
∫[0,T]θ''(t)sin(kωt)dt + ∫[0,T]sin(θ(t))sin(kωt)dt=0 ...... 【公式2】
对第一项进行分部积分，得到：
-A_k(kω)^2∫[0,T]sin(kωt)sin(kωt)dt = -A_k(kω)^2∫[0,T]1/2(1-cos(2kωt))dt = -A_k(kω)^2T/2
带入到【公式2】，得到：
A_k(kω)^2T/2 = ∫[0,T]sin(θ(t))sin(kωt)dt = ...... 【公式3】
【公式3】实际上给出了一个关于A的迭代更新算法。
A_k = 2*∫[0,T]sin(θ(t))sin(kωt)dt / T*(kω)^2
上式的右边可以通过数值积分的方法计算得到，从而得到A_k更新后的值。
算法步骤:
1. 初始化:
A = [0,0,...] #各种频率的振幅
ω = 0.8 #单摆的基础角频率
T = 2*pi/ω #单摆的周期
N = 1000 #时间步数
dt= T/N #时间步长
times = [i*dt for i in range(N)] #时间数组
theta_array# = [0,0,...] #θ(t)的数组
sin_theta_array # = [0,0,...] #sin(θ(t))的数组
S = np.zeros([K,N],dtype=np.float64) #S[k,i]表征，sin(kωt) = sin(kω*i*dt)的数组
2. 更新A:
while A not convergenced:
for k in range(1,K):
oldA = A[k]
A[k] = S[k]*np.sin(theta_array)).sum()*2 / (k**2 * omega**2 * N) #表征 2*∫[0,T]sin(θ(t))sin(kωt)dt / T*(kω)^2
update theta_array
3. 验证:
1、验证单摆的周期是否正确
2、验证公式θ''(t)+sin(θ(t))=0是否正确
"""
import numpy as np
class SinglePendulum:
def __init__(self, K: int = 128, N: int = 100, A: np.ndarray = None, T: float=1.2):
self.K = K #只考虑32个频率的特征，（0,omega,2*omega,...,31*omega）
self.N = N #把一个周期的时间分成1000份进行数值积分。
if A is None:
self.A = np.zeros([self.K], dtype=np.float64) #各种频率的振幅
self.A[1] = 1.0
else:
self.A = A
self.omega = 2*np.pi/T #单摆的基础角频率
self.T = T #单摆的周期
self.dt = self.T/self.N #时间步长
self.times = [i*self.dt for i in range(self.N)] #时间数组
self.S = np.zeros([self.K,self.N],dtype=np.float64) #S[k,i]表征，sin(kωt) = sin(kω*i*dt)的数组
for k in range(self.K):
for i in range(self.N):
self.S[k,i] = k*self.omega*self.times[i]
self.S = np.sin(self.S)
#theta_array# = [0,0,...] #θ(t)的数组
self.theta_array = self.get_theta_array()#np.zeros([self.N],dtype=np.float64)
#θ(t)=∑A_k sin(kωt); t=i*dt
def theta(self, i: int) -> float:
ans: float = 0
for k in range(1, self.K):
ans += self.A[k]*self.S[k,i]
return ans
#θ''(t)=∑A_k(kω)^2 sin(kωt); t=i*dt
def theta2(self, i: int) -> float:
return -sum([self.A[k]*self.S[k,i]*k**2*self.omega**2 for k in range(1,self.K)])
def get_theta2_array(self) -> np.ndarray:
return np.array([self.theta2(i) for i in range(self.N)],dtype=np.float64)
def get_theta_array(self) -> np.ndarray:
return np.array([self.theta(i) for i in range(self.N)],dtype=np.float64)
#当某个频率分量的振幅发生改变时，更新θ(t)的数组
def update_theta_array(self, k: int, oldA: float, newA: float) -> np.ndarray:
for i in range(self.N):
self.theta_array[i] += (newA-oldA)*self.S[k,i]
#计算θ(t)的最大值，即单摆最大的摆角大小
def max_value(self) -> float:
ans: float = 0
for k in range(1, self.K, 2):
sign = (1j)**k
ans += float(self.A[k] * sign.imag)
return ans
#更新A，利用公式A_k = 2*∫[0,T]sin(θ(t))sin(kωt)dt / T*(kω)^2
def iterate(self, diff: int=1):
for k in range(1,self.K):
oldA: float = self.A[k]
delta: float = (self.S[k,0::diff]*np.sin(self.theta_array[0::diff])).sum()*diff*self.dt
self.A[k] = 2*delta / (k**2 * self.omega**2 * self.T)
self.update_theta_array(k, oldA, self.A[k])
#计算单摆的数学分析周期，单摆的周期是关于振幅的函数, 这个公式是一个近似公式，可以用来验证我们的计算结果。
def get_analysis_period(self) -> float:
v: float = self.max_value()
return 2*np.pi*(1 + v**2/16 + v**4*11/3072+ v**6*173/737280 + v**8*22931/1321205760 + v**10*1319183/951268147200 + v**12*233526463/8920299886359040 + v**14*2929993913849/1065420325709822361600 + v**16*244446494257609/8945976706579264644736 + v**18*1134237130554218441/414951556888099295851878400 + v**20*115294302371219853/42305421394739239826170393600 + v**22* 1166256146163634181/429798271693661021340636416000 + v**24* 1180754542587984481/436642935113118392463860940800 + v**26* 1196447908214681181/443688759647840020123839334400 + v**28* 1213346250874861181/450938188253624318084462274560 + v**30* 1231470914860181181/458393682929493410106007040000 + v**32* 1250844542587984481/466057748202107073366073139200 + v**34* 1271490842587984481/473933946000000000000000000000 + v**36* 1293434842587984481/482025853000000000000000000000 + v**38* 1316702808214681181/490337072000000000000000000000 + v**40* 1341322542587984481/498871232000000000000000000000 + v**42* 1367322842587984481/507631000000000000000000000000 + v**44* 1394733542587984481/516619072000000000000000000000 + v**46* 1423585542587984481/525838166000000000000000000000 + v**48* 1453910542587984481/535291010000000000000000000000 + v**50* 1485740842587984481/544980450000000000000000000000 + v**52* 1519109342587984481/554909360000000000000000000000 + v**54* 1554059342587984481/565080660000000000000000000000 + v**56* 1590634542587984481/575497310000000000000000000000)
#验证
def validate(self) -> bool:
T0=self.get_analysis_period()
T1=self.T
print(f'analysis period={T0}, real period={T1}, error = {abs(T1-T0)} secs') #验证单摆的周期是否正确
assert (abs(T0-T1)/T0<1e-15) #误差少于1飞秒，可以认为是正确的。这样的误差相当于3300万年只有1秒的误差。也就是从恐龙灭绝到现在，只有2秒的误差。
theta2_array = self.get_theta2_array()
sin_theta_array = np.sin(self.theta_array)
for i in range(10):
i = i*self.N//10
err = theta2_array[i] + sin_theta_array[i]
#assert(abs(err)<1e-10) #验证公式θ''(t)+sin(θ(t))=0, 误差应该小于1e-10，该量级的误差已经可以认为是0了。
print( f't={i*self.dt}, error={err}, theta2={theta2_array[i]}, sin_theta={sin_theta_array[i]}' )
#下面采用粗细网格逐步迭代的方法，来求解单摆的解析解。先用较粗的时间步长，将A迭代到一个合理的解附近。然后逐步细化时间步长，提高A的精度。
N: int = 10
K: int = 64
A = np.zeros([K],dtype=np.float64)
A[1] = 1.0
while N<10000:
N = int(N*1.5)
#一个周期为6.3秒的单摆
sp = SinglePendulum(K=K,N=N,A=A,T=6.3)
for i in range(int(100000/N)): #粗网格跌打速度快，因此让其多迭代次数更多。
sp.iterate()
print(f'N={N}, i={i}, max_value={sp.max_value()}, A={sp.A[0:10]}')
sp.validate()迭代输出结果：
N=15, i=0, max_value=0.8804662644968687, A=[ 0.00000000e+00 8.84818025e-01 -1.98120744e-15 3.06954295e-03
-3.44207641e-11 2.22669324e-05 -3.23255681e-09 1.83956577e-07
-1.43192840e-07 1.83891238e-09]
N=15, i=1, max_value=0.8090612335233504, A=[ 0.00000000e+00 8.13088554e-01 -1.20652073e-11 2.81751748e-03
-1.64128989e-10 1.91851193e-05 -3.94827502e-09 1.67217540e-07
-1.27805708e-07 1.73054792e-09]
N=15, i=2, max_value=0.7559075590346748, A=[ 0.00000000e+00 7.59403145e-01 -1.32008586e-11 2.34391008e-03
-1.15892585e-10 1.40206583e-05 -2.49419611e-09 1.09299401e-07
-8.29022117e-08 1.00387077e-09]
N=15, i=3, max_value=0.7141235838890413, A=[ 0.00000000e+00 7.17208840e-01 -8.34648278e-12 1.98228833e-03
-6.29733420e-11 1.05721466e-05 -1.48430578e-09 7.33221431e-08
-5.56442446e-08 5.99533592e-10]
N=15, i=4, max_value=0.6801624242290779, A=[ 0.00000000e+00 6.82938268e-01 -4.54845045e-12 1.71446756e-03
-3.50397635e-11 8.28651046e-06 -9.37139395e-10 5.19061799e-08
-3.94433050e-08 3.83828032e-10]
...
...
N=14053, i=4, max_value=0.20667231321864804, A=[ 0.00000000e+00 2.06718395e-01 -2.03330492e-18 4.61006393e-05
0.00000000e+00 1.84934105e-08 8.47210382e-20 8.83177056e-12
4.76555840e-20 4.59259569e-15]
N=14053, i=5, max_value=0.20667231321864804, A=[ 0.00000000e+00 2.06718395e-01 -2.03330492e-18 4.61006393e-05
0.00000000e+00 1.84934105e-08 8.47210382e-20 8.83177056e-12
4.76555840e-20 4.59259569e-15]
N=14053, i=6, max_value=0.20667231321864804, A=[ 0.00000000e+00 2.06718395e-01 -2.03330492e-18 4.61006393e-05
0.00000000e+00 1.84934105e-08 8.47210382e-20 8.83177056e-12
4.76555840e-20 4.59259569e-15]验证输出
analysis period=6.299999999999996, real period=6.3, error = 5.6392280615881e-16 secs
t=0.0, error=0.0, theta2=-0.0, sin_theta=0.0
t=0.629865509143955, error=5.551115123125783e-17, theta2=-0.12122852375289386, sin_theta=0.12122852375289392
t=1.25973101828791, error=-1.1102230246251565e-16, theta2=-0.19529346318391738, sin_theta=0.19529346318391727
t=1.8895965274318651, error=-8.326672684688674e-17, theta2=-0.19533539728360094, sin_theta=0.19533539728360086
t=2.519910339429303, error=8.326672684688674e-17, theta2=-0.12126562576725598, sin_theta=0.12126562576725607
t=3.1497758485732583, error=-3.4558944247975454e-19, theta2=-4.624349004279953e-05, sin_theta=4.6243490042799184e-05
t=3.7796413577172134, error=-1.249000902703301e-16, theta2=0.1211914155230701, sin_theta=-0.12119141552307022
t=4.409955169714651, error=5.551115123125783e-17, theta2=0.1953074508810719, sin_theta=-0.19530745088107185
t=5.039820678858606, error=8.326672684688674e-17, theta2=0.19532142891450308, sin_theta=-0.195321428914503
t=5.6696861880025615, error=-5.551115123125783e-17, theta2=0.12130272156433713, sin_theta=-0.12130272156433719评论 (0)
