Function [xf,fval,exitflag,iters]=mmfsolve(varargin)



function [xf,fval,exitflag,iters]=mmfsolve(varargin)

%MMFSOLVE Solve a Set of Nonlinear Equations.

% [Xf,Fval,ExitFlag,Iters]=MMFSOLVE(FUN,Xo,Options) finds a zero of the

% vector function FUN(X) starting from the initial guess Xo. FUN is a

% function handle or a function M-file name that evaluates FUN(X). Use

% anonymous functions or nested functions to pass additional parameters to

% FUN. Options is an optional structure that defines algorithm behavior. If

% Options are not given default options are used.

%

% Options=MMFSOLVE('Name',Value,...) sets values in Options structure

% based on the Name/Value pairs:

%

% Name Values {default} Description

% 'Display' ['on' {'off'}] Display iteration information

% 'Jacobian' {'broyden'} Broyden's Jacobian approximation

% 'Jacobian' 'finite' Finite Difference Jacobian

% 'Jacobian' @JNAME Analytic Jacobian function handle

% 'FunTol' {1e-7} NORM(FUN(X),1) stopping tolerance

% 'MaxIter' {100} Maximum number of iterations

% 'MaxStep' value Maximum step size in X allowed

% if 0, algorithm chooses MaxStep

% 'Scale' ['on' {'off'}] Scale algorithm using Xo

%

% Options=MMFSOLVE('default') returns the default options structure.

%

% Options=MMFSOLVE(Options,'Name',Value,...) updates the Options structure

% with new parameter values.

%

% Output Arguments:

% Xf = final solution approximation

% Fval = FUN(Xf)

% ExitFlag = termination code:

% 1 = normal return, 2 = two steps too small

% 3 = line search failure, 4 = too many iterations

% 5 = five steps too big, 6 = stuck at minimizer

% Iters = iteration count

% Reference: Algorithm D6.1.3 in "Numerical Methods for Unconstrained

% Optimization and Nonlinear Equations," by Dennis and Schnabel,

% Siam's Classics in Applied Mathematics, no. 16, (2)nd printing 1996.

% D.C. Hanselman, University of Maine, Orono, ME 04469

% MasteringMatlab@

% Mastering MATLAB 7

% 2005-10-10

% 2005-12-21 fixed bugs noted by Rajeev Kumar

% 2006-05-03 fixed bug noted by Jason Antenucci

% 2007-09-13 works with x having any shape now, noted by Shi Jin

if nargin==1 && strcmpi('default',varargin(1)) % set default parameters

xf.Display='off';

xf.Jacobian='broyden';

xf.MaxIter=100;

xf.MaxStep=0; % to be determined by the algorithm

xf.Scale='off';

xf.FunTol=1e-7;

return

end

if isstruct(varargin{1}) % check if Options=MMFSOLVE(Options,'Name','Value',...)

if ~isfield(varargin{1},'Jacobian')

error('Options Structure not Correct for MMFSOLVE.')

end

xf=varargin{1};

for i=2:2:nargin-1

name=varargin{i};

if ~ischar(name)

error('Parameter Names Must be Strings.')

end

name=lower(name(isletter(name)));

value=varargin{i+1};

if strncmp(name,'d',1), xf.Display=value;

elseif strncmp(name,'f',1), xf.FunTol=value(1);

elseif strncmp(name,'j',1), xf.Jacobian=value;

elseif strncmp(name,'maxi',4), xf.MaxIter=value(1);

elseif strncmp(name,'maxs',4), xf.MaxStep=value(1);

elseif strncmp(name,'s',1), xf.Scale=value;

else disp(['Unknown Parameter Name --> ' name])

end

end

return

end

if ischar(varargin{1}) % check for Options=MMFSOLVE('Name','Value',...)

Pnames=char('display','funtol','jacobian','maxiter','maxstep','scale');

if strncmpi(varargin{1},Pnames,length(varargin{1}))

xf=mmfsolve('default'); % get default values

xf=mmfsolve(xf,varargin{:});

return

end

end

% MMFSOLVE(FUN,Xo,Options,P1,P2,...)

FUN=varargin{1};

if ~(isvarname(FUN) || isa(FUN,'function_handle'))

error('FUN Must be a Function Handle or M-file Name.')

end

xc=varargin{2}(:);

if nargin>2 && isstruct(varargin{3})

options=varargin{3};

else

options=mmfsolve('default');

end

if ~isstruct(options) || ~isfield(options,'Jacobian')

error('Options Structure not Correct for MMFSOLVE.')

end

% Steps 2 and 3

[options,errmsg]=local_neinck(FUN,xc,options);

xf=xc;

error(errmsg)

% Step 4

iters=0;

% Step 5

[fc,FVc]=local_nefn(FUN,xc,options);

% Steps 6 and 7

consecmax=0;

