| 查看: 364 | 回复: 1 | ||
| 【奖励】 本帖被评价1次，作者fspdlh增加金币 1 个 | ||
| 当前主题已经存档。 | ||
[资源]
【原创】三点法或五点法数值求导程序
|
||
|
把以下代码拷到m文件里就可以了 function y=NumDiff(X,Y,x,method_fit,point_min_ins,method_diff) % Get the differential coefficient of data base on X and Y % y=NumDiff(X,Y) % y=NumDiff(X,Y,x) % y=NumDiff(X,Y,x,method_fit) % y=NumDiff(X,Y,x,method_fit,point_min_ins) % y=NumDiff(X,Y,x,method_fit,point_min_ins,method_diff) % X Y must be vectors of the same length,raw data % x points output % method_fit the method of fit,'spline' or 'pchip' % point_min_ins numbers of points input to the nearest step of X % method_diff 3 or 5, 3 points method or 5 points method %-------------------------------------------------------------------------- %paraments initialization if nargin<2 error('Not enough paraments!'); end if nargin<3 x=X; end if nargin<4 method_fit='pchip'; end if nargin<5 point_min_ins=5; end point_min_ins=point_min_ins+1; if nargin<6 method_diff=5; elseif method_diff~=3 && method_diff~=5 error('method error'); end [m_X n_X]=size(X); if m_X==1 && n_X>1 X=X'; temp=m_X; m_X=n_X; n_X=temp; elseif m_X>1 && n_X>1 error('X must be a vector'); end [m_Y n_Y]=size(Y); if m_Y~=m_X && n_Y==m_X Y=Y'; elseif m_Y~=m_X && n_Y~=m_X error('X and Y must have the same groups of data!'); end %-------------------------------------------------------------------------- %step N=ceil( ( max(X)-min(X) )/ ... ( min(abs(diff(X))/point_min_ins ) ) )+1; h=(max(X)-min(X))./(N-1); X1=[min(X):h:max(X)]'; Y1=interp1(X,Y,X1,method_fit); %-------------------------------------------------------------------------- %coefficient coef=zeros(N,N); if method_diff==5 coef(1:2,1:5)=[-25 48 -36 16 -3;-3 -10 18 -6 1]; coef(3:N+1:N*(N-4))=1; coef(N+3:N+1:N*(N-3))=-8; coef(3*N+3:N+1:N*(N-1))=8; coef(4*N+3:N+1:N*N)=-1; coef(N-1:N,N-4:N)=[-1 6 -18 10 3;3 -16 36 -48 25 ]; coef=coef./12; elseif method_diff==3 coef(1,1:3)=[-3 4 -1]; coef(2:N+1:N*(N-2))=-1; coef(2*N+2:N+1:N*N)=1; coef(N,N-2:N)=[1 -4 3]; coef=coef./2; end %-------------------------------------------------------------------------- %output dif=coef*Y1./h; y=interp1(X1,dif,x,method_fit); [ Last edited by fspdlh on 2009-4-16 at 16:41 ] |
回复此楼» 猜你喜欢
filecode
已经有3人回复
-
我的国基提前知道中了，可是同事的操作让我实在接受不了，怎么会有这样的人
已经有10人回复
-
静等基金结果
已经有13人回复
-
关于Filecode分析方法
已经有3人回复
-
国自然结果
已经有8人回复
-
FileCode能看出啥？
已经有8人回复
-
应该是93bebmhtak前后十一个字符比较关键
已经有22人回复
-
奇怪，两个人的filecode固定段从头到尾一模一样
已经有4人回复
-
8月时间戳变的，举个手。玩一下，释放压力
已经有13人回复
-
好奇怪的filecode
已经有5人回复
2楼2009-03-20 13:46:34
100