非线性单摆椭圆积分解

⾮线性单摆椭圆积分解
针对中的⾮线性单摆算例，查阅⽂献得到了其椭圆积分理论解的具体形式，具体推导过程详见。

⾮线性单摆控制⽅程
基本变量为幅⾓\theta，设重⼒加速度为g，单摆长度为l，则⽆阻尼的⾮线性单摆控制⽅程为：
\frac{\rm d^{2} \theta}{\rm d t^{2}}+\omega_{0}^{2} \sin \theta=0
其中，\omega_0=\sqrt{\dfrac{g}{l}}。

椭圆积分解
设初始条件为：
\theta(0)=\theta_{0} \quad\left(\frac{\rm d \theta}{\rm d t}\right)_{t=0}=0
在此初始条件下，单摆的运动范围为[-\theta_0,\theta_0]。

其解为：
\theta(t)=2 \arcsin \left\{\sin \frac{\theta_{0}}{2} \rm {sn}\left[K\left(\sin ^{2} \frac{\theta_{0}}{2}\right)-\omega_0 t|\sin^2\theta\right]\right\}
其中，\rm{sn}(u|m)为雅可⽐椭圆函数，K(m)=\int_{0}^{1} \frac{\rm d z}{\sqrt{\left(1-z^{2}\right)\left(1-m z^{2}\right)}}为第⼀类完全椭圆积分，Python求解
在中给出了\rm{sn}(u|m)和K(m)的求解函数。

#单摆精确解，参考⽂献见上。

Duhamel项的精细积分⽅法在⾮线性微分⽅程数值求解中的应⽤_谭述君内算例
def myWrite(arr,filePath):
import os
with open(filePath,'w') as f:
flag=False
for temp in arr:
if flag==True:
f.write("\n"+str(temp))
else:
flag=True
f.write(str(temp))
print("写⼊完成")
f.close()
from scipy import special as S
import numpy as np
from math import *
import matplotlib.pyplot as plt
t=np.linspace(0, 5, num=5/0.01+1) #Time interval
g=9.80665
l=1
theta0=1.0472 #Initial condition
w0=sqrt(g/l)
k=np.sin(theta0/2)**2
#函数K为第⼀类完全椭圆积分函数 complete elliptical integral of the first kind$K(k)=\int_{0}^{\pi / 2} \frac{\mathrm{d} \theta}{\sqrt{1-k^{2} \sin ^{2} \theta}}$，⽤ellipk求解。

K=S.ellipk(k)
T=4/w0*K#⾮线性单摆的周期
sn=S.ellipj(K-w0*t,k)[0]
#雅可⽐椭圆函数(Jacobi elliptic function)，ellipj: u[0]:sn() u[1]:cn u[2]:dn u[3]:phi
theta=2*np.arcsin(sqrt(k)*sn)
plt.plot(t,theta)
plt.show()
myWrite(theta,"theta.txt")
myWrite(t,"t.txt")
print(theta[100:600:100])# Output results at t=1,2,3,4,5 s。

给出的结果theta[100:600:100]=[-1.022******* 0.9484651763 -0.8279887606 0.6653140605 -0.4673292729]，与中表1的Reference解有效数字均⼀致。

Processing math: 0%。