copy & paste로 매틀랩 파일로 저장한후 실행시키세요. 함수로 짜여있으니 파일명은 derivatives.m으로 하시고 실행은 derivatives(1) ~ derivatives(4) 중 하나 선택하셔서 해보세요. 1번은 베셀, 2번은 사인, 3번은 포물선, 4번은 원의 1,2,3차 미분값을 보여준답니다. :)
=======================================================



function []=Derivatives(input_num)
% made by Dong-Youn Shin, UMIST
% 2002-11-08
% You can use this function free but you must keep the above lines

close all
clc
figure;

switch input_num
   
    case 1
        disp('Bessel');
        % ================== Bessel ==================
        x=linspace(0,20,50);y_temp=besselj(0,x);
        xx=linspace(min(x),max(x),max(size(x))*2-1);
        h=xx(2)-xx(1);
        % Would you like to put some noise on your data?
        % Then input a small number instead of 0.0 below
        y=y_temp+0.0*(rand(size(x))-.5);
       
        for loop=1:1:max(size(xx))
            a_dydx(loop)=-besselj(1,xx(loop));
            a_d2ydx2(loop)=-besselj(0,xx(loop))+besselj(1,xx(loop))/xx(loop);
            a_d3ydx3(loop)=besselj(1,xx(loop))-2*besselj(1,xx(loop))/xx(loop)^2+besselj(0,xx(loop))/xx(loop);
        end
       
    case 2
        disp('Sinusoidal');
        % ================== sinusoidal==================
        x=linspace(0,2*pi,50);y_temp=sin(x);
        xx=linspace(min(x),max(x),max(size(x))*2-1);
        h=xx(2)-xx(1);
        % Would you like to put some noise on your data?
        % Then input a small number instead of 0.0 below
        y=y_temp+0.0*(rand(size(x))-.5);
       
        for loop=1:1:max(size(xx))
            a_dydx(loop)=cos(xx(loop));
            a_d2ydx2(loop)=-sin(xx(loop));
            a_d3ydx3(loop)=-cos(xx(loop));
        end
       
    case 3
        disp('Parabolic');
        % ================== Parabolic ==================
        r=10;
        Multiplier=1;
        x=linspace(0,r*0.99,100);
        xx=linspace(min(x),max(x),max(size(x))*2-1);
        h=xx(2)-xx(1);
        for loop=1:1:max(size(x));
            y(loop)=sqrt(r-x(loop));
        end
       
        for loop=1:1:max(size(xx));
            a_dydx(loop)=-1/(2*sqrt(r-xx(loop)));
            a_d2ydx2(loop)=1/(4*(xx(loop)-r)*sqrt(r-xx(loop)));
            a_d3ydx3(loop)=-3/(8*(xx(loop)-r)^2*sqrt(r-xx(loop)));
        end
       
    case 4       
        % ================== Circle ==================
        r=10;
        Multiplier=1;
        x=linspace(0,2*r*0.99,100);
        xx=linspace(min(x),max(x),max(size(x))*2-1);
        h=xx(2)-xx(1);
       
        for loop=1:1:max(size(x));
            if x(loop)<r
                y(loop)=r;
            else
                y(loop)=sqrt(r^2-(x(loop)-r)^2);
            end
        end
       
        for loop=1:1:max(size(xx));
            if xx(loop)<r
                a_dydx(loop)=0;
                a_d2ydx2(loop)=0;
                a_d3ydx3(loop)=0;
            else
                m_x=xx(loop)-r;
                a_dydx(loop)=-m_x/sqrt(r^2-m_x^2);
                a_d2ydx2(loop)=r^2/((m_x^2-r^2)*sqrt(r^2-m_x^2));
                a_d3ydx3(loop)=-3*m_x*r^2/((m_x^2-r^2)^2*sqrt(r^2-m_x^2));
            end
        end
       
    otherwise
        disp('You must choose 1 to 4');
        close all;
        return;
end


start=cputime;
y1=interp1(x,y,xx,'pchip');

subplot(2,2,1);
hold on;
plot(xx,y1,'-or','MarkerSize',2,'MarkerFaceColor','r');
plot(x,y,'ob','MarkerSize',2,'MarkerFaceColor','b');
legend('pchip - Noisy',0);
hold off;

finish=cputime-start;
disp('interp1 takes'),disp(finish)


start=cputime;
y1=spline(x,y,xx);

subplot(2,2,2);
hold on;
plot(xx,y1,'-or','MarkerSize',2,'MarkerFaceColor','r');
plot(x,y,'ob','MarkerSize',2,'MarkerFaceColor','b');
legend('spline - Noisy',0);
hold off;

finish=cputime-start;
disp('spline takes'),disp(finish)


start=cputime;
[y1,y2]=f_CurveSmoothing(x,y,xx,60,0);

subplot(2,2,3);
hold on;
plot(xx,y2,'-or','MarkerSize',2,'MarkerFaceColor','r');
plot(x,y,'ob','MarkerSize',2,'MarkerFaceColor','b');
legend('CurveSmoothing - Noisy',0);
hold off;

finish=cputime-start;
disp('CurveSmoothing takes'),disp(finish)


