【编程实践】使用matlab求椭圆弧长

1 椭圆的表达式

Standard equation:
$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$$
Parametric equation:
$$\begin{array}{r}x=acos\alpha \\ y=bsin\alpha \end{array}$$
General form equation:
$$A{x}^2+Bxy+C{y}^2+Dx+Ey+F=0$$
thereinto

$$\begin{array}{r}A=a^{2}si{n}^2\theta +b^{2}co{s}^2\theta \\ B=2(b^2-a^2)sin\theta cos\theta \\ C=a^{2}co{s}^2\theta +b^{2}si{n}^2\theta \\ D=-2A{x}_0-B{y}_0\\ E=-B{x}_0-2C{y}_0\\ F=-\frac{1}{2}(D{x}_0+E{y}_0)-a^2{b}^2\end{array}$$
$x_0,y_0$为椭圆圆心坐标，$a,b$分别为椭圆的长短半轴，$\theta$为长半轴与x轴的夹角。

2 椭圆的弧长计算理论公式

$$l=a{\int }_{{\theta }_1}^{{\theta }_2}\sqrt{1-e^{2}si{n}^2\theta }，{\theta }_1,{\theta }_2\in [0,\frac{\pi }{2}]$$
thereinto, $e=\sqrt{1-\frac{b^2}{a^2}}$. Here, the angle we talk about is in the first quadrant.

$$C=4l$$
Using the Taylor series expansion, the circumference of the ellipse is
$$C\approx \pi (a+b)$$

3 以具体案例计算

$$\frac{x^2}{3^2}+\frac{y^2}{2^2}=1$$

// 参数方程
theta = 0:0.01:2*pi %起止角度和步长
X = 3*cos(theta)
Y = 2*sin(theta)
plot(X,Y)

绘图输出显示🦢：

// 标准形式
theta = 0:0.01:2*pi
f =@(X,Y)X.^2/9+Y.^2/4-1
plot(X,Y)

绘图输出显示🦏：

// An highlighted block

% 参数设定
a = 3
b = 2
e = sqrt(1-b^2/a^2)

theta_1 = 0 %起始角度
theta_2 = pi/2 %终止角度
theta_x = theta_1:0.01:theta_2 % 被积区间


y= a * sqrt(1-e.^2*sin(theta_x).^2) % 直接以被积函数声明
syms theta_vx
y2=@(theta_vx)a * sqrt(1-e.^2*sin(theta_vx).^2) % 匿名函数，被积函数句柄


trapz(theta_x,y) % 梯度积分，trapz(X,Y)根据X指定的坐标或标量间距对Y进行积分，因椭圆第一象限曲线为凸曲线，则梯度法计算结果略小于精确值
quadgk(y2,theta_1,theta_2) % 高斯-勒让德积分，使用高阶全局自适应积分和默认误差容限在区间内对函数句柄求积分
integral(y2,theta_1,theta_2) % 数值积分，使用全局自适应积分和模型误差容限在区间内以数值形式为函数求积分

运行结果🐳：

另，$[0,\frac{\pi }{2}]$区间四分之一椭圆弧长不精确计算结果：
$$ans=\pi \ast (a+b)/4=3.1415\ast (2+3)/4=3.9268$$

相同角度范围，椭圆弧长结果对比：
----------🌴----------$[\frac{\pi }{4},\frac{\pi }{2}]$间距下，

----------🌴----------$[0,\frac{\pi }{4}]$间距下，