numerical_jacobian

 NUMERICAL_JACOBIAN Determine Jacobian matrix via finite difference at
 certain variable values.

 J = NUMERICAL_JACOBIAN(FUNC,X) computes the Jacobian matrix of the ODE 
 y'=FUNC(y) (i.e. the partial derivatives of FUNC) at variable value 
 column vector X. FUNC is a function handle and depends only on the 
 variable values. Implemented is the finite difference approach in that
 df/dx_k is approximated by 1/e_k * (FUNC(X+h*e_k/2)-FUNC(X-h*e_k/2)) 
 (centered differences) with sufficiently small stepsize h. The size of 
 the matrix will be length(FUNC(X)) times length(X), and is usually 
 quadratic for ODE systems.
 
 J = NUMERICAL_JACOBIAN(FUNC,X,EPSILON,METHOD) allows to tune the stepsize
 EPSILON (default: 1e-6) and to deviate from the default of centered 
 differences by setting METHOD to 'forward' for forward differences 
 (FUNC(X+EPSILON*e_k)-FUNC(X)) or to 'backward' for backward differences 
 in the derivative computation (FUNC(X)-FUNC(X-EPSILON*e_k)). Note that 
 centered differences are considered to be more exact in general than the 
 other two. However, if variable values are close to or at zero, forward 
 differences could be more appropriate to capture the function behavior. 
 
 Example
   J = numerical_jacobian(@(x)func_POSm4(1,x,klin,knonlin),[1 2 1 1])

 See also: numerical_jacobian_complex()

0001 function[J]=numerical_jacobian(func,x,epsilon,method)
0002 % NUMERICAL_JACOBIAN Determine Jacobian matrix via finite difference at
0003 % certain variable values.
0004 %
0005 % J = NUMERICAL_JACOBIAN(FUNC,X) computes the Jacobian matrix of the ODE
0006 % y'=FUNC(y) (i.e. the partial derivatives of FUNC) at variable value
0007 % column vector X. FUNC is a function handle and depends only on the
0008 % variable values. Implemented is the finite difference approach in that
0009 % df/dx_k is approximated by 1/e_k * (FUNC(X+h*e_k/2)-FUNC(X-h*e_k/2))
0010 % (centered differences) with sufficiently small stepsize h. The size of
0011 % the matrix will be length(FUNC(X)) times length(X), and is usually
0012 % quadratic for ODE systems.
0013 %
0014 % J = NUMERICAL_JACOBIAN(FUNC,X,EPSILON,METHOD) allows to tune the stepsize
0015 % EPSILON (default: 1e-6) and to deviate from the default of centered
0016 % differences by setting METHOD to 'forward' for forward differences
0017 % (FUNC(X+EPSILON*e_k)-FUNC(X)) or to 'backward' for backward differences
0018 % in the derivative computation (FUNC(X)-FUNC(X-EPSILON*e_k)). Note that
0019 % centered differences are considered to be more exact in general than the
0020 % other two. However, if variable values are close to or at zero, forward
0021 % differences could be more appropriate to capture the function behavior.
0022 %
0023 % Example
0024 %   J = numerical_jacobian(@(x)func_POSm4(1,x,klin,knonlin),[1 2 1 1])
0025 %
0026 % See also: numerical_jacobian_complex()
0027 
0028 %set default value for the perturbation
0029 if nargin<3
0030     epsilon= 1e-6;
0031 end
0032 if nargin<4
0033     method='centered';
0034 end
0035 
0036 
0037     %for df/dxi:
0038 for i=1:length(x)
0039     xplus=x;
0040     xminus=x;
0041     if strcmp(method, 'centered')
0042         xplus(i)=xplus(i)+epsilon/2;
0043         xminus(i)=xminus(i)-epsilon/2;
0044     end
0045     if strcmp(method, 'backward')
0046         xminus(i)=xminus(i)-epsilon;
0047     end
0048     if strcmp(method,'forward')
0049         xplus(i)=xplus(i)+epsilon;
0050     end
0051     J_temp=(feval(func,xplus)-feval(func,xminus))/epsilon;
0052     if i==1 %during the first run - initialize J
0053         J=zeros(length(J_temp),length(x));
0054     end
0055     J(:,i)=J_temp; %fill the Jacobian for df/dxi
0056 end
0057