start=cputime;
[y1,p] = csaps(x,y,0.99,xx);


subplot(2,2,4);
hold on;
plot(xx,y1,'-or','MarkerSize',2,'MarkerFaceColor','r');
plot(x,y,'ob','MarkerSize',2,'MarkerFaceColor','b');
legend('csaps - Noisy',0);
hold off;

finish=cputime-start;
disp('csaps takes'),disp(finish)


% Polynomial fitting numerical derivatives
% If there is a noise, then activate the lines below deactivated by '%'
n_dydx=f_deriv_new(xx,y1,1);
% [n_dydx,p]=csaps(xx,n_dydx,0.99,xx);
n_d2ydx2=f_deriv_new(xx,n_dydx,1);
% [n_d2ydx2,p]=csaps(xx,n_d2ydx2,0.99,xx);
n_d3ydx3=f_deriv_new(xx,n_d2ydx2,1);
% [n_d3ydx3,p]=csaps(xx,n_d3ydx3,0.99,xx);


% 5 point numerical derivatives
% If there is a noise, then activate the lines below deactivated by '%'
o_dydx=f_deriv_old(y1,h);
% [o_dydx,p]=csaps(xx,o_dydx,0.99,xx);
o_d2ydx2=f_deriv_old(o_dydx,h);
% [o_d2ydx2,p]=csaps(xx,o_d2ydx2,0.99,xx);
o_d3ydx3=f_deriv_old(o_d2ydx2,h);
% [o_d3ydx3,p]=csaps(xx,o_d3ydx3,0.99,xx);





figure;
subplot(3,1,1);
hold on;
plot(xx,n_dydx,'-or','MarkerSize',2,'MarkerFaceColor','r');
plot(xx,o_dydx,'-ob','MarkerSize',2,'MarkerFaceColor','b');
plot(xx,a_dydx,'-ok','MarkerSize',2,'MarkerFaceColor','k');
legend('Polynomial fitted','5 point','Analytic',0);
title('1st derivative');
hold off;

subplot(3,1,2);
hold on;
plot(xx,n_d2ydx2,'-or','MarkerSize',2,'MarkerFaceColor','r');
plot(xx,o_d2ydx2,'-ob','MarkerSize',2,'MarkerFaceColor','b');
plot(xx,a_d2ydx2,'-ok','MarkerSize',2,'MarkerFaceColor','k');
legend('Polynomial fitted','5 point','Analytic',0);
title('2nd derivative');
hold off;

subplot(3,1,3);
hold on;
plot(xx,n_d3ydx3,'-or','MarkerSize',2,'MarkerFaceColor','r');
plot(xx,o_d3ydx3,'-ob','MarkerSize',2,'MarkerFaceColor','b');
plot(xx,a_d3ydx3,'-ok','MarkerSize',2,'MarkerFaceColor','k');
legend('Polynomial fitted','5 point','Analytic',0);
title('3rd derivative');
hold off;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Numerical derivative using polynomial curve fitting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydx=f_deriv_new(x,y,m_Multiplier)
% Numerical derivative using polynomial curve fitting
% made by Dong-Youn Shin, UMIST
% 2002-11-08
% You can use this function free but you must keep the above
% lines

x=x*m_Multiplier;
y=y*m_Multiplier;

% 5 point polynomial curve fitting
N=5;
Break=1;
array_size=max(size(y));
dydx=zeros(1,array_size);

for loop=3:3:array_size
    p=polyfit(x(loop-2:loop+2),y(loop-2:loop+2),2);
    k=polyder(p);
    dydx_temp=polyval(k,x(loop-2:loop+2));

    if loop==3
        dydx(loop-2:loop+1)=dydx_temp(loop-2:loop+1);
    elseif loop+2==array_size
        dydx(loop-1:loop+1)=dydx_temp(2:4);
       
        p=polyfit(x(loop-2:loop+2),y(loop-2:loop+2),2);
        k=polyder(p);
        dydx_temp=polyval(k,x(loop-2:loop+2));

        dydx(loop+2)=dydx_temp(5);
        break;
    elseif loop+4==array_size
        dydx(loop-1:loop+1)=dydx_temp(2:4);
       
        loop=loop+2;
        p=polyfit(x(loop-2:loop+2),y(loop-2:loop+2),2);
        k=polyder(p);
        dydx_temp=polyval(k,x(loop-2:loop+2));

        dydx(loop:loop+2)=dydx_temp(3:5);
        break;
    elseif loop+3==array_size
        dydx(loop-1:loop+1)=dydx_temp(2:4);
       
        loop=loop+1;
        p=polyfit(x(loop-2:loop+2),y(loop-2:loop+2),2);
        k=polyder(p);
        dydx_temp=polyval(k,x(loop-2:loop+2));

        dydx(loop+1:loop+2)=dydx_temp(4:5);
        break;
    else
        dydx(loop-1:loop+1)=dydx_temp(2:4);
    end
end
   
y=y/m_Multiplier;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Numerical derivative using 5 points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  V = f_deriv_old(U,h)
% get derivatives
% U = y values
% h = uniform step size


