home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 4.ddi / OPTIM.DI$ / FSOLVE2.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  7.6 KB  |  283 lines
  1. function [xf,termcode,path] = fsolve2(fvec,x0,details,fparam,jac,scale)
  2. %FSOLVE2 Solution to a system of nonlinear equations.
  3. %    X = FSOLVE2('f',X0) starts at X0 and produces a new vector X which
  4. %    solves for f(x) = 0.  'f' is a string containing the name of the
  5. %    function to be solved, normally an M-file.  For example, to
  6. %    find the x, y, and z that solve the simultaneous equations
  7. %
  8. %         sin(x) + y^2 + log(z) - 7 = 0
  9. %         3*x + 2^y - z^3 + 1 = 0
  10. %         x + y + z - 5 = 0
  11. %
  12. %    Use X = FSOLVE2('xyz',[1 1 1]') where XYZ is the M-file
  13. %
  14. %      function q = xyz(p)
  15. %      x = p(1); y = p(2); z = p(3);
  16. %      q = zeros(3,1);
  17. %      q(1) = sin(x) + y^2 + log(z) - 7;
  18. %      q(2) = 3*x + 2^y - z^3 + 1;
  19. %      q(3) = x + y + z - 5;
  20. %
  21. %    FSOLVE2 can take many other optional parameters; see the M-file
  22. %    for more information.
  23.  
  24. %    Copyright (c) 1988 by the MathWorks, Inc.
  25. %    Coded in MATLAB by Richard T. Behrens, April 1988.
  26. %    Revised 11/27/88 JNL
  27.  
  28. %  [XF,TERMCODE,PATH] = FSOLVE(FVEC,X0,DETAILS,FPARAM,JAC,SCALE)
  29. %
  30. %  This function is for solving systems of nonlinear equations.  For more
  31. %  information see a users guide.
  32. %
  33. %  *** INPUTS ***
  34. %  FVEC     - The name of a function which maps n-vectors into n-vectors.
  35. %  X0       - An initial guess of a solution (starting point for iterations).
  36. %  DETAILS  - (optional) A vector whose elements select various algorithmic
  37. %             options and specify various tolerances.
  38. %  FPARAM   - (optional) A set of parameters (constants) which, if nonempty,
  39. %             is passed on as a second argument to FVEC and JAC.
  40. %  JAC      - (optional) The name of a function which implements the Jacobian
  41. %             matrix (if available) of the function named in FVEC.
  42. %  SCALE    - (optional) 'Typical' values of X (1st column) and F (2nd col.)
  43. %             for scaling purposes (no zeros in SCALE, please).
  44. % *** OUTPUTS ***
  45. %  XF       - The final approximation of the solution.
  46. %  TERMCODE - (optional) Indicates the stopping reason (equals 1 for normal
  47. %             termination).
  48. %  PATH     - (optional) Returns the sequence of iterates.
  49. %
  50. %  --> NOTE:  All vectors must be column-vectors.
  51.  
  52. %  Based on Algorithm D6.1.3:  Part of the modular software system from the
  53. %  appendix of the book "Numerical Methods for Unconstrained Optimization and
  54. %  Nonlinear Equations" by Dennis & Schnabel, 1983.  Please refer to that book
  55. %  if you need more detailed information about the algorithms.
  56.  
  57. % Initialization.
  58. %
  59. if (nargin < 3)
  60.    details = zeros(16,1);
  61. end
  62. i = length(details);
  63. if (i < 16)
  64.    details((i+1):16,1) = zeros(16-i,1);  % For unspecified details.
  65. end
  66. if (nargin < 4)
  67.    details(15) = 0;                      % No parameters to pass.
  68. elseif ~length(fparam)
  69.    details(15) = 0;
  70. else
  71.    details(15) = 1;
  72. end
  73. if (details(15) == 0)
  74.    fparam = [];
  75. end
  76. if (nargin < 5)
  77.    details(4) = 0;                       % No analytic Jacobian.
  78. elseif ~length(jac)
  79.    details(4) = 0;
  80. end
  81. if (details(4) == 0)
  82.    jac = '';
  83. end
  84. if (nargin < 6)
  85.    scale = [];                           % No scaling given.
  86.    if (details(16) == 2)
  87.       details(16) = 1;
  88.    end
  89. end
  90. if (nargout < 3)
  91.    details(14) = 0;                      % No path output.
  92. else
  93.    details(14) = 1;
  94. end
  95. if details(14), path = x0.'; end
  96. if (details(1) == 2), btrack = []; end
  97. nofun = 0;        % Number of function evaluations.
  98.  
  99.  
  100.  
  101. %
  102. % Algorithm step 2.
  103. %
  104. if (details(16) > 0)     % Might need F(x0) for scaling.
  105.    if details(15)
  106.       [fc,FVplus,nofun] = nefn(x0,ones(length(x0),1),fvec,nofun,fparam);
  107.    else
  108.       [fc,FVplus,nofun] = nefn(x0,ones(length(x0),1),fvec,nofun);
  109.    end
  110. else
  111.    FVplus = zeros(length(x0),1);
  112. end
  113. [details,Sx,SF,termcode] = neinck(x0,FVplus,details,scale);
  114.  
  115. %
  116. % Algorithm step 3.
  117. %
  118. if (termcode < 0)
  119.    xf = x0;
  120.    if details(14), path = [path;xf.']; end
  121.    return
  122. end
  123.  
  124. %
  125. % Algorithm step 4.
  126. %
  127. itncount = 0;
  128.  
  129. %
  130. % Algorithm step 5.
  131. %
  132. if details(15)
  133.    [fc,FVplus,nofun] = nefn(x0,SF,fvec,nofun,fparam);
  134. else
  135.    [fc,FVplus,nofun] = nefn(x0,SF,fvec,nofun);
  136. end
  137.  
  138. %
  139. % Algorithm step 6.
  140. %
  141. %[termcode,consecmax] = nestop0(x0,FVplus,SF,details(8));
  142. consecmax = 0;
  143. if (max(SF .* abs(FVplus)) <= 1E-2 * details(8))
  144.    termcode = 1;
  145. else
  146.    termcode = 0;
  147. end
  148.  
  149. %
  150. % Algorithm step 7.
  151. %
  152. if (termcode > 0)
  153.    xf = x0;
  154.    if (details(14)), path = [path;xf.']; end
  155. else
  156.    if details(4)
  157.       if details(15)
  158.          [Jc,tmp] = feval(jac,x0,fparam);
  159.       else
  160.          [Jc,tmp] = feval(jac,x0);
  161.       end           
  162.         nofun = tmp + nofun;
  163.    else
  164.       [Jc,nofun] = nefdjac(fvec,FVplus,x0,Sx,details,nofun,fparam);
  165.    end
  166.    gc = Jc' * (FVplus .* (SF.^2));
  167.    FVc = FVplus;
  168. end
  169.  
  170. %
  171. % Algorithm step 8.
  172. %
  173. xc = x0;
  174.  
  175. %
  176. % Algorithm step 9.
  177. %
  178. restart = 1;
  179. norest = 0;
  180.  
  181. %
  182. % Algorithm step 10 (iteration).
  183. %
  184. if (details(1) > 0), clc, end
  185. while (termcode == 0)
  186.    itncount = itncount + 1;
  187.    if (details(4) | details(3) | (1-details(5)))
  188.       [M,Hc,sN] = nemodel(FVc,Jc,gc,SF,Sx,details(2));
  189.    else
  190.       error('Factored model not implemented.')
  191. %     [] = nemodfac();
  192.    end
  193.    if (details(2) == 1)
  194.       [retcode,xplus,fplus,FVplus,maxtaken,nofun,btrack] = ...
  195.       nelnsrch(xc,fc,fvec,gc,sN,Sx,SF,details,nofun,btrack,fparam);
  196.    elseif (details(2) == 2)
  197.       error('Hookstep not implemented.')
  198. %     [] = nehookst();
  199.    else
  200.       error('Dogleg not implemented.')
  201. %     [] = nedogdrv();
  202.    end
  203.    if ((retcode ~= 1) | (restart) | (details(4)) | (details(3)))
  204.       if (details(4))
  205.          if details(15)
  206.             [Jc,tmp] = feval(jac,xplus,fparam);
  207.          else
  208.             [Jc,tmp] = feval(jac,xplus);
  209.          end
  210.             nofun = tmp + nofun;
  211.       elseif (details(3))
  212.          [Jc,nofun] = nefdjac(fvec,FVplus,xplus,Sx,details,nofun,fparam);
  213.       elseif (details(5))
  214.          error('Factored secant method not implemented.')
  215. %        [] = nebroyf();
  216.       else
  217.          Jc = nebroyuf(Jc,xc,xplus,FVc,FVplus,Sx,details(13));  % Broyden update
  218.       end
  219.       if (details(5))
  220.          error('Gradient calculation for factored method not implemented.')
  221.          % Calculate gc using QR factorization (see book).
  222.       else
  223.          gc = Jc' * (FVplus .* (SF.^2));
  224.       end
  225.       [consecmax,termcode] = nestop(xc,xplus,FVplus,fplus,gc,Sx,SF,retcode, ...
  226.                                 details,itncount,maxtaken,consecmax);
  227.    end
  228.    if (((retcode == 1) | (termcode == 2)) & (1-restart) & ...
  229.                                      (1-details(4)) & (1-details(3)))
  230.       [Jc,nofun] = nefdjac(fvec,FVc,xc,Sx,details,nofun,fparam);
  231.       gc = Jc' * (FVc .* (SF.^2));
  232.       if ((details(2) == 2) | (details(2) == 3))
  233.          details(7) = -1;
  234.       end
  235.       restart = 1;
  236.       norest = norest + 1;
  237.    else
  238.       if (termcode > 0)
  239.          xf = xplus;
  240.          if (details(14)), path = [path;xf.']; end
  241.       else
  242.          restart = 0;
  243.          if (details(14)), path = [path;xplus.']; end
  244.       end
  245.       xc = xplus;
  246.       fc = fplus;
  247.       FVc = FVplus;
  248.       if (details(1) > 0)
  249.          home
  250.          disp('The current iteration is: ')
  251.          xc
  252.       end
  253.    end
  254. end
  255.  
  256. if (details(1) == 2)
  257.    disp('Press CR to see statistics . . .')
  258.    fprintf([' ',7])
  259.    pause
  260.  
  261.    clc
  262.    format compact
  263.    disp('Function: ')
  264.    fvec
  265.    disp('Starting point: ')
  266.    x0.'
  267.    disp('Termination condition: ')
  268.    termcode
  269.    disp('Number of iterations: ')
  270.    itncount
  271.    disp('Total number of function evaluations: ')
  272.    nofun
  273.    disp('Final norm(F(x)): ')
  274.    norm(FVc)
  275.    if ((1-details(3)) & (1-details(4)))
  276.       disp('Number of restarts for secant methods: ')
  277.       norest
  278.    end
  279.    disp('Backtrack information: ')
  280.    btrack
  281.    pause
  282. end
  283.