numerical_jacobian_complex

 NUMERICAL_JACOBIAN_COMPLEX Determine Jacobian matrix via finite 
 difference with complex step at certain variable values.

 Derivative approximation using a complex step approach (a small imaginary
 number is added) delivers more exact values than normal finite 
 differences that use real steps (Martins et al., 2003). It is exploited 
 that the partial deriavtive df/dx_k can be approximated by 
 Imag(f(x+h*e_k*i))/h for h close to zero. Thereby, no difference has to 
 be computed avoiding subtractive cancellation errors and thus improving 
 computation accuracy.

 J = NUMERICAL_JACOBIAN_COMPLEX(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 using the complex step approach. FUNC is a function
 handle and depends only on the variable values. The size of the matrix 
 is length(FUNC(X)) times length(X), and is quadratic for ODE systems.
 
 J = NUMERICAL_JACOBIAN(FUNC,X,EPSILON) allows to tune the stepsize
 EPSILON (default: 1e-6).
 
 Example
   klin=[165,0.044,0.27,550,5000,78,4.4,5.1];
   knonlin=[0.3,2];
   J = numerical_jacobian_complex(@(x)func_POSm4(1,x,klin,knonlin),...
       [1 2 1 1])

 See also: numerical_jacobian()

 References: Martins JRRA, Sturdza P, Alonso JJ. The complex-step
 derivative approximation. ACM Trans Math Softw. 2003;29(3):245–62.

0001 function[J]=numerical_jacobian_complex(func,x,epsilon)
0002 % NUMERICAL_JACOBIAN_COMPLEX Determine Jacobian matrix via finite
0003 % difference with complex step at certain variable values.
0004 %
0005 % Derivative approximation using a complex step approach (a small imaginary
0006 % number is added) delivers more exact values than normal finite
0007 % differences that use real steps (Martins et al., 2003). It is exploited
0008 % that the partial deriavtive df/dx_k can be approximated by
0009 % Imag(f(x+h*e_k*i))/h for h close to zero. Thereby, no difference has to
0010 % be computed avoiding subtractive cancellation errors and thus improving
0011 % computation accuracy.
0012 %
0013 % J = NUMERICAL_JACOBIAN_COMPLEX(FUNC,X) computes the Jacobian matrix of
0014 % the ODE y'=FUNC(y) (i.e. the partial derivatives of FUNC) at variable
0015 % value column vector X using the complex step approach. FUNC is a function
0016 % handle and depends only on the variable values. The size of the matrix
0017 % is length(FUNC(X)) times length(X), and is quadratic for ODE systems.
0018 %
0019 % J = NUMERICAL_JACOBIAN(FUNC,X,EPSILON) allows to tune the stepsize
0020 % EPSILON (default: 1e-6).
0021 %
0022 % Example
0023 %   klin=[165,0.044,0.27,550,5000,78,4.4,5.1];
0024 %   knonlin=[0.3,2];
0025 %   J = numerical_jacobian_complex(@(x)func_POSm4(1,x,klin,knonlin),...
0026 %       [1 2 1 1])
0027 %
0028 % See also: numerical_jacobian()
0029 %
0030 % References: Martins JRRA, Sturdza P, Alonso JJ. The complex-step
0031 % derivative approximation. ACM Trans Math Softw. 2003;29(3):245–62.
0032 
0033 
0034 %set default value for the perturbation
0035 if nargin<3
0036     epsilon= 1e-6;
0037 end
0038 
0039     %for df/dxi:
0040 for k=1:length(x)
0041     xplus=x;
0042     xplus(k)=xplus(k)+1i*epsilon;
0043     J_temp=imag(feval(func,xplus))/epsilon;
0044     if k==1 %during the first run - initialize J
0045         J=zeros(length(J_temp),length(x));
0046     end
0047     J(:,k)=J_temp; %fill the Jacobian for df/dxi
0048 end
0049