基于matlab计算不等间距样本的一阶和二阶导数

%Demo for 'derivativeNE.m'%After 5 clicks on 'run' you get the results for sine curves%No plot for dy1c(x1c), since dy1c is nearly on top of dy1clcclose all%Create parbola or sineif exist('curveTypeCounter','var')  curveTypeCounter=curveTypeCounter+1;  if curveTypeCounter>10    curveTypeCounter=1;  endelse  clear  curveTypeCounter=1;endcurveType=(curveTypeCounter>5);%  curveType=1;%0=force parabola, 1=force sine%Create parbola or sinen=round(rand(1)*150+10);%Number of samples, min=10, max=160% n=20;%Force n samplesif curveType==0  %Parabola with spike at x=0  xr=rand(n,1)*2;  x=sort(xr)-1;  %create spike  indxL=find(x<-0.1,1,'last');  indxR=find(x>0.15,1,'first');  x=[x(1:indxL);-0.1;0;0.15;x(indxR:end)];  x(1)=-1;  x(end)=1;  y=x.^2;  y(indxL+2)=0.05;%Spike amplitude  legendText1=sprintf('y=x^{2} with spike at x=0');else  %Sine  xr=rand(n,1)*2*pi;  x=sort(xr);  x(1)=0;  x(end)=2*pi;  y=sin(x);  legendText1='y=sin(x)';end[dy1,dy2,x1c,dy1c]=derivativeNE(x,y);%Calculate derivatives%Plot resultsnl=char(10);%new line, compatible to older Matlab releaseshfig=figure(1);clfhfig.Position(3)=hfig.Position(4)*1.8;yyaxis leftp1=plot(x,y,'.-g');% original curvexlabel('x')if curveType==0  %Parabola  ylabel('y = x^2 and spike at x=0')  xticks([-1,-.5,-.1,0,.15,.5,1])  text(1.27,0.5,['The exact derivates' nl 'at x=-0.1, 0, 0.15' nl 'are not defined as numbers.' nl ...    'They are for the 1st derivative' nl 'Heaviside functions' nl ...    'and for the 2nd derivative' nl 'Dirac functions.'],'BackgroundColor',[1,.85,.85],'FontSize',8)  yyaxis right  xexact=-1:0.01:1;  dy1exact=2*xexact;  p2=plot(xexact(1:91),dy1exact(1:91),'color',[.7,.7,1],'LineWidth',4);%left part exact 1st derivative  hold on  indxx0=find(x==0);  dyL=(y(indxx0)-y(indxx0-1))/(x(indxx0)-x(indxx0-1));  dyR=(y(indxx0+1)-y(indxx0))/(x(indxx0+1)-x(indxx0));  plot(xexact([91,91,101,101]),[dy1exact(91),dyL,dyL,dyR],'-','color',[.7,.7,1],'LineWidth',4);%left part of spike exact 1st derivative  plot(xexact([101,116,116]),[dyR,dyR,dy1exact(116)],'-','color',[.7,.7,1],'LineWidth',4);%right part of spike exact 1st derivative  plot(xexact(116:end),dy1exact(116:end),'-','color',[.7,.7,1],'LineWidth',4);%right part exact 1st derivative  p3=plot(xexact([1,90]),[2,2],'-','color',[1,.7,.8],'LineWidth',4);%left part exact 2nd derivative  plot(xexact([92,115]),[0,0],'-','color',[1,.7,.8],'LineWidth',4);%left and right part of spike exact 2nd derivative  plot(xexact([117,end]),[2,2],'-','color',[1,.7,.8],'LineWidth',4);%right part exact 2nd derivative  %Plot Diracs  quiver(xexact(91),0,0,8,'-','Color',[1,0.7,0.7],'LineWidth',2.5,'MaxHeadSize',.06)  quiver(xexact(101),0,0,-6,'-','Color',[1,0.7,0.7],'LineWidth',2.5,'MaxHeadSize',.08)  quiver(xexact(116),0,0,8,'-','Color',[1,0.7,0.7],'LineWidth',2.5,'MaxHeadSize',.06)    p4=plot(x,dy1,'b-');%1st numerical derivative interpolated  % p4=plot(x1s,dy1c,'.y');%1st derivative centered  p5=plot(x,dy2,'+r-');%2nd numerical derivative  legend([p1,p4,p2,p5,p3],{sprintf('y=x^{2} with spike at x=0'),'1st derivative numerical','1st derivative exact',...    '2nd derivative numerical','2nd derivative exact'},'Location','northeastoutside')  ylim([-8,12])    %Calculate Mean Absolute Errors and max errors without range of spike  %1st derivative at x  x4MAEdy1=[x(1:indxx0-2);x(indxx0+2:end)];%     x for MAE  dy14MAE= [dy1(1:indxx0-2);dy1(indxx0+2:end)];% dy1 for MAE  dy1AD=abs(dy14MAE-2*x4MAEdy1);%                Absolute differences  [dy1Max,dy1Maxindx]=max(dy1AD);%               Max absolute difference of dy1-dy1exact  dy1MAE= sum(dy1AD)/numel(x4MAEdy1);%           MAE dy1    %1st derivative at x1c including spike  x4MAEdy1c=[x1c(1:indxx0-2);x1c(indxx0+1:end)];%    x for MAE  dy1c4MAE= [dy1c(1:indxx0-2);dy1c(indxx0+1:end)];%  dy1c for MAE  dy1cAD=abs(dy1c4MAE-2*x4MAEdy1c);%                 Absolute differences  [dy1cMax,dy1cMaxindx]=max(dy1cAD);%                Max absolute difference of dy1c-dy1exact  dy1cMAE= sum(dy1cAD)/numel(x4MAEdy1c);%            MAE dy1c    %2nd derivative  x4MAEdy2=[x(2:indxx0-3);x(indxx0+3:end-1)];%      x for dy2 MAE  dy24MAE=[dy2(2:indxx0-3);dy2(indxx0+3:end-1)];%   dy2 for MAE  dy2AD=abs(dy24MAE-2);%                            Absolute differences  [dy2Max,dy2Maxindx]=max(dy2AD);%                  Max absolute difference of dy2-dy2exact  dy2MAE=sum(dy2AD)/numel(x4MAEdy2);%               MAE dy2    %Output text for MAE  text(1.27,-6,['Mean Abs. Error     Max abs. error' nl ...     'dy1:  ' num2str(dy1MAE,'%.1e') '         ' num2str(dy1Max,2) ' @x=' num2str(x4MAEdy1(dy1Maxindx),2) nl ...     'dy1c: ' num2str(dy1cMAE,'%.1e') '         ' num2str(dy1cMax,2) ' @x=' num2str(x4MAEdy1c(dy1cMaxindx),2) nl ...    'dy2:  ' num2str(dy2MAE,'%.1e') '         ' num2str(dy2Max,2) ' @x=' num2str(x4MAEdy2(dy2Maxindx),2) nl ...    'Range around spike is excluded'], ...    'BackgroundColor',[0.98,.98,.98],'FontSize',8)else  %Sine  ylim([-1.05,1.05])  ylabel('y = sin(x)')  xticks(0:pi/2:2*pi+eps)  set(gca,'TickLabelInterpreter','latex','XTickLabel',{'0','$\frac{\pi}{2}$','$\pi$','$\frac{3\pi}{2}$','$2\pi$'})%   set(gca,'XMinorGrid','on')%Only for test  yyaxis right  xexact=linspace(0,2*pi,200);  p2=plot(xexact,cos(xexact),'color',[.8,.8,1],'LineWidth',4);%         exact 1st derivative  hold on  p3=plot(xexact,-sin(xexact),'-','color',[1,.8,0.8],'LineWidth',4);%   exact 2nd derivative  p4=plot(x,dy1,'b-');%                                                 1st numerical derivative interpolated  % p4=plot(x1s,dy1c,'y-');%1st derivative centered  p5=plot(x,dy2,'+r-');%                                                2nd numerical derivative  legend([p1,p4,p2,p5,p3],{'y = sin(x)','1st derivative numerical','1st derivative exact',...    '2nd derivative numerical','2nd derivative exact'},'Location','northeastoutside')  ylim([-1.05,1.05])  xlim([0,2*pi])    %Calculate Mean Absolute Errors and max errors   %1st derivative at x  dy1AD=abs(dy1-cos(x));%            Absolute differences  [dy1Max,dy1Maxindx]=max(dy1AD);%   Max absolute difference of dy1-dy1exact  dy1MAE= sum(dy1AD)/numel(x);%      MAE dy1    %1st derivative at x1c   dy1cAD=abs(dy1c-cos(x1c));%             Absolute differences  [dy1cMax,dy1cMaxindx]=max(dy1cAD);%     Max absolute difference of dy1c-dy1exact  dy1cMAE= sum(dy1cAD)/numel(x1c);%       MAE dy1c    %2n derivative  x4MAEdy2=x(2:end-1);%                x for dy2 MAE  dy24MAE=dy2(2:end-1);%               dy2 for MAE  dy2AD=abs(dy24MAE+sin(x4MAEdy2));%   Absolute differences  [dy2Max,dy2Maxindx]=max(dy2AD);%     Max absolute difference of dy2-dy2exact  dy2MAE=sum(dy2AD)/numel(x4MAEdy2);%  MAE dy2      %Output text for MAE    text(7.2,0,['Mean Abs. Error     Max abs. error' nl ...     'dy1:  ' num2str(dy1MAE,'%.1e') '         ' num2str(dy1Max,'%.1e') ' @x=' num2str(x(dy1Maxindx),2) nl ...     'dy1c: ' num2str(dy1cMAE,'%.1e') '         ' num2str(dy1cMax,'%.1e') ' @x=' num2str(x1c(dy1cMaxindx),2) nl ...    'dy2:  ' num2str(dy2MAE,'%.1e') '         ' num2str(dy2Max,'%.1e') ' @x=' num2str(x4MAEdy2(dy2Maxindx),2)], ...    'BackgroundColor',[0.98,.98,.98],'FontSize',8)endylabel('1st and 2nd derivatives')hold offgrid on

