平面曲线曲率的计算-MATLAB

定义

• 曲率：曲线上某个点的切线方向角对弧长的变化率，表明曲线在该点的弯曲程度。
• 设有光滑曲线 $C$（即曲线 $C$ 二阶可导），在曲线 $C$ 上选定一点 $M$ 作为度量的基点，点 $M$ 对应于弧长 $s$ 且切线方向角为 $\alpha$，曲线 $C$ 上另一点 $M^{\prime }$ 对应于弧 $s+\Delta s$ 且切线方向角为 $\alpha +\Delta \alpha$。那么，弧 $\stackrel{⌢}{M{M}^{\prime }}$ 的长度为 $|\Delta s|$，从点 $M$ 到 $M^{\prime }$ 切线转过的角度为 $|\Delta \alpha |$。
• 平均曲率：采用单位弧长上切线转过的角度来表达弧 $\stackrel{⌢}{M{M}^{\prime }}$ 的平均弯曲程度，称为弧 $\stackrel{⌢}{M{M}^{\prime }}$ 的平均曲率，即：$$\overline{K}=\frac{|\Delta \alpha |}{|\Delta s|}$$
• 曲率：当 $\Delta s\to 0$ 时，即 $M^{\prime }\to M$ 时，平均曲率存在极限，称为曲线 $C$ 在点 $M$ 处的曲率，即：$$K=\underset{\Delta s\to 0}{\lim }\frac{|\Delta \alpha |}{|\Delta s|}=|\underset{\Delta s\to 0}{\lim }\frac{\Delta \alpha }{\Delta s}|$$

定义公式

设光滑曲线 $C$ 表示为 $y=f(x)$，则根据上述定义，有曲率计算公式：$$K=\frac{\ddot{y}}{(1+{\dot{y}}^2{)}^{\frac{3}{2}}}$$ 推导过程请参考相关教材或文献。

差分公式

然而，对于参数函数 $y=y(t),x=x(t)$ 或者离散数据，无法直接采用上述定义公式，上述定义公式经推导得到差分计算公式：$$K=\frac{|\dot{x}\ddot{y}-\ddot{x}\dot{y}|}{[{\dot{x}}^2+{\dot{y}}^2{]}^{\frac{3}{2}}}$$

MATLAB程序

%% 正弦函数
x = 0:0.01:2*pi;
y = sin(x);
figure(1);
plot(x,y);

%% 定义公式
j = abs(sin(x)) ./ ((1+cos(x).^2).^(3/2));
figure(2);
plot(x,j);

%% 差分公式
x1 = diff(x);	% 一阶导
x2 = diff(x1);	% 二阶导
y1 = diff(y);
y2 = diff(y1);
x2(length(x1)) = x2(end);	% 使数组维度一致
y2(length(y1)) = y2(end);
k = abs(x1.*y2-x2.*y1) ./ (x1.^2+y1.^2).^(3/2);
k(length(x)) = k(end);
figure(3);
plot(x,k);

运行结果如下：

图1 正弦曲线

图2 定义公式计算得曲率

图3 差分公式计算得曲率