octave_packages/ocs-0.1.3/nls/nls_newton_raphson.m

   1 ## Copyright (C) 2006-2009  Carlo de Falco
   2 ##
   3 ## This file is part of:
   4 ## OCS - A Circuit Simulator for Octave
   5 ##
   6 ## OCS is free software; you can redistribute it and/or modify
   7 ## it under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation.
   9 ##
  10 ## This program is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with this program (see the file LICENSE); if not,
  17 ## see <http://www.gnu.org/licenses/>.
  18 ##
  19 ## author: Carlo de Falco <cdf _AT_ users.sourceforge.net>
  20 ## author: Massimiliano Culpo <culpo@math.uni-wuppertal.de>
  21
  22 ## -*- texinfo -*-
  23 ## @deftypefn{Function File}{[@var{y},@var{numit},@var{resnrm}] =} @
  24 ## nls_newton_raphson (@var{y0},@var{RES},@var{JAC},@var{tol},@
  25 ## @var{maxit},@var{verbosity},@var{update});
  26 ##
  27 ## Solve a non-linear system of equations using the Newton-Raphson
  28 ## method with damping and return the computed solution vector @var{y}.
  29 ##
  30 ## The initial guess for the algorithm is set to @var{y0}.
  31 ##
  32 ## The Jacobian and residual at each step are computed via the function
  33 ## handles @var{RES} and @var{JAC}.
  34 ##
  35 ## The variables @var{tol} and @var{maxit} are the relative tolerance on the
  36 ## error of the computed solution and the maximum number of iterations to be
  37 ## performed by the algorithm.
  38 ##
  39 ## The optional parameter @var{verbosity} produce verbose output if non-zero.
  40 ##
  41 ## The optional function handle @var{update} may be used to provide
  42 ## a faster mean to update Jacobian and residual at runtime.
  43 ##
  44 ## @var{numit} is the number of performed iterations while @var{resnrm}
  45 ## is a vector containing the residual norm at each step.
  46 ##
  47 ## @seealso{nls_stationary,tst_backward_euler,tst_theta_method,tst_daspk,tst_odepkg}
  48 ## @end deftypefn
  49
  50 function [y,ii,resnrm] = nls_newton_raphson(y0,RES,JAC,tol,maxit,\
  51                                             verbosity,update);
  52
  53   ## Check input
  54   ## FIXME: add input check!
  55   if ((nargin < 5) || (nargin > 7))
  56     error("nls_newton_raphson: wrong number of input parameters.");
  57   endif
  58
  59   if ~exist("verbosity")
  60     verbosity = 0;
  61   endif
  62
  63   if ~exist("update")
  64     update = @(x) ({});
  65   endif
  66
  67   jjtot = 0;
  68   y     = y0;
  69
  70   uptodate  = update(y);
  71   res_y     = RES(y,uptodate{:});
  72   resnrm(1) = norm(res_y,inf);
  73
  74   for ii = 1:maxit
  75
  76     jac_y      = JAC(y,uptodate{:});
  77     ynew       = jac_y\(-res_y+jac_y*y);
  78     uptodate   = update(ynew);
  79     res_y      = RES(ynew,uptodate{:});
  80
  81     resnrm(ii+1) = norm(res_y,inf);
  82
  83     jj = 0;
  84     while ((resnrm(ii+1)>resnrm(ii))&&(jj<10))
  85       jj++;
  86       damp = 2^(-jj);
  87       ynew = y*(1-damp) + ynew*damp;
  88       uptodate = update(ynew);
  89       res_y    = RES(ynew,uptodate{:});
  90       resnrm(ii+1) = norm(res_y,inf);
  91     endwhile
  92
  93     jjtot += jj;
  94     y      = ynew;
  95
  96     if resnrm(ii+1)<tol
  97       if (verbosity)
  98         fprintf(1,"Converged in %d newton iterations and ",ii);
  99         fprintf(1,"%d damping iterations.\n",jjtot);
 100       endif
 101       break
 102     elseif ii==maxit
 103       if(verbosity)
 104         fprintf(1,"Not converged, nrm=%g.\n",resnrm(maxit))
 105       endif
 106       break
 107     endif
 108   endfor
 109
 110 endfunction