jacobian.m

来自「matlab算法集 matlab算法集」· M 代码 · 共 46 行

function J = jacobian (x,f)
%-----------------------------------------------------------------------
% Usage:       J = jacobian (x,f)
%
% Description: Estimate the Jacobain matrix, J = df(x)/dx, numerically
%              using central differences.
%
% Inputs:      x = n by 1 vector at which J is to be evaluated
%              f = string containing name of user-defined function 
%                  which specifies f(x).  The form of f is:
%
%                  function y = f(x)
%
%                  When f is called, it must take the n by 1 vector x
%                  and compute the m by 1 vector y = f(x).
% Outputs:     J = m by n Jacobian matrix of f(x):
%
%                  J(k,j) = df(k)/x(j) , k = [1:m], j = [1:n] 
%-----------------------------------------------------------------------

% Initialize
                       
   chkvec (x,1,'jacobian');
   chkfun (feval(f,x),2,'jacobian');
   h = sqrt(eps);
   n = length (x);
   m = length (feval(f,x));
   x0 = zeros (n,1);
   x1 = zeros (m,1);
   x2 = zeros (m,1);
   x0 = x;

% Compute J = df/dx

   for j = 1 : n
      x0(j) = x(j) + h;  	
      x1 = feval(f,x0);
      x0(j) = x(j) - h;
      x2 = feval(f,x0);
      x0(j) = x(j);
      for i = 1 : m
         J(i,j) = (x1(i) - x2(i))/(2*h);
      end
   end
%-----------------------------------------------------------------------