Calculate the first and second derivatives of unequal spacing samples based on matlab

?Author’s brief introduction: A Matlab simulation developer who loves scientific research. He cultivates his mind and technology at the same time. Matlab project cooperation can be privately messaged.

Personal homepage: Matlab Research Studio

Personal credo: Investigate things to gain knowledge.

For more Matlab simulation content click

Intelligent optimization algorithm Neural network prediction Radar communication Wireless sensor Power system

Signal Processing Image Processing Path Planning Cellular Automata UAV

? Content introduction

1.1 Grammar

[dy1,dy2,x1c,dy1c]=derivative NE(x,y)

1.2 Input

x as column or row vector, more than 2 values

y as a column or row vector with the same number of values as x

1.3 output

dy1: first derivative

dy2: second derivative

The first dy2 at x(1) and the last dy2 at x(end) do not exist and are set to NaN

x1c and dy1c return the exact first derivative dy1c(x1c)

If x or y is a row vector, all output is returned as a row vector.

On error, all output returns NaN.

1.4 Grammar Example

For the first derivative:

dy1=derivative NE([1;3;4;5],[1;9;16;25]);

For the second derivative:

[~,dy2]=derivative NE([1;3;4;5],[1;9;16;25]);

For more accurate first derivatives:

[~,~,x1c,y1c]=derivativeNE([1;3;4;5],[1;9;16;25]);

? Complete code

