企业官网建设流程全解析

数学表达式

椭圆积分是一个非常著名的函数家族，和椭圆关系最近的是第二类，其形式为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^2+y^2)^2=a^2(x^2-y^2)(x2+y2)2=a2(x2−y2)，其弧长为s = a ∫ 0 x d t 1 − t 4 s=a\int^x_0\frac{\mathrm dt}{\sqrt{1-t^4}}s=a∫0x​1−t4​dt​，做代换t 2 = sin ⁡ θ t^2=\sin\thetat2=sinθ即可得到得到m mm为定值的第一类椭圆积分。

第一类椭圆积分最直观的物理意义，是大角度单摆的无量纲时间。设摆长为L LL，最大摆角θ 0 \theta_0θ0​，根据能量守恒定律，可以推导出单摆从θ = 0 \theta=0θ=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}}dt=2gL​​cosθ−cosθ0​​dθ​

引入无量纲时间τ = t g L \tau=t\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}m=sin22θ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 m=1m=1附近的K ( 1 − p ) K(1-p)K(1−p)。

代码如下

importnumpyasnpfromscipy.specialimportellipk,ellipkinc,ellipkm1importmatplotlib.pyplotasplt plt.rcParams['font.family']=['Times New Roman']xs=np.linspace(0,1,100)ys=ellipk(xs)plt.plot(xs,ellipk(xs),label="ellipk")plt.plot(xs,ellipkm1(xs),label="ellipkm1")plt.legend()plt.grid()plt.show()

第一类不完全椭圆积分可视化

第二类不完全椭圆积分可视化效果如下

代码为

phi=np.linspace(0,np.pi,100)m=np.linspace(0,1,100)Phi,M=np.meshgrid(phi,m)E=ellipkinc(Phi,M)ax=plt.subplot(projection='3d')ax.plot_surface(Phi,M,E,cmap='jet',alpha=0.9)ax.set_xlabel(r'$\phi$')ax.set_ylabel(r'$m$')ax.set_zlabel(r'$E(\phi, m)$')plt.show()