function [dy1,dy2,x1c,dy1c]=derivativeNE(x,y)%Calculate the 1st and 2nd derivative for Not Equal spaced samples.%%INPUT%  x as column or row vector, more than 2 values%  y as column or row vector, same amount of values as for x%%OUTPUT%  dy1: 1st derivative, y'(x)%       %  dy2: 2nd derivative, y"(x)%       The first dy2 at x(1) and the last dy2 at x(end) don't exist and are set to NaN%  x1c and dy1c return a more accurate (about 3 times better) 1st derivative dy1c(x1c), but with new x-values%%  If x or y is a row vector, all outputs are returned as row vectors.%  On error NaN is returned%%Syntax examples:%  dy1=derivativeNE([1;3;4;5],[1;9;16;25]);           % 1st derivative y'(x)%  [~,dy2]=derivativeNE([1;3;4;5],[1;9;16;25]);       % 2nd derivative y"(x)%  [~,~,x1c,y1c]=derivativeNE([1;3;4;5],[1;9;16;25]); % 1st derivative more accurate (about 3 times), y'(x1c)%%Peter Seibold, August 2023%Check inputsif size(x,2)>1 || size(y,2)>1  %Either x or y is a row vector  isRowVector=true;  x=x(:);%convert to column vector  y=y(:);else  isRowVector=false;endif numel(x)~=numel(y) || numel(x)<3  dy1=NaN; dy2=NaN; x1c=NaN; dy1c=NaN;%Return all output as NaN  disp('x and y must have the same size and at least 3 elements!')  returnend%1st derivativedx=diff(x);dy=diff(y);dy1c=dy./dx;%          1st derivative for centered xx1c=x(1:end-1)+dx/2;%  center each x (shift each x by half distance to next x)%interpolate/extrapolate dy1c(x1c) to match x, result: dy1(x)y2=[dy1c(2);dy1c(2:end);dy1c(end)];   %repeat border values for extrapolationy1=[dy1c(1);dy1c(1:end-1);dy1c(end-1)];x1=[x1c(1);x1c(1:end-1);x1c(end-1)];x2=[x1c(2);x1c(2:end);x1c(end)];dy1=(y2-y1)./(x2-x1).*(x-x2)+y2;      %interpolate/extrapolate%2nd derivativey32x32=dy1c(2:end);                   %1st derivatives at even positionsy21x21=dy1c(1:end-1);                 %1st derivatives at odd positionsdx31=y32x32;                          %create vector with dummy values, with two elements less than number of x-elementsdx31(1:2:end)=x(3:2:end)-x(1:2:end-2);%x-difference odd pairsdx31(2:2:end)=x(4:2:end)-x(2:2:end-2);%x-difference even pairsdy2c=2*(y32x32-y21x21)./dx31;         %2nd derivatives dy2 at center of the two surrounding samplesx2c=dy2c;                             %create vector with dummy values, with two elements less than number of x-elementsx2c(1:2:end)=(x(1:2:end-2)+x(3:2:end))/2;%x for dyc at center of the two surrounding samples, (overwrite dummy values)x2c(2:2:end)=(x(2:2:end-2)+x(4:2:end))/2;%dito for even x centered, now we have all center x2c for 2nd derivatives%interpolate dy2c(x2c) to match x, result: dy2(x)if numel(x)>3  dy2=[NaN;interp1(x2c,dy2c,x(2:end-1),'linear','extrap');NaN];%interpolate and pad unavailable border valueselse  %Interpolation needs at least two dy2c samples  dy2=[NaN;dy2c;NaN];%    output only the single dy2(x)end%Convert back to user formatif isRowVector  %Either x or y input was a row vector  %Convert to row vector  dy1=dy1';  dy2=dy2';  x1c=x1c';  dy1c=dy1c';end