?mo for 'derivativeNE.m'</code><code>?ter 5 clicks on 'run' you get the results for sine curves</code><code>%No plot for dy1c(x1c), since dy1c is nearly on top of dy1</code><code>clc</code><code>close all</code><code>%Create parbola or sine</code><code>if exist('curveTypeCounter', 'var')</code><code> curveTypeCounter=curveTypeCounter + 1;</code><code> if curveTypeCounter>10</code><code> curveTypeCounter=1;</code><code> end</code><code>else</code><code> clear</code><code> curveTypeCounter=1;</code><code>end</code><code>curveType=(curveTypeCounter>5);</code><code>% curveType=1;%0=force parabola, 1=force sine</code><code>%Create parbola or sine</code><code>n=round(rand(1)*150 + 10 );%Number of samples, min=10, max=160</code><code>% n=20;%Force n samples</code><code>if curveType==0</code><code> % Parabola with spike at x=0</code><code> xr=rand(n,1)*2;</code><code> x=sort(xr)-1;</code><code> %create spike</code><code> indxL=find(x<-0.1,1,'last');</code><code> indxR=find(x>0.15,1,'first');</code> <code> x=[x(1:indxL);-0.1;0;0.15;x(indxR:end)];</code><code> x(1)=-1;</code><code> x(end)=1;</code><code> y=x.^2;</code><code> y(indxL + 2)=0.05;%Spike amplitude</code><code> legendText1=sprintf ('y=x^{2} with spike at x=0');</code><code>else</code><code> %Sine</code><code> xr=rand(n,1) *2*pi;</code><code> x=sort(xr);</code><code> x(1)=0;</code><code> x(end)=2*pi;</code><code> y=sin(x);</code><code> legendText1='y=sin(x)';</code><code>end</code><code>?</code><code>[dy1,dy2,x1c,dy1c]=derivativeNE(x,y);?lculate derivatives</code><code>?</code><code>%Plot results</code><code>nl =char(10);%new line, compatible to older Matlab releases</code><code>hfig=figure(1);</code><code>clf</code><code>hfig.Position(3) =hfig.Position(4)*1.8;</code><code>yyaxis left</code><code>p1=plot(x,y,'.-g');% original curve</code><code>xlabel('x')</code><code>if curveType==0</code><code> %Parabola</code><code> ylabel('y = x^2 and spike at x=0' )</code><code> xticks([-1,-.5,-.1,0,.15,.5,1])</code><code> text(1.27,0.5,['The exact derivates' nl 'at x=-0.1, 0, 0.15' nl 'are not defined as numbers.' nl ...</code><code> 'They are for the 1st derivative' nl 'Heaviside functions' nl .. .</code><code> 'and for the 2nd derivative' nl 'Dirac functions.'],'BackgroundColor',[1,.85,.85],'FontSize',8)</code><code> yyaxis right</code><code> xexact=-1:0.01:1;</code><code> dy1exact=2*xexact;</code><code> p2=plot(xexact(1:91),dy1exact (1:91),'color',[.7,.7,1],'LineWidth',4);%left part exact 1st derivative</code><code> hold on</code><code> indxx0 =find(x==0);</code><code> dyL=(y(indxx0)-y(indxx0-1))/(x(indxx0)-x(indxx0-1));</code> <code> dyR=(y(indxx0 + 1)-y(indxx0))/(x(indxx0 + 1)-x(indxx0));</code><code> plot(xexact([91,91,101,101]) ,[dy1exact(91),dyL,dyL,dyR],'-','color',[.7,.7,1],'LineWidth',4);%left part of spike exact 1st derivative</code><code> plot(xexact([101,116,116]),[dyR,dyR,dy1exact(116)],'-','color',[.7,.7,1],'LineWidth',4);% right part of spike exact 1st derivative</code><code> plot(xexact(116:end),dy1exact(116:end),'-','color',[.7,.7,1],'LineWidth ',4);%right part exact 1st derivative</code><code> p3=plot(xexact([1,90]),[2,2],'-','color',[1,.7 ,.8],'LineWidth',4);%left part exact 2nd derivative</code><code> plot(xexact([92,115]),[0,0],'-','color',[1 ,.7,.8],'LineWidth',4);%left and right part of spike exact 2nd derivative</code><code> plot(xexact([117,end]),[2,2],' -','color',[1,.7,.8],'LineWidth',4);%right part exact 2nd derivative</code><code> %Plot Diracs</code><code> quiver(xexact (91),0,0,8,'-','Color',[1,0.7,0.7],'LineWidth',2.5,'MaxHeadSize',.06)</code><code> quiver(xexact( 101),0,0,-6,'-','Color',[1,0.7,0.7],'LineWidth',2.5,'MaxHeadSize',.08)</code><code> quiver(xexact( 116),0,0,8,'-','Color',[1,0.7,0.7],'LineWidth',2.5,'MaxHeadSize',.06)</code><code> </code><code> p4=plot(x,dy1,'b-');%1st numerical derivative interpolated</code><code> % p4=plot(x1s,dy1c,'.y');%1st derivative centered</code><code> p5=plot(x,dy2,' + r-');%2nd numerical derivative</code><code> legend([p1,p4,p2,p5,p3],{sprintf('y= x^{2} with spike at x=0'),'1st derivative numerical','1st derivative exact',...</code><code> '2nd derivative numerical','2nd derivative exact'},' Location','northeastoutside')</code><code> ylim([-8,12])</code><code> </code><code> ?lculate Mean Absolute Errors and max errors without range of spike</code><code> %1st derivative at x</code><code> x4MAEdy1=[x(1:indxx0-2);x(indxx0 + 2:end)];% x for MAE</code><code> dy14MAE= [dy1(1:indxx0-2);dy1(indxx0 + 2:end)];% dy1 for MAE</code><code> dy1AD=abs(dy14MAE-2*x4MAEdy1);% Absolute differences</code><code> [dy1Max,dy1Maxindx]=max(dy1AD);% Max absolute difference of dy1-dy1exact</code><code> dy1MAE= sum(dy1AD)/numel(x4MAEdy1);% MAE dy1</code> <code> </code><code> %1st derivative at x1c including spike</code><code> x4MAEdy1c=[x1c(1:indxx0-2);x1c(indxx0 + 1:end)];% x for MAE </code><code> dy1c4MAE= [dy1c(1:indxx0-2);dy1c(indxx0 + 1:end)];% dy1c for MAE</code><code> dy1cAD=abs(dy1c4MAE-2*x4MAEdy1c) ;% Absolute differences</code><code> [dy1cMax,dy1cMaxindx]=max(dy1cAD);% Max absolute difference of dy1c-dy1exact</code><code> dy1cMAE= sum(dy1cAD)/numel(x4MAEdy1c);% MAE dy1c</code><code> </code><code> %2nd derivative</code><code> x4MAEdy2=[x(2:indxx0-3);x(indxx0 + 3:end-1)]; % x for dy2 MAE</code><code> dy24MAE=[dy2(2:indxx0-3);dy2(indxx0 + 3:end-1)];% dy2 for MAE</code><code> dy2AD=abs (dy24MAE-2);% Absolute differences</code><code> [dy2Max,dy2Maxindx]=max(dy2AD);% Max absolute difference of dy2-dy2exact</code><code> dy2MAE=sum(dy2AD)/numel (x4MAEdy2);% MAE dy2</code><code> </code><code> %Output text for MAE</code><code> text(1.27,-6,['Mean Abs. Error Max abs. error ' nl ...</code><code> 'dy1: ' num2str(dy1MAE,'%.1e') ' ' num2str(dy1Max,2) ' @x=' num2str(x4MAEdy1(dy1Maxindx),2) nl . ..</code><code> 'dy1c: ' num2str(dy1cMAE,'%.1e') ' ' num2str(dy1cMax,2) ' @x=' num2str(x4MAEdy1c(dy1cMaxindx),2) nl ...</code><code> 'dy2: ' num2str(dy2MAE,'%.1e') ' ' num2str(dy2Max,2) ' @x=' num2str(x4MAEdy2(dy2Maxindx),2) nl ...</code> <code> 'Range around spike is excluded'], ...</code><code> 'BackgroundColor',[0.98,.98,.98],'FontSize',8)</code><code>? </code><code>else</code><code> %Sine</code><code> ylim([-1.05,1.05])</code><code> ylabel('y = sin(x)' )</code><code> xticks(0:pi/2:2*pi + eps)</code><code> set(gca,'TickLabelInterpreter','latex','XTickLabel',{'0', '$\frac{\pi}{2}$','$\pi$','$\frac{3\pi}{2}$','$2\pi$'})</code><code>% set(gca,'XMinorGrid','on')%Only for test</code><code> yyaxis right</code><code> xexact=linspace(0,2*pi,200);</code><code> p2=plot(xexact,cos(xexact),'color',[.8,.8,1],'LineWidth',4);% exact 1st derivative</code><code> hold on</code><code> p3=plot(xexact,-sin(xexact),'-','color',[1,.8,0.8],'LineWidth',4);% exact 2nd derivative</code> <code> p4=plot(x,dy1,'b-');% 1st numerical derivative interpolated</code><code> % p4=plot(x1s,dy1c,'y-');% 1st derivative centered</code> code><code> p5=plot(x,dy2,' + r-');% 2nd numerical derivative</code><code> legend([p1,p4,p2,p5,p3],{'y = sin (x)','1st derivative numerical','1st derivative exact',...</code><code> '2nd derivative numerical','2nd derivative exact'},'Location','northeastoutside')</code><code> ylim([-1.05,1.05])</code><code> xlim([0,2*pi])</code><code> </code><code> ?lculate Mean Absolute Errors and max errors </code><code> %1st derivative at x</code><code> dy1AD=abs(dy1-cos(x));% Absolute differences</code><code> [dy1Max,dy1Maxindx]= max(dy1AD);% Max absolute difference of dy1-dy1exact</code><code> dy1MAE= sum(dy1AD)/numel(x);% MAE dy1</code><code> </code><code> % 1st derivative at x1c </code><code> dy1cAD=abs(dy1c-cos(x1c));% Absolute differences</code><code> [dy1cMax,dy1cMaxindx]=max(dy1cAD);% Max absolute difference of dy1c -dy1exact</code><code> dy1cMAE= sum(dy1cAD)/numel(x1c);% MAE dy1c</code><code> </code><code> %2n derivative</code><code> x4MAEdy2= x(2:end-1);% x for dy2 MAE</code><code> dy24MAE=dy2(2:end-1);% dy2 for MAE</code><code> dy2AD=abs(dy24MAE + sin (x4MAEdy2));% Absolute differences</code><code> [dy2Max,dy2Maxindx]=max(dy2AD);% Max absolute difference of dy2-dy2exact</code><code> dy2MAE=sum(dy2AD)/numel( x4MAEdy2);% MAE dy2</code><code> </code><code> </code><code> %Output text for MAE</code><code> text(7.2,0,['Mean Abs. Error Max abs. error' nl ...</code><code> 'dy1: ' num2str(dy1MAE,'%.1e') ' ' num2str(dy1Max,'%.1e') ' @x=' num2str( x(dy1Maxindx),2) nl ...</code><code> 'dy1c: ' num2str(dy1cMAE,'%.1e') ' ' num2str(dy1cMax,'%.1e') ' @x=' num2str (x1c(dy1cMaxindx),2) nl ...</code><code> 'dy2: ' num2str(dy2MAE,'%.1e') ' ' num2str(dy2Max,'%.1e') ' @x=' num2str(x4MAEdy2(dy2Maxindx),2)], ...</code><code> 'BackgroundColor',[0.98,.98,.98],'FontSize',8)</code><code>?</code><code>end</code><code>ylabel('1st and 2nd derivatives')</code><code>hold off</code><code>grid on</code><code>?< /pre>
<pre>function [dy1,dy2,x1c,dy1c]=derivativeNE(x,y)</code><code>?lculate the 1st and 2nd derivative for Not Equal spaced samples.</code><code>%</code><code> code><code>%INPUT</code><code>% x as column or row vector, more than 2 values</code><code>% y as column or row vector, same amount of values as for x</code><code> code><code>%</code><code>%OUTPUT</code><code>% dy1: 1st derivative, y'(x)</code><code>% </code><code>% dy2 : 2nd derivative, y"(x)</code><code>% The first dy2 at x(1) and the last dy2 at x(end) don't exist and are set to NaN</code><code> % x1c and dy1c return a more accurate (about 3 times better) 1st derivative dy1c(x1c), but with new x-values</code><code>%</code><code>% If x or y is a row vector, all outputs are returned as row vectors.</code><code>% On error NaN is returned</code><code>%</code><code>%Syntax examples:</code><code>% dy1=derivativeNE([1;3;4;5],[1;9;16;25]); % 1st derivative y'(x)</code><code>% [~,dy2]=derivativeNE([ 1;3;4;5],[1;9;16;25]); % 2nd derivative y"(x)</code><code>% [~,~,x1c,y1c]=derivativeNE([1 ;3;4;5],[1;9;16;25]); % 1st derivative more accurate (about 3 times), y'(x1c)</code><code>%</code><code> %Peter Seibold, August 2023</code><code>?</code><code>%Check inputs</code><code>if size(x,2)>1 || size(y,2)>1 </code><code> %Either x or y is a row vector</code><code> isRowVector=true;</code><code> x=x(:);%convert to column vector</code> <code> y=y(:);</code><code>else</code><code> isRowVector=false;</code><code>end</code><code>if numel(x)~ =numel(y) || numel(x)<3</code><code> dy1=NaN; dy2=NaN; x1c=NaN; dy1c=NaN;%Return all output as NaN</code><code> disp ('x and y must have the same size and at least 3 elements!')</code><code> return</code><code>end</code><code>?</code><code>% 1st derivative</code><code>dx=diff(x);</code><code>dy=diff(y);</code><code>dy1c=dy./dx;% 1st derivative for centered x </code><code>x1c=x(1:end-1) + dx/2;% center each x (shift each x by half distance to next x)</code><code>?</code><code>%interpolate/extrapolate dy1c(x1c) to match x, result: dy1(x)</code><code>y2=[dy1c(2);dy1c(2:end);dy1c(end)]; %repeat border values for extrapolation</code><code>y1=[dy1c(1);dy1c(1:end-1);dy1c(end-1)];</code><code>x1=[x1c(1) ;x1c(1:end-1);x1c(end-1)];</code><code>x2=[x1c(2);x1c(2:end);x1c(end)];</code> <code>dy1=(y2-y1)./(x2-x1).*(x-x2) + y2; %interpolate/extrapolate</code><code>?</code><code>%2nd derivative</code><code>y32x32=dy1c(2:end); %1st derivatives at even positions</code><code>y21x21=dy1c(1:end-1); %1st derivatives at odd positions</code><code>dx31=y32x32; %create vector with dummy values, with two elements less than number of x-elements</code><code>dx31(1:2:end)=x(3:2:end)-x( 1:2:end-2);%x-difference odd pairs</code><code>dx31(2:2:end)=x(4:2:end)-x(2:2:end-2) ;%x-difference even pairs</code><code>dy2c=2*(y32x32-y21x21)./dx31; %2nd derivatives dy2 at center of the two surrounding samples</code><code>x2c=dy2c; % create vector with dummy values, with two elements less than number of x-elements</code><code>x2c(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)</code><code>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</code><code>?</code><code>%interpolate dy2c(x2c) to match x, result: dy2(x)</code><code>if numel(x)>3</code><code> dy2=[NaN;interp1(x2c,dy2c,x(2:end-1) ,'linear','extrap');NaN];%interpolate and pad unavailable border values</code><code>else</code><code> %Interpolation needs at least two dy2c samples</code><code> dy2=[NaN;dy2c;NaN];% output only the single dy2(x)</code><code>end</code><code>?</code><code>%Convert back to user format</code> code><code>if isRowVector</code><code> %Either x or y input was a row vector</code><code> %Convert to row vector</code><code> dy1=dy1';</code> code><code> dy2=dy2';</code><code> x1c=x1c';</code><code> dy1c=dy1c';</code><code>end</code><code>? </code><code>?</code><code>?

