diff3p.m

数值微分三点公式法

%利用三点公式计算数值（偏）导数。
%说明：三点公式diff3p()的计算精度比diff()和gradient()相对更高些
%(1):程序在diff3p.m,调用格式如下
%DF=DIFF3P(P,H)   计算一阶数值导数
%DF=DIFF3P(P,H,O,D)  计算D指定的变量的O阶数值（偏）导数
%（2）H是步长，O-order,指定导数的阶数，D=1表示对X求导，D=2对Y求导，它们的默认值均为1。
%F是函数F(X)或F(X,Y)在离散数据点的函数值数组，DF与F是同规模的数组。
%（3）当F是向量时，DF是函数F(X)在离散数据点的数值导数。当F是矩阵时，若D=1则DF的每一行是函数F(X)的数值导数或F(X,Y)
%对X的数值偏导数，D=2则DF的每一列是F(X,Y)对Y的数值偏导数
%（4）适合计算一阶或二阶导数
 
 
% 三点公式：F(X)在三点x0,x1=x0+h,x2=x0+2h处的一阶数值导数为
%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处的二阶数值导数为(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