%利用三點公式計算數(shù)值（偏）導(dǎo)數(shù)。
%說明：三點公式diff3p()的計算精度比diff()和gradient()相對更高些
%(1):程序在diff3p.m,調(diào)用格式如下
%DF=DIFF3P(P,H)   計算一階數(shù)值導(dǎo)數(shù)
%DF=DIFF3P(P,H,O,D)  計算D指定的變量的O階數(shù)值（偏）導(dǎo)數(shù)
%（2）H是步長，O-order,指定導(dǎo)數(shù)的階數(shù)，D=1表示對X求導(dǎo)，D=2對Y求導(dǎo)，它們的默認(rèn)值均為1。
%F是函數(shù)F(X)或F(X,Y)在離散數(shù)據(jù)點的函數(shù)值數(shù)組，DF與F是同規(guī)模的數(shù)組。
%（3）當(dāng)F是向量時，DF是函數(shù)F(X)在離散數(shù)據(jù)點的數(shù)值導(dǎo)數(shù)。當(dāng)F是矩陣時，若D=1則DF的每一行是函數(shù)F(X)的數(shù)值導(dǎo)數(shù)或F(X,Y)
%對X的數(shù)值偏導(dǎo)數(shù)，D=2則DF的每一列是F(X,Y)對Y的數(shù)值偏導(dǎo)數(shù)
%（4）適合計算一階或二階導(dǎo)數(shù)
 
 
% 三點公式：F(X)在三點x0,x1=x0+h,x2=x0+2h處的一階數(shù)值導(dǎo)數(shù)為
%f'(x0)=(1/2h)[-3f(x0)+4f(x1)-f(x2)]
%f'(x1)=(1/2h)[-f(x0)+f(x2)]
%f'(x2)=(1/2h)[f(x0)-4f(x1)+3f(x2)]
%f(x)在三點x0,x1,x2處的二階數(shù)值導(dǎo)數(shù)為(i=1,2,3)
%f''(xi)=(1/h)[f(x0)-2f(x1)+f(x2)]
function df=diff3p(f,h,o,d)
if nargin<2
    h=1;
end
if nargin<3
    o=1;
end
if nargin<4
    d=1;
end
dim=size(f);
if(dim(2)==1&d~=1)
    strf=setstr(inputname(1));
    strd=setstr(inputname(4));
    fprintf('the dimension of %s doesn''t match with %s \n',strf,strd);
    df=[];
    return;
end
if dim(2)==1
    f=f.';
end
if all(dim>[1,1])
    if d==1
        dd=dim(2);
    else
        dd=dim(1);
    end
elseif dim(1)==1
    dd=dim(2);
else
    dd=dim(2);
end
if d==2
    f=f.';
end
df(:,1)=-3.*f(:,1)+4.*f(:,2)-f(:,3);
df(:,dd)=f(:,dd-2)-4.*f(:,dd-1)+3.*f(:,dd);
for k=2:dd-1
    df(:,k)=f(:,k+1)-f(:,k-1);
end
df=df./(2*h);
if dim(2)==1
    f=f.';
end
if d==2
    df=df.';
    f=f.';
end
if o~=1
    for kk=2:o
        df==diff3p(df,h,o-1,d);
    end
end