fsolvenew.m (17595B)
1 function [x,FVAL,EXITFLAG,OUTPUT,JACOB] = fsolvenew(FUN,x,options,varargin) 2 %FSOLVE Solves nonlinear equations by a least squares method. 3 % 4 % FSOLVE solves equations of the form: 5 % 6 % F(X)=0 where F and X may be vectors or matrices. 7 % 8 % X=FSOLVE(FUN,X0) starts at the matrix X0 and tries to solve the 9 % equations in FUN. FUN accepts input X and returns a vector (matrix) of 10 % equation values F evaluated at X. 11 % 12 % X=FSOLVE(FUN,X0,OPTIONS) minimizes with the default optimization 13 % parameters replaced by values in the structure OPTIONS, an argument 14 % created with the OPTIMSET function. See OPTIMSET for details. Used 15 % options are Display, TolX, TolFun, DerivativeCheck, Diagnostics, Jacobian, 16 % JacobMult, JacobPattern, LineSearchType, LevenbergMarquardt, MaxFunEvals, 17 % MaxIter, DiffMinChange and DiffMaxChange, LargeScale, MaxPCGIter, 18 % PrecondBandWidth, TolPCG, TypicalX. Use the Jacobian option to specify that 19 % FUN also returns a second output argument J that is the Jacobian matrix at 20 % the point X. If FUN returns a vector F of m components when X has length n, 21 % then J is an m-by-n matrix where J(i,j) is the partial derivative of F(i) 22 % with respect to x(j). (Note that the Jacobian J is the transpose of the 23 % gradient of F.) 24 % 25 % X=FSOLVE(FUN,X0,OPTIONS,P1,P2,...) passes the problem-dependent 26 % parameters P1,P2,... directly to the function FUN: FUN(X,P1,P2,...). 27 % Pass an empty matrix for OPTIONS to use the default values. 28 % 29 % [X,FVAL]=FSOLVE(FUN,X0,...) returns the value of the objective function 30 % at X. 31 % 32 % [X,FVAL,EXITFLAG]=FSOLVE(FUN,X0,...) returns a string EXITFLAG that 33 % describes the exit condition of FSOLVE. 34 % If EXITFLAG is: 35 % > 0 then FSOLVE converged to a solution X. 36 % 0 then the maximum number of function evaluations was reached. 37 % < 0 then FSOLVE did not converge to a solution. 38 % 39 % [X,FVAL,EXITFLAG,OUTPUT]=FSOLVE(FUN,X0,...) returns a structure OUTPUT 40 % with the number of iterations taken in OUTPUT.iterations, the number of 41 % function evaluations in OUTPUT.funcCount, the algorithm used in OUTPUT.algorithm, 42 % the number of CG iterations (if used) in OUTPUT.cgiterations, and the first-order 43 % optimality (if used) in OUTPUT.firstorderopt. 44 % 45 % [X,FVAL,EXITFLAG,OUTPUT,JACOB]=FSOLVE(FUN,X0,...) returns the 46 % Jacobian of FUN at X. 47 % 48 % Examples 49 % FUN can be specified using @: 50 % x = fsolve(@myfun,[2 3 4],optimset('Display','iter')) 51 % 52 % where MYFUN is a MATLAB function such as: 53 % 54 % function F = myfun(x) 55 % F = sin(x); 56 % 57 % FUN can also be an inline object: 58 % 59 % fun = inline('sin(3*x)'); 60 % x = fsolve(fun,[1 4],optimset('Display','off')); 61 % 62 % See also OPTIMSET, LSQNONLIN, @, INLINE. 63 64 % Copyright 1990-2002 The MathWorks, Inc. 65 % $Revision: 1.39 $ $Date: 2002/05/14 19:21:40 $ 66 % Andy Grace 7-9-90. 67 68 % Grandfathered FSOLVE call for Optimization Toolbox versions prior to 2.0: 69 % [X,OPTIONS]=FSOLVE(FUN,X0,OPTIONS,GRADFUN,P1,P2,...) 70 % 71 % ------------Initialization---------------- 72 73 defaultopt = struct('Display','final','LargeScale','off',... 74 'NonlEqnAlgorithm','dogleg',... 75 'TolX',1e-6,'TolFun',1e-6,'DerivativeCheck','off',... 76 'Diagnostics','off',... 77 'Jacobian','off','JacobMult',[],...% JacobMult set to [] by default 78 'JacobPattern','sparse(ones(Jrows,Jcols))',... 79 'MaxFunEvals','100*numberOfVariables',... 80 'DiffMaxChange',1e-1,'DiffMinChange',1e-8,... 81 'PrecondBandWidth',0,'TypicalX','ones(numberOfVariables,1)',... 82 'MaxPCGIter','max(1,floor(numberOfVariables/2))', ... 83 'TolPCG',0.1,'MaxIter',400,... 84 'LineSearchType','quadcubic','LevenbergMarquardt','off'); 85 % If just 'defaults' passed in, return the default options in X 86 if nargin==1 & nargout <= 1 & isequal(FUN,'defaults') 87 x = defaultopt; 88 return 89 end 90 91 if nargin < 2, error('FSOLVE requires two input arguments');end 92 if nargin < 3, options=[]; end 93 94 % These are added so that we can have the same code as in lsqnonlin which 95 % actually has upper and lower bounds. 96 LB = []; UB = []; 97 98 %[x,FVAL,EXITFLAG,OUTPUT,JACOB] = fsolve(FUNin,x,options,varargin) 99 % Note: don't send varargin in as a comma separated list!! 100 numargin = nargin; numargout = nargout; 101 [calltype, GRADFUN, varargin] = parse_call(FUN,options,numargin,numargout,varargin); 102 103 if isequal(calltype,'new') % fsolve version 2.* 104 105 xstart=x(:); 106 numberOfVariables=length(xstart); 107 108 large = 'large-scale'; 109 medium = 'medium-scale: line search'; 110 dogleg = 'trust-region dogleg'; 111 112 switch optimget(options,'Display',defaultopt,'fast') 113 case {'off','none'} 114 verbosity = 0; 115 case 'iter' 116 verbosity = 2; 117 case 'final' 118 verbosity = 1; 119 case 'testing' 120 verbosity = Inf; 121 otherwise 122 verbosity = 1; 123 end 124 diagnostics = isequal(optimget(options,'Diagnostics',defaultopt,'fast'),'on'); 125 gradflag = strcmp(optimget(options,'Jacobian',defaultopt,'fast'),'on'); 126 % 0 means large-scale trust-region, 1 means medium-scale algorithm 127 mediumflag = strcmp(optimget(options,'LargeScale',defaultopt,'fast'),'off'); 128 switch optimget(options,'NonlEqnAlgorithm',defaultopt,'fast') 129 case 'dogleg' 130 algorithmflag = 1; 131 case 'lm' 132 algorithmflag = 2; 133 case 'gn' 134 algorithmflag = 3; 135 otherwise 136 algorithmflag = 1; 137 end 138 mtxmpy = optimget(options,'JacobMult',defaultopt,'fast'); 139 if isequal(mtxmpy,'atamult') 140 warnstr = sprintf('%s\n%s\n%s\n', ... 141 'Potential function name clash with a Toolbox helper function:',... 142 'Use a name besides ''atamult'' for your JacobMult function to',... 143 'avoid errors or unexpected results.'); 144 warning(warnstr) 145 end 146 147 % Convert to inline function as needed 148 if ~isempty(FUN) % will detect empty string, empty matrix, empty cell array 149 [funfcn, msg] = fprefcnchk(FUN,'fsolve',length(varargin),gradflag); 150 else 151 errmsg = sprintf('%s\n%s', ... 152 'FUN must be a function name, valid string expression, or inline object;', ... 153 ' or, FUN may be a cell array that contains these type of objects.'); 154 error(errmsg) 155 end 156 157 JAC = []; 158 x(:) = xstart; 159 switch funfcn{1} 160 case 'fun' 161 fuser = feval(funfcn{3},x,varargin{:}); 162 f = fuser(:); 163 nfun=length(f); 164 case 'fungrad' 165 [fuser,JAC] = feval(funfcn{3},x,varargin{:}); 166 f = fuser(:); 167 nfun=length(f); 168 case 'fun_then_grad' 169 fuser = feval(funfcn{3},x,varargin{:}); 170 f = fuser(:); 171 JAC = feval(funfcn{4},x,varargin{:}); 172 nfun=length(f); 173 otherwise 174 error('Undefined calltype in FSOLVE'); 175 end 176 177 if gradflag 178 % check size of JAC 179 [Jrows, Jcols]=size(JAC); 180 if isempty(mtxmpy) 181 % Not using 'JacobMult' so Jacobian must be correct size 182 if Jrows~=nfun | Jcols ~=numberOfVariables 183 errstr = sprintf('%s\n%s%d%s%d\n',... 184 'User-defined Jacobian is not the correct size:',... 185 ' the Jacobian matrix should be ',nfun,'-by-',numberOfVariables); 186 error(errstr); 187 end 188 end 189 else 190 Jrows = nfun; 191 Jcols = numberOfVariables; 192 end 193 194 YDATA = []; caller = 'fsolve'; 195 196 % large-scale method and enough equations (as many as variables) 197 if ~mediumflag & nfun >= numberOfVariables 198 OUTPUT.algorithm = large; 199 200 % large-scale method and not enough equations -- 201 % switch to medium-scale algorithm 202 elseif ~mediumflag & nfun < numberOfVariables 203 warnstr = sprintf('%s\n%s\n', ... 204 'Large-scale method requires at least as many equations as variables; ',... 205 ' switching to line-search method instead.'); 206 warning(warnstr); 207 OUTPUT.algorithm = medium; 208 209 % medium-scale and no bounds 210 elseif mediumflag & isempty(LB) & isempty(UB) 211 if algorithmflag == 1 & nfun == numberOfVariables % dogleg method 212 OUTPUT.algorithm = dogleg; 213 else 214 if algorithmflag == 1 & nfun ~= numberOfVariables 215 warning('Optimization:fsolve:NonSquareSystem', ... 216 ['Default trust-region dogleg method of FSOLVE cannot\n handle non-square systems; ', ... 217 'switching to Gauss-Newton method.']); 218 algorithmflag = 3; 219 end 220 OUTPUT.algorithm = medium; 221 if algorithmflag == 2 222 % Calling nlsq which looks at LevenbergMarquardt option 223 % (changing option "unsafely" for speed; users should use optimset) 224 options.LevenbergMarquardt = 'on'; 225 else % algorithmflag == 3 226 % (changing option "unsafely" for speed; users should use optimset) 227 options.LevenbergMarquardt = 'off'; 228 end 229 end 230 231 % medium-scale and bounds and enough equations, switch to trust region 232 elseif mediumflag & (~isempty(LB) | ~isempty(UB)) & nfun >= numberOfVariables 233 warnstr = sprintf('%s\n%s\n', ... 234 'Line-search method does not handle bound constraints; ',... 235 ' switching to large-scale trust-region method instead.'); 236 warning(warnstr); 237 OUTPUT.algorithm = large; 238 239 % can't handle this one: 240 elseif mediumflag & (~isempty(LB) | ~isempty(UB)) & nfun < numberOfVariables 241 errstr = sprintf('%s\n%s\n%s\n', ... 242 'Line-search method does not handle bound constraints ',... 243 ' and trust-region method requires at least as many equations as variables; ',... 244 ' aborting.'); 245 error(errstr); 246 247 end 248 249 if diagnostics > 0 250 % Do diagnostics on information so far 251 constflag = 0; gradconstflag = 0; non_eq=0;non_ineq=0;lin_eq=0;lin_ineq=0; 252 confcn{1}=[];c=[];ceq=[];cGRAD=[];ceqGRAD=[]; 253 hessflag = 0; HESS=[]; 254 msg = diagnose('fsolve',OUTPUT,gradflag,hessflag,constflag,gradconstflag,... 255 mediumflag,options,defaultopt,xstart,non_eq,... 256 non_ineq,lin_eq,lin_ineq,LB,UB,funfcn,confcn,f,JAC,HESS,c,ceq,cGRAD,ceqGRAD); 257 258 end 259 260 % Execute algorithm 261 if isequal(OUTPUT.algorithm, large) 262 if ~gradflag 263 Jstr = optimget(options,'JacobPattern',defaultopt,'fast'); 264 if ischar(Jstr) 265 if isequal(lower(Jstr),'sparse(ones(jrows,jcols))') 266 Jstr = sparse(ones(Jrows,Jcols)); 267 else 268 error('Option ''JacobPattern'' must be a matrix if not the default.') 269 end 270 end 271 else 272 Jstr = []; 273 end 274 computeLambda = 0; 275 [x,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT,msg]=... 276 snls(funfcn,x,LB,UB,verbosity,options,defaultopt,f,JAC,YDATA,caller,... 277 Jstr,computeLambda,varargin{:}); 278 % elseif isequal(OUTPUT.algorithm, dogleg) 279 % % trust region dogleg method 280 % Jstr = []; 281 % [x,FVAL,JACOB,EXITFLAG,OUTPUT,msg]=... 282 % trustnleqn(funfcn,x,verbosity,gradflag,options,defaultopt,f,JAC,... 283 % YDATA,Jstr,varargin{:}); 284 else 285 % line search (Gauss-Newton or Levenberg-Marquardt) 286 [x,FVAL,JACOB,EXITFLAG,OUTPUT,msg] = ... 287 nlsq(funfcn,x,verbosity,options,defaultopt,f,JAC,YDATA,caller,varargin{:}); 288 end 289 290 Resnorm = FVAL'*FVAL; % assumes FVAL still a vector 291 if EXITFLAG > 0 % if we think we converged: 292 if Resnorm > sqrt(optimget(options,'TolFun',defaultopt,'fast')) 293 if verbosity > 0 294 disp('Optimizer appears to be converging to a minimum that is not a root:') 295 disp('Sum of squares of the function values is > sqrt(options.TolFun).') 296 disp('Try again with a new starting point.') 297 end 298 EXITFLAG = -1; 299 else 300 if verbosity > 0 301 disp(msg); 302 end 303 end 304 else 305 if verbosity > 0 306 disp(msg); 307 end 308 end 309 310 % Reset FVAL to shape of the user-function output, fuser 311 FVAL = reshape(FVAL,size(fuser)); 312 313 % end FSOLVE 2.* 314 else % version 1.5 FSOLVE 315 316 if length(options)<5; 317 options(5)=0; 318 end 319 % Switch methods making Gauss Newton the default method. 320 if options(5)==0; options(5)=1; else options(5)=0; end 321 322 % Convert to inline function as needed. 323 if ~isempty(FUN) 324 [funfcn, msg] = fcnchk(FUN,length(varargin)); 325 if ~isempty(msg) 326 error(msg); 327 end 328 else 329 error('FUN must be a function name or valid expression.') 330 end 331 332 if ~isempty(GRADFUN) 333 [gradfcn, msg] = fcnchk(GRADFUN,length(varargin)); 334 if ~isempty(msg) 335 error(msg); 336 end 337 else 338 gradfcn = []; 339 end 340 341 [x,options] = nlsqold(funfcn,x,options,gradfcn,varargin{:}); 342 343 if options(8)>10*options(3) & options(1)>0 344 disp('Optimizer is stuck at a minimum that is not a root') 345 disp('Try again with a new starting guess') 346 end 347 348 % Set the second output argument FVAL to be options as in old calling syntax 349 FVAL = options; 350 351 % end fsolve version 1.5.* 352 353 end 354 355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 356 function [calltype, GRADFUN, otherargs] = parse_call(FUN,options,numargin,numargout,otherargs) 357 % PARSE_CALL Determine which calling syntax is being used: the FSOLVE prior to 2.0, or 358 % in version 2.0 or later of the Toolbox. 359 % old call: [X,OPTIONS]=FSOLVE(FUN,X0,OPTIONS,GRADFUN,varargin) 360 % new call: [X,FVAL,EXITFLAG,OUTPUT,JACOB]=FSOLVE(FUN,X0,OPTIONS,varargin) 361 362 if numargout > 2 % [X,FVAL,EXITFLAG,...]=FSOLVE (...) 363 calltype = 'new'; 364 GRADFUN = []; 365 elseif isa(FUN,'cell') % FUN == {...} 366 calltype = 'new'; 367 GRADFUN = []; 368 elseif ~isempty(options) & isa(options,'double') % OPTIONS == scalar or and array 369 calltype = 'old'; 370 if length(otherargs) > 0 371 GRADFUN = otherargs{1}; 372 otherargs = otherargs(2:end); 373 else 374 GRADFUN = []; 375 end 376 elseif isa(options,'struct') % OPTIONS has fields 377 calltype = 'new'; 378 GRADFUN = []; 379 else % Ambiguous 380 warnstr = sprintf('%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n',... 381 'Cannot determine from calling sequence whether to use new',... 382 ' (2.0 or later) FSOLVE function or grandfathered FSOLVE function.',... 383 ' Assuming new syntax; if call was grandfathered FSOLVE syntax, ',... 384 ' this may give unexpected results.',... 385 ' To force new syntax and avoid warning, add options structure argument:',... 386 ' x = fsolve(@sin,3,optimset(''fsolve''));',... 387 ' To force grandfathered syntax and avoid warning, add options array argument:',... 388 ' x = fsolve(@sin,3,foptions);'); 389 warning(warnstr) 390 calltype = 'new'; 391 GRADFUN = []; 392 end 393 394 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 395 function [allfcns,msg] = fprefcnchk(funstr,caller,lenVarIn,gradflag) 396 %PREFCNCHK Pre- and post-process function expression for FUNCHK. 397 % [ALLFCNS,MSG] = PREFUNCHK(FUNSTR,CALLER,lenVarIn,GRADFLAG) takes 398 % the (nonempty) expression FUNSTR from CALLER with LenVarIn extra arguments, 399 % parses it according to what CALLER is, then returns a string or inline 400 % object in ALLFCNS. If an error occurs, this message is put in MSG. 401 % 402 % ALLFCNS is a cell array: 403 % ALLFCNS{1} contains a flag 404 % that says if the objective and gradients are together in one function 405 % (calltype=='fungrad') or in two functions (calltype='fun_then_grad') 406 % or there is no gradient (calltype=='fun'), etc. 407 % ALLFCNS{2} contains the string CALLER. 408 % ALLFCNS{3} contains the objective function 409 % ALLFCNS{4} contains the gradient function (transpose of Jacobian). 410 % 411 % NOTE: we assume FUNSTR is nonempty. 412 % Initialize 413 msg=''; 414 allfcns = {}; 415 funfcn = []; 416 gradfcn = []; 417 418 if gradflag 419 calltype = 'fungrad'; 420 else 421 calltype = 'fun'; 422 end 423 424 % {fun} 425 if isa(funstr, 'cell') & length(funstr)==1 426 % take the cellarray apart: we know it is nonempty 427 if gradflag 428 calltype = 'fungrad';0 429 end 430 [funfcn, msg] = fcnchk(funstr{1},lenVarIn); 431 if ~isempty(msg) 432 error(msg); 433 end 434 435 % {fun,[]} 436 elseif isa(funstr, 'cell') & length(funstr)==2 & isempty(funstr{2}) 437 if gradflag 438 calltype = 'fungrad'; 439 end 440 [funfcn, msg] = fcnchk(funstr{1},lenVarIn); 441 if ~isempty(msg) 442 error(msg); 443 end 444 445 % {fun, grad} 446 elseif isa(funstr, 'cell') & length(funstr)==2 % and ~isempty(funstr{2}) 447 448 [funfcn, msg] = fcnchk(funstr{1},lenVarIn); 449 if ~isempty(msg) 450 error(msg); 451 end 452 [gradfcn, msg] = fcnchk(funstr{2},lenVarIn); 453 if ~isempty(msg) 454 error(msg); 455 end 456 calltype = 'fun_then_grad'; 457 if ~gradflag 458 warnstr = ... 459 sprintf('%s\n%s\n%s\n','Jacobian function provided but OPTIONS.Jacobian=''off'';', ... 460 ' ignoring Jacobian function and using finite-differencing.', ... 461 ' Rerun with OPTIONS.Jacobian=''on'' to use Jacobian function.'); 462 warning(warnstr); 463 calltype = 'fun'; 464 end 465 466 elseif ~isa(funstr, 'cell') %Not a cell; is a string expression, function name string or inline object 467 [funfcn, msg] = fcnchk(funstr,lenVarIn); 468 if ~isempty(msg) 469 error(msg); 470 end 471 if gradflag % gradient and function in one function/M-file 472 gradfcn = funfcn; % Do this so graderr will print the correct name 473 end 474 else 475 errmsg = sprintf('%s\n%s', ... 476 'FUN must be a function name, valid string expression, or inline object;', ... 477 ' or, FUN may be a cell array that contains these type of objects.'); 478 error(errmsg) 479 end 480 481 allfcns{1} = calltype; 482 allfcns{2} = caller; 483 allfcns{3} = funfcn; 484 allfcns{4} = gradfcn; 485 allfcns{5}=[]; 486