Code covered by the BSD License  

Highlights from
Complex fsolve

Complex fsolve

by

 

Computes numerically complex solutions of a system of equations

[cfsol fsol fval exitflag output jacob ]=cfsolve(funhandle, gradhandle, hessianhandle, InitReal,InitIm,varargin)
function [cfsol fsol fval exitflag output jacob ]=cfsolve(funhandle, gradhandle, hessianhandle, InitReal,InitIm,varargin)
% complex solver quasi-newton solver for matlab
% Calling: [cfsol fsol fval exitflag output jacob ]= cfsolve(funhandle, gradhandle, hessianhandle, InitReal,InitIm,varargin)
%  
% Computes numerically a complex solution of the function in funhandle. funhandle can be an
% anonymous function, a function handle, or a symbolic that computes a set of
% equations. Transforms internally a symbolic function to matlab function. Uses
% persistent function to store these for subsequent calls. If any of the
% functions are changed it will recreate the matlab functions of the symbolics them.
% 
% INPUT: 
%
%       funhandle: a symbolic or anonymous function , can be matrix or vector
%       valued .i.e. f:C^m --> C^{nxn}
%
%       gradhandle: gradhandle is the jacobian, 
%            
%       hessianhandle: the hessian of funhandle
%
%       InitReal: The initialisation of the real part. A vector in R^m where m
%       is the number of variables in funhandle.
%
%       InitIm: The initialisation of the complex part. Setting this to empty
%       will cause the solver to search only for real solutions.
%       
%       options: options are passed as in the case of fsolve
%
% OUTPUT:
%
%       cfsol: The complex solutions in the form cfsol(k)=a(k)+b(k)*j
%
%       fsol:The complex solutions in in parts , i.e. cfsol(1,k)=a(k) ,
%       cfsol(2,k)=b(k)
%
%       fval: Objective function termination value. The objective is
%       sum(funhandle.^2)
%
%       exitflag: exitvalues as in fsolve
%
%       output: as in fsolve, a structure describing the termination reasons
%
%       jacob: The jacobian at the solution
%


persistent funhandle_ funhandle_sym
persistent gradhandle_ gradhandle_sym
persistent hessianhandle_ hessianhandle_sym



sir=size(InitReal);
sim=size(InitIm);

if ~isempty(gradhandle)
  SOLVER_=1;
  if ~isempty(hessianhandle)
    SOLVER_=2;
  end
else SOLVER_=0;
end


if any(sir~=sim) && any(sim~=[ 0 0 ])
  error('cfsolve:Invalidinput','Init dimensions are not compatible');
else
  if sim(1)>sim(2)
    warning('cfsolve:transposing','Transposing InitIm (arg3) to match dimensions of InitReal (arg2)');
    InitIm=InitIm.';
  end
  if sir(1)>sir(2)
    warning('cfsolve:transposing','Trasposing InitReal (arg2) to match dimensions of InitIm');
    InitReal=InitReal.';
  end
end

InitVal=[InitReal; InitIm];
options=optimset(varargin{:});
if SOLVER_
  options=optimset(options,'GradObj','on');
  if SOLVER_==2
    options=optimset(option,'Hessian','on');
  end
end

