jacobian copy.m

MDPSAS工具箱是马里兰大学开发的

function [J, R, V, template] = jacobian(A);

% Use finite-differences to compute Jacobian

% Copyright (c) by Raymond A. Adomaitis, 1998-2005

Akey = keys(get(A,'var'));
Aval = values(get(A,'var'));
[V, template] = columnnae(A,'var');

Aresjacflag = get(A,'resjacflag');

% Finite difference scaling is based on Solving Nonlinear Equations
% with Newton's Method, C. T. Kelly, SIAM, 2003.

Epsil      = sqrt(eps)*max([abs(V),ones(size(V))],[],2).*sgn(V);
VPlusEpsil = V + Epsil;

if ~Aresjacflag     % default - residual is not set up to
                    % compute the Jacobian by finite differences

% check if residual in object A must be updated

   if ~get(A,'residflag'), A = residual(A); end
   R = columnnae(A,'resid');

   n = 0;
   B = A;

   for i = 1:length(Akey)
      for j = 1:prod(template{i})

         n = n+1;
         B = set(B,'var', ...
                values(setval(get(A,'var'),VPlusEpsil(n),Akey{i},j)));
         B = residual(B);
            
         RPlusEpsil = columnnae(B,'resid');
         J(:,n) = ( RPlusEpsil - R )/Epsil(n);
      
         end
      end

else % residual is set up to compute Jacobian by finite differences

   [A, J, R] = residual(A,Akey,Aval,V,template, Epsil,VPlusEpsil);
   end

% RAA 9/15/2004 changed from: 
%
%    Epsil = 1e-4*mean(abs(V));
%    if Epsil < sqrt(eps), Epsil = sqrt(eps); end
%    VDiff = Epsil*V;
%    I = find( abs(VDiff) < Epsil );
%    VPlusEpsil(I) = V(I) + Epsil;
%    VDiff(I) = Epsil;
%    J(:,n) = ( RPlusEpsil - R )/VDiff(n);

% -------------------- sgn.m -------------------- %

function B = sgn(A)

% return the sign of A if nonzero; return 1 if zero

I = find( A==0 );
if ~isempty(I)
   B(I,1) = 1;
   end

J = find( A~=0 );
if ~isempty(J)
   B(J,1) = A(J,1)./abs(A(J,1));
   end