x0 = linspace(0.1,2,100);%x0,y0驗(yàn)證函數(shù)離散點(diǎn)，可以非等間隔
y0 = 1./x0;
h1 = abs(diff(x0)) ;
h = [h1 h1(end)];
ht = h;
yapp1 = gradient(y0)./ht; %matlab數(shù)值近似
yapp2 = del2(y0)./ht; %matlab數(shù)值近似
k2 = abs(yapp2)./(1+yapp1.^2).^(3/2);
figure
plot(k2)
title('曲率曲線')
[~,maxFlag] = max(k2);%曲率最大位置
x_max = x0(maxFlag);
y_max = y0(maxFlag);
%畫(huà)出圖像 標(biāo)注曲率最大點(diǎn)
figure
plot(x0,y0,'.-');
hold on;
plot(x_max,y_max,'rp')
title('標(biāo)注最大曲率點(diǎn)')
xlabel('log10((norm(B*Xk-L)))')
ylabel('log10((norm(Xk)))')

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

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

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