switch SOLVER_
  case 1
    if isa(funhandle,'sym') % Symbolic handler      
      if isempty(funhandle_sym) || any(size(funhandle)~=size(funhandle_sym)) || any(any(funhandle~=funhandle_sym))
        funhandle_sym=funhandle;
        funhandle=matlabFunction(funhandle);
        funhandle_=funhandle;
      else
        funhandle=funhandle_;
      end      
    end
    if isa(gradhandle,'sym')
      if isempty(gradhandle_sym) || any(size(gradhandle)~=size(gradhandle_sym)) || any(any(gradhandle~=gradhandle_sym))
        gradhandle_sym=gradhandle;
        gradhandle=matlabFunction(gradhandle);
        gradhandle_=gradhandle;
      else
        gradhandle=gradhandle_;
      end            
    end
    [fsol fval exitflag output jacob]=fsolve(@(x) thisfun_G(funhandle,gradhandle,x), InitVal, options);
  case 2
    if isa(funhandle,'sym') % Symbolic handler      
      if isempty(funhandle_sym) || any(size(funhandle)~=size(funhandle_sym)) || any(any(funhandle~=funhandle_sym))
        funhandle_sym=funhandle;
        funhandle=matlabFunction(funhandle);
        funhandle_=funhandle;
      else
        funhandle=funhandle_;
      end      
    end
    if isa(gradhandle,'sym')
      if isempty(gradhandle_sym) || any(size(gradhandle)~=size(gradhandle_sym)) || any(any(gradhandle~=gradhandle_sym))
        gradhandle_sym=gradhandle;
        gradhandle=matlabFunction(gradhandle);
        gradhandle_=gradhandle;
      else
        gradhandle=gradhandle_;
      end            
    end
    if isa(hessianhandle,'sym')
      if isempty(hessianhandle_sym) || any(size(hessianhandle)~=size(hessianhandle_sym)) || any(any(hessianhandle~=hessianhandle_sym))
        hessianhandle_sym=hessianhandle;
        hessianhandle=matlabFunction(hessianhandle);
        hessianhandle_=hessianhandle;
      else
        hessianhandle=hessianhandle_;
      end            
    end
    [fsol fval exitflag output jacob]=fsolve(@(x) thisfun_GH(funhandle,gradhandle,hessianhandle,x), InitVal, options);
  otherwise
    if isa(funhandle,'sym')
      if isempty(funhandle_sym) || any(size(funhandle)~=size(funhandle_sym)) || any(any(funhandle~=funhandle_sym))
        funhandle_sym=funhandle;
        funhandle=matlabFunction(funhandle);
        funhandle_=funhandle;
      else
        funhandle=funhandle_;
      end      
    end
    [fsol fval exitflag output jacob]=fsolve(@(x) thisfun(funhandle,x), InitVal, options);
end

if size(fsol,1)==2
  cfsol=fsol(1,:)+fsol(2,:)*1i;
else
  cfsol=fsol;
end

end

% Internal function block

function fsolsep=thisfun(funhandle,z)

z=realorcomplex_(z);
if nargin(funhandle)==numel(z)
  zc=dealz_(z);
  fval=feval(funhandle,zc{:});
else
  fval=feval(funhandle,z);
end

fsolsep=[real(fval) imag(fval)];

end

function [fval gval]=thisfun_G(funhandle,gradhandle,z)

z=realorcomplex_(z);
if nargin(funhandle)==numel(z)
  zc=dealz_(z);
  fval_=feval(funhandle,zc{:});
  gval_=feval(gradhandle,zc{:});
else
  fval_=feval(funhandle,z);
  gval_=feval(gradhandle,z);
end
fval=[real(fval_) imag(fval_)];
gval=[real(gval_) imag(gval_)];
end

function [fval gval hval]=thisfun_GH(funhandle,gradhandle,hessianhandle,z)

z=realorcomplex_(z);
if nargin(funhandle)==numel(z)
  zc=dealz_(z);
  fval_=feval(funhandle,zc{:});
  gval_=feval(gradhandle,zc{:});
  hval_=feval(hessianhandle,zc{:});
else
  fval_=feval(funhandle,z);
  gval_=feval(gradhandle,z);
  hval_=feval(hessianhandle,z);
end

fval=[real(fval_) imag(fval_)];
gval=[real(gval_) imag(gval_)];
hval=[real(hval_) imag(hval_)];

end

function [zc z]=dealz_(z)

for i=1:size(z,2)
  zc{i}=z(:,i);
end

end

function z=realorcomplex_(z)

if size(z,1)==2
  z=z(1,:)+1i*z(2,:);
end

end

Contact us