if max(options.Sf.*abs(FVc))1/sqrt(eps)

H=Jc'*diag(options.Sf);

H=H*H';

Hnorm=options.iSx(1)*abs(H(1,:))*options.iSx;

for i=2:n

Hper=abs(H(:,i))'*options.iSx + abs(H(i,:))*options.iSx;

Hper=options.iSx(i)*Hper;

Hnorm=max(Hnorm,Hper);

end

H=H+sqrt(n*eps)*Hnorm*(diag(options.Sx)^2);

M=local_choldecomp(H);

Sn=-M'\(M\gc);

else

for j=2:n

M(1:j-1,j)=M(1:j-1,j)*options.Sx(j);

end

M2=M2.*options.Sx;

Sn=-options.Sf.*FVc;

Sn=local_qrsolve(M,M1,M2,Sn);

end

%---------------------------------------------------------------------------

function [L,maxadd]=local_choldecomp(H)

minl=0;

maxoffl=sqrt(max(diag(H)));

minl2=sqrt(eps)*maxoffl;

maxadd=0;

n=size(H,1);

L=zeros(n);

for j=1:n

if j==1

L(j,j)=H(j,j);

else

L(j,j)=H(j,j)-L(j,1:j-1)*L(j,1:j-1)';

end

minljj=0;

for i=j+1:n

if (j==1)

L(i,j)=H(j,i);

else

L(i,j)=H(j,i)-L(i,1:j-1)*L(j,1:j-1)';

end

minljj=max(abs(L(i,j)),minljj);

end

minljj=max(minljj/maxoffl,minl);

if (L(j,j)>minljj^2) % Normal Cholesky

L(j,j)=sqrt(L(j,j));

else % Augment H(j,j)

minljj=max(minl2,minljj);

maxadd=max(maxadd,(minljj^2-L(j,j)));

L(j,j)=minljj;

end

for i=j+1:n

L(i,j)=L(i,j)/L(j,j);

end

end

%---------------------------------------------------------------------------

function est=local_condest(M,M2)

n=length(M2);

p=zeros(n,1);

pm=p;

x=p;

est=norm(triu(M)-diag(diag(M))+diag(M2),1);

x(1)=1/M2(1);

p(2:n)=M(1,2:n)*x(1);

for j=2:n

xp=(+1-p(j))/M2(j);

xm=(-1-p(j))/M2(j);

temp=abs(xp);

tempm=abs(xm);

for i=j+1:n

pm(i)=p(i)+M(j,i)*xm;

tempm=tempm+abs(pm(i))/abs(M2(i));

p(i)=p(i)+M(j,i)*xp;

temp=temp+abs(p(i))/abs(M2(i));

end

if temp>=tempm

x(j)=xp;

else

x(j)=xm; p(j+1:n)=pm(j+1:n);

end

end

est=est/norm(x,1);

x=local_rsolve(M,M2,x);

est=est*norm(x,1);

%---------------------------------------------------------------------------

function [M,M1,M2,sing]=local_qrdecomp(M)

n=size(M,1);

M1=zeros(n,1);

M2=M1;

sing=0;

for k=1:n-1

eta=max(abs(M(k:n,k)));

if eta==0

M1(k)=0;

M2(k)=0;

sing=1;

else

M(k:n,k)=M(k:n,k)/eta;

sigma=local_sign(M(k,k))*norm(M(k:n,k));

M(k,k)=M(k,k)+sigma;

M1(k)=sigma*M(k,k);

M2(k)=-eta*sigma;

tau=(M(k:n,k)'*M(k:n,(k+1):n))/M1(k);

M(k:n,(k+1):n)=M(k:n,(k+1):n)-M(k:n,k)*tau;

end

end

M2(n)=M(n,n);

if M(n,n)==0

sing=1;

end

%---------------------------------------------------------------------------

function b=local_qrsolve(M,M1,M2,b)

n=length(M2);

for j=1:n-1

tau=(M(j:n,j)'*b(j:n))/M1(j);

b(j:n)=b(j:n)-tau*M(j:n,j);

end

b=local_rsolve(M,M2,b);

%---------------------------------------------------------------------------

function b=local_rsolve(M,M2,b)

n=length(M2);

b(n)=b(n)/M2(n); % rsolve

for i=n-1:-1:1

b(i)=(b(i)-M(i,i+1:n)*b((i+1):n))/M2(i);

end

%---------------------------------------------------------------------------

function Jc=local_broyunfac(Jc,xc,xplus,FVc,FVplus,options)

s=xplus-xc;

Sx=options.Sx;

denom=norm(Sx.*s)^2;

tempi=FVplus-FVc-Jc*s;

i=find(abs(tempi)1

xplus=xc+lambda*p;

[fplus,FVplus]=local_nefn(FUN,xplus,options,varargin{:});

if fplus0.99*options.MaxStep);

elseif lambda ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download