欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

MATLAB 如何求取离散点的曲率最大值

程序员文章站 2022-03-25 12:49:20
我就废话不多说了,大家还是直接看代码吧~x0 = linspace(0.1,2,100);%x0,y0验证函数离散点,可以非等间隔y0 = 1./x0;h1 = abs(diff(x0)) ;h =...

我就废话不多说了,大家还是直接看代码吧~

x0 = linspace(0.1,2,100);%x0,y0验证函数离散点,可以非等间隔
y0 = 1./x0;
h1 = abs(diff(x0)) ;
h = [h1 h1(end)];
ht = h;
yapp1 = gradient(y0)./ht; %matlab数值近似
yapp2 = del2(y0)./ht; %matlab数值近似
k2 = abs(yapp2)./(1+yapp1.^2).^(3/2);
figure
plot(k2)
title('曲率曲线')
[~,maxflag] = max(k2);%曲率最大位置
x_max = x0(maxflag);
y_max = y0(maxflag);
%画出图像 标注曲率最大点
figure
plot(x0,y0,'.-');
hold on;
plot(x_max,y_max,'rp')
title('标注最大曲率点')
xlabel('log10((norm(b*xk-l)))')
ylabel('log10((norm(xk)))')

MATLAB 如何求取离散点的曲率最大值

补充:matlab 插值+计算离散点曲率

思路:点足够密的话直接用 diff、gradient 求曲率,稀疏的话先插值再算曲率。

公式:

MATLAB 如何求取离散点的曲率最大值

点密的情况 输入曲线坐标(1-2)求一、二阶导数(4-9)通过公式求得曲率(10)

x = 0:0.01:7;
y = cos(x*0.5*pi); 
h1 = abs(diff(x));
h = [h1 h1(end)];
ht = h; 
y1 = gradient(y)./ht;
y2 = gradient(y1)./ht;
curv = abs(y2)./sqrt((1+y1.^2).^3); 
plot(x,y,'-',x,curv,'--r);
legend('raw data, 'curvature','location',"best");
grid on

图像与下文理论值图像相同

点稀疏的情况

1、输入散点坐标(1-2)

2、用样条曲线(b-spline)等方法插值得到拟合曲线(3-4)

3、diff、gradient 函数求拟合曲线的一、二阶导数(6-11)

4、通过公式求得曲率(12)

例:余弦函数取 8 个点,用 b-spline 插值

x = 0:1:7;
y = cos(x*0.5*pi);
xx = 0:0.01:7;
yy = spline(x,y,xx); 
h1 = abs(diff(xx));
h = [h1 h1(end)];
ht = h; 
yy1 = gradient(yy)./ht;
yy2 = gradient(yy1)./ht;
curv = abs(yy2)./sqrt((1+yy1.^2).^3); 
plot(xx,yy,'-',xx,curv,'--r',x,y,'o-');
legend('b-spline', 'curvature','raw data','location',"best");
grid on

MATLAB 如何求取离散点的曲率最大值

补充用法

求最大曲率并在图中标出

[max_val,max_ind]=max(curv);
hold on
plot(xx(max_ind),yy(max_ind),'*r');

与理论值(余弦函数曲线)对比

曲线对比

MATLAB 如何求取离散点的曲率最大值

曲率对比

几种插值方法对比

列举四种方法,分别为:分段线性插值、三次样条曲线(b-spline)插值、三次 hermite 插值(pchip)、修正 akima 分段三次 hermite 插值(akima)

case 1: 三维螺线

MATLAB 如何求取离散点的曲率最大值

三维螺线散点

MATLAB 如何求取离散点的曲率最大值

插值

MATLAB 如何求取离散点的曲率最大值

俯视

MATLAB 如何求取离散点的曲率最大值

侧视

case 2:二维梯形波

MATLAB 如何求取离散点的曲率最大值

二维梯形波

case 3:三维不规则折线

MATLAB 如何求取离散点的曲率最大值

三维不规则折线(不等间距)

对比可得:

case 1:b-spline>akima>pchip>linear

case 2:linear>pchip>akima>b-spline

case 3:linear≈pchip≈akima>b-spline

故在插值的时候需要选择适合的计算方法

以上为个人经验,希望能给大家一个参考,也希望大家多多支持。如有错误或未考虑完全的地方,望不吝赐教。