diffctr.m

来自「本人收集的一些有关matlab的代码程序设计 也不知道改选什么分类」· M 代码 · 共 38 行

function dy = diffctr(y, h, n, key)
%DIFFCTR 中心差分算法求取高阶数值微分。
%
%   dy = diffctr(y, h, n, key)
%
%  y 为给定的数据向量，h 为步长， n 为微分阶次 (<=4), 
%  key 是差分算法编号 (1 或 2), 返回的 dy 为数值微分。

%Designed by Prof D Xue (c) 2000
yx1=[y 0 0 0 0 0]; yx2=[0 y 0 0 0 0]; yx3=[0 0 y 0 0 0];
yx4=[0 0 0 y 0 0]; yx5=[0 0 0 0 y 0]; yx6=[0 0 0 0 0 y];
switch n
case 1
   if key==1,dy=(diff(yx1)+diff(yx2))/(2*h); L0=3;
   else, 
      dy = (-diff(yx1)+7*diff(yx2)+7*diff(yx3)-diff(yx4))/(12*h); L0=4;
   end
case 2
   if key==1, dy=(diff(yx1)-diff(yx2))/(h^2); L0=3; 
   else, 
      dy = (-diff(yx1)+15*diff(yx2)-15*diff(yx3)+diff(yx4))/(12*h^2); L0=4;
   end
case 3
   if key==1, dy=(diff(yx1)-diff(yx2)-diff(yx3)+diff(yx4))/(2*h^3); L0=4;
   else, 
      dy = (-diff(yx1)+7*diff(yx2)-6*diff(yx3)-6*diff(yx4)+7*diff(yx5)-...
         diff(yx6))/(8*h^3); L0=5;
   end
case 4
   if key==1,
      dy=(diff(yx1)-3*diff(yx2)+3*diff(yx3)-diff(yx4))/(h^4); L0=4; 
   else, 
      dy = (-diff(yx1)+11*diff(yx2)-28*diff(yx3)+28*diff(yx4)-...
         11*diff(yx5)+diff(yx6))/(6*h^4);L0=5;
   end
end
L1=L0; if key==2, L1=L1+1; end, dy=dy(L1:end-L0);