针对 谭述君, 高强, 钟万勰. 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 解有效数字均一致.
来源: https://www.cnblogs.com/hwzhou/p/12301902.html