www.gusucode.com > Matlab动力系统和时间序列分析工具箱 > Matlab动力系统和时间序列分析工具箱/lab432/toolbox/get_derivative.m
function [Yr,state]=get_derivative(method,order,X,Y,Xr); % method - now available only 'poly^4' % order- order of derivative % given: Y=f(X) % for order=1 function compute Yr=DY/DX=f'(Xr) % for order=2 function compute Yr=D2Y/DX2=f''(Xr) % for order=0 function compute Yr=f(Xr) using poly^4 interpolation % state-{1-OK, 0-Error} % last modified 4.10.04 state=0; switch method case 'poly^4' Yr=ones(size(Xr)); LR=length(Xr); MaxL=length(X); MinL=1; for i=1:LR % R=find(X(MinL:MaxL)<=Xr(i)); tp=MinL; R=0; while tp<=MaxL if X(tp)<=Xr(i) R=R+1; else break end tp=tp+1; end if R==0 return end Cp=max(R)+MinL-1; if Cp<3 if Cp==2 XJ=X(Cp-1:Cp+3); YJ=Y(Cp-1:Cp+3); else XJ=X(Cp:Cp+4); YJ=Y(Cp:Cp+4); end elseif Cp>MaxL-2 if Cp==MaxL-1 XJ=X(Cp-3:Cp+1); YJ=Y(Cp-3:Cp+1); else XJ=X(Cp-4:Cp); YJ=Y(Cp-4:Cp); end else XJ=X(Cp-2:Cp+2); YJ=Y(Cp-2:Cp+2); end for deg=4:-1:2 % deg=4; [p,ok] = lab432_polyfit(XJ,YJ,deg); if ok break end end switch order case 0 Yr(i)=polyval(p,Xr(i)); case 1 Yr(i)=polyval(polyder(p),Xr(i)); case 2 Yr(i)=polyval(polyder(polyder(p)),Xr(i)); otherwise disp('Invalid order. Order must be 0,1 or 2'); end MinL=Cp; end end state=1; function [p,ok] = lab432_polyfit(x,y,n) %POLYFIT Fit polynomial to data. % POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of % degree N that fits the data, P(X(I))~=Y(I), in a least-squares sense. ok=1; x = x(:); y = y(:); if ~isequal(size(x),size(y)) error('X and Y vectors must be the same size.') end V(:,n+1) = ones(length(x),1); for j = n:-1:1 V(:,j) = x.*V(:,j+1); end % Solve least squares problem [Q,R] = qr(V,0); ws = warning('off','all'); p = R\(Q'*y); % Same as p = V\y; warning(ws); if size(R,2) > size(R,1) ok=0; elseif condest(R) > 1.0e10 ok=0; end p = p.';