inv_sol.sci

来自「用来实现三维阻抗及光学断层成像重建的matlab程序」· SCI 代码 · 共 91 行

function [solf,solp]=inv_sol(I,voltage,tol,mat_ref,vtx,%simp,elec,no_pl,zc,sym,gnd_ind,tfac,Reg,it)
solf=[];solp=[];
//function [solf,solp] = inv_sol(J,I,voltage,tol,mat_ref,vtx,simp,elec,no_pl,zc,sym,gnd_ind,tfac,Reg,it);
// 
//Calculates a Newton non-linear inverse solution.
// 
// 
// 
//solf    = The inverse solution
//mat_ref = Initial guess on the solution
//I       = The current patterns 
//elec    = The electrode faces
//zc      = The contact impedance vector
//voltage = The measurements
//tfac    = The regularisation parameter
//Reg     = The regularisation matrix 
//it      = Number of iterations
//vtx     = The vertices matrix
//simp    = The simplices matrix
//gnd_ind = The ground index (node/electrode)
//no_pl   = The number of planes
 
 
tol = 0.0001;
//Inverse calculations error tolerance. Change accodringly 
 
sol_upd = mat_ref;
//Initial estimate - homogeneous background
 
solp = zeros(size(%simp,1),1);
//Total change (plotting purposes)
 
 
//! mtlb(end) can be replaced by end() or end whether end is an m-file or not
Ib = I(size(vtx,1)+1:mtlb(end),:);
 
el_no = size(elec,1);
 
for i = 1:it
   
  [E,pp] = fem_master_full(vtx,%simp,sol_upd,gnd_ind,elec,zc,sym);
   
  IntGrad = integrofgrad(vtx,%simp,sol_upd);
   
  if i==1 then
    mtlb_sprintf('Current fields for iteration %d',i);
    V = forward_solver(vtx,E,I,tol,pp);
    [viH,viV,indH,indV,df] = get_3d_meas(elec,vtx,V,Ib,no_pl);
     
    //! mtlb(end) can be replaced by end() or end whether end is an m-file or not
    dfv = df(1:2:mtlb(end));
    vi = viH;
    mtlb_sprintf('Measurement fields for iteration %d',i);
    v_f = m_3d_fields(vtx,el_no,indH,E,tol,gnd_ind);
  else
    mtlb_sprintf('Current fields for iteration %d',i);
    V = forward_solver(vtx,E,I,tol,pp,V);
    [viH,viV,indH,indV,df] = get_3d_meas(elec,vtx,V,Ib,no_pl);
     
    //! mtlb(end) can be replaced by end() or end whether end is an m-file or not
    dfv = df(1:2:mtlb(end));
    vi = viH;
    mtlb_sprintf('Measurement fields for iteration %d',i);
    v_f = m_3d_fields(vtx,el_no,indH,E,tol,gnd_ind,v_f);
  end
   
  J = jacobian_3d(I,elec,vtx,%simp,gnd_ind,sol_upd,zc,IntGrad,v_f,dfv,tol,sym);
   
   
  sol = (J'*J+tfac*Reg'*Reg)\(J'*(vi-voltage));
  // - tfac*Reg'*Reg*sol);
  //sol = pinv([J;sgrt(tfac)*Reg],tol) * [sqrt(tfac)*Reg*sol_upd; (ref-voltage)];
  sol_upd = sol_upd+sol;
  solp = solp+sol;
   
  mtlb_sprintf('Error norm at iteration %d is %f',i,norm(vi-voltage));
   
end
//for it
 
solf = mat_ref+sol;
 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// This is part of the EIDORS suite.
// Copyright (c) N. Polydorides 2001
// Copying permitted under terms of GNU GPL
// See enclosed file gpl.html for details.
// EIDORS 3D version 1.0
// MATLAB version 5.3 R11
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%