第一类椭圆积分及其可视化

文章目录数学表达式第一类完全椭圆积分可视化第一类不完全椭圆积分可视化数学表达式椭圆积分是一个非常著名的函数家族和椭圆关系最近的是第二类其形式为E ( ϕ , m ) ∫ 0 ϕ 1 − m sin ⁡ 2 t d t E(\phi, m)\int^\phi_0\sqrt{1-m\sin^2t}\mathrm dtE(ϕ,m)∫0ϕ​1−msin2t​dt表示椭圆的弧长其中的被积函数有时候被称作椭圆根式成了椭圆家族最明显的标识。后来勒让德对椭圆积分进行分类将性质最好的划作第一类。第一类不完全椭圆积分的表达式为K ( ϕ , m ) ∫ 0 ϕ d t 1 − m sin ⁡ 2 t K(\phi, m)\int^\phi_0\frac{\mathrm dt}{\sqrt{1-m\sin^2t}}K(ϕ,m)∫0ϕ​1−msin2t​dt​这个积分的起源其实和椭圆关系不大最早出现在双扭线弧长的计算中双扭线的表达为( x 2 y 2 ) 2 a 2 ( x 2 − y 2 ) (x^2y^2)^2a^2(x^2-y^2)(x2y2)2a2(x2−y2)其弧长为s a ∫ 0 x d t 1 − t 4 sa\int^x_0\frac{\mathrm dt}{\sqrt{1-t^4}}sa∫0x​1−t4​dt​做代换t 2 sin ⁡ θ t^2\sin\thetat2sinθ即可得到得到m mm为定值的第一类椭圆积分。第一类椭圆积分最直观的物理意义是大角度单摆的无量纲时间。设摆长为L LL最大摆角θ 0 \theta_0θ0​根据能量守恒定律可以推导出单摆从θ 0 \theta0θ0摆动到任意角度θ \thetaθ所需的时间微元d t L 2 g d θ cos ⁡ θ − cos ⁡ θ 0 \mathrm dt \sqrt{\frac{L}{2g}}\frac{\mathrm d\theta}{\sqrt{\cos\theta-\cos\theta_0}}dt2gL​​cosθ−cosθ0​​dθ​引入无量纲时间τ t g L \taut\sqrt{\frac{g}{L}}τtLg​​, 令sin ⁡ ϕ sin ⁡ θ / 2 sin ⁡ θ 0 / 2 \sin\phi\frac{\sin\theta/2}{\sin\theta_0/2}sinϕsinθ0​/2sinθ/2​则上式变为K ( ϕ , m ) ∫ 0 ϕ d θ 1 − m sin ⁡ 2 θ K(\phi, m) \int^\phi_0\frac{\mathrm d\theta}{\sqrt{1-m\sin^2\theta}}K(ϕ,m)∫0ϕ​1−msin2θ​dθ​其中m sin ⁡ 2 θ 0 2 m\sin^2\frac{\theta_0}{2}msin22θ0​​。当ϕ π 2 \phi\frac{\pi}{2}ϕ2π​时上式为第一类完全椭圆积分K ( m ) ∫ π 2 ϕ d θ 1 − m sin ⁡ 2 θ K(m) \int^\phi_\frac{\pi}{2}\frac{\mathrm d\theta}{\sqrt{1-m\sin^2\theta}}K(m)∫2π​ϕ​1−msin2θ​dθ​第一类完全椭圆积分可视化scipy.special中提供了第一类椭圆积分其中ellipk和ellipkm1是完全积分ellipkinc是不完全积分。考虑到完全椭圆积分是不完全椭圆积分的特例下面先演示ellipk和ellipkm1。这两个函数的区别是前者与上文K ( m ) K(m)K(m)的表达式相同后者主攻m 1 m1m1附近的K ( 1 − p ) K(1-p)K(1−p)。代码如下importnumpyasnpfromscipy.specialimportellipk,ellipkinc,ellipkm1importmatplotlib.pyplotasplt plt.rcParams[font.family][Times New Roman]xsnp.linspace(0,1,100)ysellipk(xs)plt.plot(xs,ellipk(xs),labelellipk)plt.plot(xs,ellipkm1(xs),labelellipkm1)plt.legend()plt.grid()plt.show()第一类不完全椭圆积分可视化第二类不完全椭圆积分可视化效果如下代码为phinp.linspace(0,np.pi,100)mnp.linspace(0,1,100)Phi,Mnp.meshgrid(phi,m)Eellipkinc(Phi,M)axplt.subplot(projection3d)ax.plot_surface(Phi,M,E,cmapjet,alpha0.9)ax.set_xlabel(r$\phi$)ax.set_ylabel(r$m$)ax.set_zlabel(r$E(\phi, m)$)plt.show()