www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregusermod/CalcJacob.m

    function J= CalcJacob(U,x)
%CALCJACOB calculate jacobian
%
% J= CalcJacob(U,x)

%  Copyright 2000-2007 The MathWorks, Inc. and Ford Global Technologies, Inc.

try
    J= feval(U.funcName,U,'jacobian',x);
    if isempty(J)
        % Jacobian is not defined in User function so use a numerical Jacobian
        J = iNumericalJacobian(U,x);
    end
    if ~all(isfinite(J(:)))
        error(message('mbc:xregusermod:InvalidValue'))
    end
catch ME
    if ~strcmp(ME.identifier,'MATLAB:UndefinedFunction') || ~isempty(strfind(ME.message,'i_jacobian'))
        % Jacobian is not defined in User function so use a numerical Jacobian
        J = iNumericalJacobian(U,x);
    else
        ReportError(U,ME)
        % don't make jacobian NaN as sometimes we want to do a QR decomp and
        % this can be bad
        J= zeros(size(x,1),length(U.parameters));
    end
end

function J = iNumericalJacobian(U,x)

try
    FTY= eval(U,x);
    Thresh = iCalcThreshold(U);
    J = numjac(@fbeval,U,double(U),FTY,Thresh,[],0,[],[],x);
    if ~all(isfinite(J(:)))
        error(message('mbc:xregusermod:InvalidValue'))
    end
catch ME
    ReportError(U,ME)
    J= zeros(size(x,1),length(U.parameters));
end

function Thresh = iCalcThreshold(U)

b = double(U);    
Threshold = max(abs(b*1e-8),1e-6);
% set typical value to nonzero where b is 0 otherwise numjac is poor

TypicalY = initial(U);
TypicalY(TypicalY==0) = 1;

Thresh = {Threshold,TypicalY};