perturb_jacobian_test

 Perturbation Jacobians
 $Id: perturb_jacobian_test.m 3315 2012-07-01 17:32:54Z bgrychtol $

0001 % Perturbation Jacobians
0002 % $Id: perturb_jacobian_test.m 3315 2012-07-01 17:32:54Z bgrychtol $
0003 
0004 % imdl= mk_common_model('c2c2',16);
0005   imdl= mk_common_model('a3cr',16);
0006   imdl= mk_common_model('n3r2',16);
0007   imdl.fwd_model.nodes = imdl.fwd_model.nodes*.25;
0008   img= calc_jacobian_bkgnd(imdl);
0009 
0010   img.fwd_model = mdl_normalize(img.fwd_model,0);
0011 
0012   img.fwd_model.jacobian=   @np_calc_jacobian;
0013   img.fwd_model.system_mat= @np_calc_system_mat;
0014   img.fwd_model.solve=      @np_fwd_solve;
0015   J_np= calc_jacobian( img );
0016 
0017   img.fwd_model.jacobian=   @jacobian_perturb;
0018   img.fwd_model.system_mat= @np_calc_system_mat;
0019   img.fwd_model.solve=      @np_fwd_solve;
0020   J_np_p= calc_jacobian( img );
0021 
0022   img.fwd_model.jacobian=   @jacobian_adjoint;
0023   img.fwd_model.system_mat= @system_mat_1st_order;
0024   img.fwd_model.solve=      @fwd_solve_1st_order;
0025   J_aa= calc_jacobian( img ); % 2 for bug in my code
0026 
0027   img.fwd_model.jacobian=   @jacobian_perturb;
0028   img.fwd_model.system_mat= @system_mat_1st_order;
0029   img.fwd_model.solve=      @fwd_solve_1st_order;
0030   J_aa_p= calc_jacobian( img ); % 2 for bug in my code
0031 
0032   norm(J_aa - J_aa_p,'fro')/norm(J_aa,'fro')
0033   norm(J_np - J_np_p,'fro')/norm(J_np,'fro')
0034   norm(J_np - J_aa_p,'fro')/norm(J_np,'fro')
0035   norm(J_np - J_aa  ,'fro')/norm(J_np,'fro')
0036 
0037 % Demo model with EIDORS3D
0038   imdl= mk_common_model('n3r2',16);
0039   img= calc_jacobian_bkgnd(imdl);
0040   for i=1:16; imdl.fwd_model.electrode(i).z_contact= .01;end
0041 
0042 % Calculate the normal Jacobian
0043   img.fwd_model.jacobian=   @np_calc_jacobian;
0044   img.fwd_model.system_mat= @np_calc_system_mat;
0045   img.fwd_model.solve=      @np_fwd_solve;
0046   J_np= calc_jacobian( img );
0047 
0048 % Calculate the perturbation Jacobian
0049   vh= fwd_solve(img);
0050 
0051   for el= 1:50:size(img.fwd_model.elems,1);
0052      fprintf('\nelem#%03d: ',el);
0053      for delta= [1e-4,1e-6,1e-8] % delta is perturb
0054         img_i = img;
0055         img_i.elem_data(el)= img_i.elem_data(el) + delta;
0056         vi= fwd_solve(img_i);
0057         J_p = (vi.meas - vh.meas)/delta; % perturb Jacobian
0058         fprintf(' %8.6f', norm(J_p - J_np(:,el))/norm(J_np(:,el)));
0059      end
0060   end
0061