% DERIV Column-wise derivative estimation.
%        Computes 5-point discrete derivative estimates for each column
%        of the input matrix U.
%        V = DERIV(U)  Returns a matrix V of the same size as U each
%        element containing the estimates of the column-wise derivative of
%        a function sampled at unit intervals in U.
%
%        DERIV Uses 5 points both in the interior and at the edges.
%        If size of the matrix U is less than 5 use simpler 3-point
%        or 2-point estimate.
%
%  If you want to go really fancy, you may even use a 7-point formulae:
%  [-1/60  3/20  -3/4  0    3/4  -3/20  1/60] - in the interior,
%  [-49/20  6    -7.5  20/3  -15/4  6/5    1/6]  - for edge point
%
%        See also DIFF, GRADIENT, DEL2
%  Kirill K. Pankratov,  kirill@plume.mit.edu
%  March 22, 1994

 % 5-point coefficients for derivatives ..................................
 % Edges (1st, 2nd pts):
Cde = [-25/12 4 -3 4/3 -1/4; -1/4 -5/6 3/2 -1/2 1/12];
 % Interior
Cd = [1/12 -2/3 0 2/3 -1/12];

 % Determine the size of the input and make column if vector .............
ist = 0;
ly = size(U,1);
if ly==1, ist = 1; U = U(:); ly = length(U); end

 % If only 2 points - simple difference ..................................
if ly==2
  V(1,:) = U(2,:)-U(1,:);
  V(2,:) = V(1,:);
  if ist, V = V'; end    % Transpose output if necessary
  return
end

 % Now if more than 2 points - more complicated procedure ...............
if ly<5        % If less than 5 points

  V(2:ly-1,:) = (U(3:ly,:)-U(1:ly-2,:))/2;        % First
  V(1,:) = 2*U(2,:)-1.5*U(1,:)-.5*U(3,:);        % Last
  V(ly,:) = 1.5*U(ly,:)-2*U(ly-1,:)+.5*U(ly-2,:); % Interior

else            % If 5 points or more - regular, more accurate estimation

  % Edges ............
 V(1:2,:) = Cde*U(1:5,:);
  V([ly ly-1],:) = -fliplr(Cde)*U(ly-4:ly,:);

  % Interior .........
  V(3:ly-2,:) = Cd(1)*U(1:ly-4,:)+Cd(2)*U(2:ly-3,:)+Cd(3)*U(3:ly-2,:);
  V(3:ly-2,:) = V(3:ly-2,:)+Cd(4)*U(4:ly-1,:)+Cd(5)*U(5:ly,:);

end

V=V/h;

if ist, V = V'; end    % Transpose output if necessary



% Curve smoothing function
function [y,yy] = f_CurveSmoothing(x,y,xx,m_Error,m_drawoption)
% Curve smoothing function
% made by Dong-Youn Shin, UMIST
% 2002-11-08
% You can use this function free but you must keep the above lines

Break=1;   
while Break==1
    yold=y;
   
    for loop=3:3:max(size(x))-2
        p=polyfit(x(loop-2:loop+2),y(loop-2:loop+2),2);
        ytemp=polyval(p,x(loop-2:loop+2));
        y(loop-1:loop+1)=ytemp(2:4);
       
        for loop1=1:1:max(size(xx))
            if x(loop-2)<=xx(loop1)
                Min=loop1;
                break;
            end
        end
        for loop1=max(size(xx)):-1:1
            if x(loop+2)>=xx(loop1)
                Max=loop1;
                break;
            end
        end
       
        ytemp=polyval(p,xx(Min:Max));
        yy(Min:Max)=ytemp;
                       
        if m_drawoption==1
            set(gcf,'DoubleBuffer','on');
            plot(x,y,'-ob',x(loop-1:loop+1),y(loop-1:loop+1),'-*r');
            drawnow;
            pause(0.1);
        end
    end
   
    for loop=max(size(x))-2:-2:3
        p=polyfit(x(loop-2:loop+2),y(loop-2:loop+2),2);
        ytemp=polyval(p,x(loop-2:loop+2));
        y(loop-1:loop+1)=ytemp(2:4);
       
        for loop1=1:1:max(size(xx))
            if x(loop-2)<=xx(loop1)
                Min=loop1;
                break;
            end
        end
        for loop1=max(size(xx)):-1:1
            if x(loop+2)>=xx(loop1)
                Max=loop1;
                break;
            end
        end
       
        ytemp=polyval(p,xx(Min:Max));
        yy(Min:Max)=ytemp;
       
       
        if m_drawoption==1
            set(gcf,'DoubleBuffer','on');
            plot(x,y,'-ob',x(loop-1:loop+1),y(loop-1:loop+1),'-*r');
            drawnow;
            pause(0.1);
        end
    end
   
    Break=0;
    for loop=1:1:max(size(x))
       
        if yold(loop)~=0
            Error=abs((y(loop)-yold(loop))/yold(loop)*100);
        elseif y(loop)~=0
            Error=abs((yold(loop)-y(loop))/y(loop)*100);
        else
            Error=0;
        end
       
        if Error>m_Error*100
            Break=1;
        end
    end
 
end