www.gusucode.com > funfun工具箱matlab源码程序 > funfun/private/ntrp23.m
function [yinterp,ypinterp] = ntrp23(tinterp,t,y,~,~,h,f,idxNonNegative) %NTRP23 Interpolation helper function for ODE23. % YINTERP = NTRP23(TINTERP,T,Y,TNEW,YNEW,H,F,IDX) uses data computed in ODE23 % to approximate the solution at time TINTERP. TINTERP may be a scalar % or a row vector. % The arguments TNEW and YNEW do not affect the computations. They are % required for consistency of syntax with other interpolation functions. % Any values entered for TNEW and YNEW are ignored. % % [YINTERP,YPINTERP] = NTRP23(TINTERP,T,Y,TNEW,YNEW,H,F,IDX) returns also the % derivative of the polynomial approximating the solution. % % IDX has indices of solution components that must be non-negative. Negative % YINTERP(IDX) are replaced with zeros and the derivative YPINTERP(IDX) is % set to zero. % % See also ODE23, DEVAL. % Mark W. Reichelt and Lawrence F. Shampine, 6-13-94 % Copyright 1984-2009 The MathWorks, Inc. BI = [ 1 -4/3 5/9 0 1 -2/3 0 4/3 -8/9 0 -1 1 ]; s = (tinterp - t)/h; yinterp = y(:,ones(size(tinterp))) + f*(h*BI)*cumprod([s;s;s]); ypinterp = []; if nargout > 1 ypinterp = f*BI*[ ones(size(s)); cumprod([2*s;3/2*s])]; end % Non-negative solution if ~isempty(idxNonNegative) idx = find(yinterp(idxNonNegative,:)<0); % vectorized if ~isempty(idx) w = yinterp(idxNonNegative,:); w(idx) = 0; yinterp(idxNonNegative,:) = w; if nargout > 1 % the derivative w = ypinterp(idxNonNegative,:); w(idx) = 0; ypinterp(idxNonNegative,:) = w; end end end