非線性單擺橢圓積分解

針對 譚述君, 高強, 鍾萬勰. Duhamel項的精細積分方法在非線性微分方程數值求解中的應用[J]. 計算力學學報, 2010(05):13-19.中的非線性單擺算例，查閱文獻獲得了其橢圓積分理論解的具體形式，具體推導過程詳見Belendez A, Pascual C, Mendez D I, et al. Exact solution for the nonlinear pendulum[J]. Revista Brasileira De Ensino De Fisica, 2007, 29(4): 645-648.。

非線性單擺控制方程

基本變量爲幅角\(\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求解

在Scipy.special中給出了\(\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.0223301784 0.9484651763 -0.8279887606 0.6653140605 -0.4673292729]，與 譚述君, 高強, 鍾萬勰. Duhamel項的精細積分方法在非線性微分方程數值求解中的應用[J]. 計算力學學報, 2010(05):13-19.中表1的Reference解有效數字均一致。