ERT电阻层析成像

这套代码很好的对ERT电阻层析成像进行了方针，它可以跟comsol实现通信 将comsol的二维图像保存成.m函数的形式

function [XrhoetaF] = cgls（Abkreorths） 
%CGLS Conjugate gradient algorithm applied implicitly to the normal equations. 
% 
% [XrhoetaF] = cgls（Abkreorths） 
% 
% Performs k steps of the conjugate gradient algorithm applied 
% implicitly to the normal equations A‘*A*x = A‘*b. 
% 
% The routine returns all k solutions stored as columns of 
% the matrix X.  The corresponding solution and residual norms 
% are returned in the vectors eta and rho respectively. 
% 
% If the singular values s are also provided cgls computes the 
% filter factors associated with each step and stores them 
% columnwise in the matrix F. 
% 
% Reorthogonalization of the normal equation residual vectors 
% A‘*（A*X（:i）-b） is controlled by means of reorth: 
%    reorth = 0 : no reorthogonalization （default） 
%    reorth = 1 : reorthogonalization by means of MGS. 
 
% References: A. Bjorck “Numerical Methods for Least Squares Problems“ 
% SIAM Philadelphia 1996. 
% C. R. Vogel “Solving ill-conditioned linear systems using the 
% conjugate gradient method“ Report Dept. of Mathematical 
% Sciences Montana State University 1987. 
  
% Per Christian Hansen IMM April 8 2001. 
 
% The fudge threshold is used to prevent filter factors from exploding. 
fudge_thr = 1e-4; 
  
% Initialization. 
if （k < 1） error（‘Number of steps k must be positive‘） end 
if （nargin==3） reorth = 0; end 
if （nargout==4 & nargin<5） error（‘Too few input arguments‘） end 
if （reorth<0 | reorth>1） error（‘Illegal reorth‘） end 
[mn] = size（A）; X = zeros（nk）; 
if （reorth==1） ATr = zeros（nk）; end 
if （nargout > 1） 
  eta = zeros（k1）; rho = eta; 
end 
if （nargin==5） 
  F = zeros（nk）; Fd = zeros（n1）; s2 = s.^2; 
end 
 
% Prepare for CG iteration. 
x = zeros（n1）; 
d = （b‘*A）‘;   % A‘*b; 
r = b; 
normr2 = d‘*d; 
if （reorth==1） ATr（:1） = d/norm（d）; end 
 
% Iterate. 
for j=1:k 
 
  % Update x and r vectors. 
  Ad = A*d; alpha = normr2/（Ad‘*Ad）; 
  x  = x + alpha*d; 
  r  = r - alpha*Ad; 
  s  = （r‘*A）‘;   % A‘*r; 
 
  % Reorthogonalize s to previous s-vectors if required. 
  if （reorth==1） 
    for i=1:j-1 s = s - （ATr（:i）‘*s）*ATr（:i）; end 
    ATr（:j） = s/norm（s）; 
  end 
 
  % Update d vector. 
  normr2_new = s‘*s; 
  beta = normr2_new/normr2; 
  normr2 = normr2_new; 
  d = s + beta*d; 
  X（:j） = x; 
   
  % Compute norms if required. 
  if （nargout>1） rho（j） = norm（r）; end 
  if （nargout>2） eta（j） = norm（x）; end 
 
  % Compute filter factors if required. 
  if （nargin==5） 
    if （j==1） 
      F（:1） = alpha*s2; 
      Fd = s2 - s2.*F（:1） + beta*s2; 
    else 
      F（:j） = F（:j-1） + alpha*Fd; 
      Fd = s2 - s2.*F（:j） + beta*Fd; 
    end 
    if （j > 2） 
      f = find（abs（F（:j-1）-1） < fudge_thr & abs（F（:j-2）-1） < fudge_thr）; 
      if （length（f） > 0） F（fj） = ones（length（f）1）; end 
    end 
  end 
 
end 

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2851  2002-08-26 16:03  ERT\cgls.m

     文件     859901  2013-08-30 20:35  ERT\cirs.mat

     文件       7445  2013-12-25 20:38  ERT\Currenteit.m

     文件      21628  2013-08-30 20:07  ERT\empty-model.mphbin

     文件     997626  2012-09-02 19:01  ERT\ert仿真操作过程_comsol.pdf

     文件        361  2014-03-01 16:01  ERT\inversecgls.m

     文件       1719  2013-12-25 21:22  ERT\JacobianERT.m

     文件       5311  2013-12-25 20:53  ERT\measure.m

     文件       3528  2013-08-30 20:56  ERT\measure1.m

     文件        610  2013-12-25 21:35  ERT\nodeeit.m

     文件       1432  2013-08-30 20:35  ERT\num.mat

     文件       5534  2013-09-01 16:08  ERT\text.m

     文件       3939  2013-08-30 19:00  ERT\text1.m

     文件       1719  2013-08-30 20:35  ERT\uel.mat

     文件       1682  2013-08-30 20:35  ERT\uref.mat

     文件        492  2013-08-30 20:36  ERT\xy.mat

     文件        130  2013-08-30 18:55  ERT\新建文本文档.txt

     目录          0  2014-11-12 18:20  ERT

----------- ---------  ---------- -----  ----

              1915908                    18