www.gusucode.com > funfun工具箱matlab源码程序 > funfun/private/ntrp23tb.m
function [yint,ypint] = ntrp23tb(tint,t,y,tnew,ynew,t2,y2,idxNonNegative) %NTRP23TB Interpolation helper function for ODE23TB. % YINT = NTRP23TB(TINT,T,Y,TNEW,YNEW,T2,Y2) uses data computed in ODE23TB % to approximate the solution at TINT. TINT may be a scalar or a row vector. % [YINT,YPINT] = NTRP23TB(TINT,T,Y,TNEW,YNEW,T2,Y2) returns also the % derivative of the polynomial approximating the solution. % % IDX has indices of solution components that must be non-negative. Negative % YINT(IDX) are replaced with zeros and the derivative YPINT(IDX) is set % to zero. % % See also ODE23TB, DEVAL. % Mark W. Reichelt, Lawrence F. Shampine, and Yanyuan Ma, 7-1-97 % Copyright 1984-2005 The MathWorks, Inc. a1 = (((tint - tnew) .* (tint - t2)) ./ ((t - tnew) .* (t - t2))); a2 = (((tint - t) .* (tint - tnew)) ./ ((t2 - t) .* (t2 - tnew))); a3 = (((tint - t) .* (tint - t2)) ./ ((tnew - t) .* (tnew - t2))); yint = y*a1 + y2*a2 + ynew*a3; ypint = []; if nargout > 1 ap1 = (((tint - tnew) + (tint - t2)) ./ ((t - tnew) .* (t - t2))); ap2 = (((tint - t) + (tint - tnew)) ./ ((t2 - t) .* (t2 - tnew))); ap3 = (((tint - t) + (tint - t2)) ./ ((tnew - t) .* (tnew - t2))); ypint = y*ap1 + y2*ap2 + ynew*ap3; end % Non-negative solution if ~isempty(idxNonNegative) idx = find(yint(idxNonNegative,:)<0); % vectorized if ~isempty(idx) w = yint(idxNonNegative,:); w(idx) = 0; yint(idxNonNegative,:) = w; if nargout > 1 % the derivative w = ypint(idxNonNegative,:); w(idx) = 0; ypint(idxNonNegative,:) = w; end end end