? Run results

? code to get follow me

Some theories are quoted from online literature. If there is any infringement, please contact the blogger to delete it

Follow me to receive massive matlab e-books and mathematical modeling materials

Simulation consultation

1 Improvement and application of various intelligent optimization algorithms

Production scheduling, economic scheduling, assembly line scheduling, charging optimization, workshop scheduling, departure optimization, reservoir scheduling, three-dimensional packing, logistics location selection, cargo space optimization, bus scheduling optimization, charging pile layout optimization, workshop layout optimization, Container ship stowage optimization, water pump combination optimization, medical resource allocation optimization, facility layout optimization, visible area base station and drone site selection optimization

2 Machine learning and deep learning

Convolutional Neural Network (CNN), LSTM, Support Vector Machine (SVM), Least Squares Support Vector Machine (LSSVM), Extreme Learning Machine (ELM), Kernel Extreme Learning Machine (KELM), BP, RBF, Width Learning, DBN, RF, RBF, DELM, XGBOOST, TCN realize wind power forecasting, photovoltaic forecasting, battery life forecasting, radiation source identification, traffic flow forecasting, load forecasting, stock price forecasting, PM2.5 concentration forecasting, battery health status forecasting, water bodies Optical parameter inversion, NLOS signal identification, accurate subway parking prediction, transformer fault diagnosis

2. Image processing

Image recognition, image segmentation, image detection, image hiding, image registration, image stitching, image fusion, image enhancement, image compression perception

3 Path planning

Traveling salesman problem (TSP), vehicle routing problem (VRP, MVRP, CVRP, VRPTW, etc.), UAV three-dimensional path planning, UAV collaboration, UAV formation, robot path planning, raster map path planning , multimodal transportation problem, vehicle cooperative UAV path planning, antenna linear array distribution optimization, workshop layout optimization

4 UAV applications

UAV path planning, UAV control, UAV formation, UAV coordination, UAV task assignment

5 Wireless sensor positioning and layout

Sensor deployment optimization, communication protocol optimization, routing optimization, target positioning optimization, Dv-Hop positioning optimization, Leach protocol optimization, WSN coverage optimization, multicast optimization, RSSI positioning optimization

6 Signal processing

Signal recognition, signal encryption, signal denoising, signal enhancement, radar signal processing, signal watermark embedding and extraction, EMG signal, EEG signal, signal timing optimization