1 ## Copyright (C) 2006-2009 Carlo de Falco
3 ## This file is part of:
4 ## OCS - A Circuit Simulator for Octave
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.
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.
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/>.
19 ## author: Carlo de Falco <cdf _AT_ users.sourceforge.net>
20 ## author: Massimiliano Culpo <culpo@math.uni-wuppertal.de>
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});
27 ## Solve a non-linear system of equations using the Newton-Raphson
28 ## method with damping and return the computed solution vector @var{y}.
30 ## The initial guess for the algorithm is set to @var{y0}.
32 ## The Jacobian and residual at each step are computed via the function
33 ## handles @var{RES} and @var{JAC}.
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.
39 ## The optional parameter @var{verbosity} produce verbose output if non-zero.
41 ## The optional function handle @var{update} may be used to provide
42 ## a faster mean to update Jacobian and residual at runtime.
44 ## @var{numit} is the number of performed iterations while @var{resnrm}
45 ## is a vector containing the residual norm at each step.
47 ## @seealso{nls_stationary,tst_backward_euler,tst_theta_method,tst_daspk,tst_odepkg}
50 function [y,ii,resnrm] = nls_newton_raphson(y0,RES,JAC,tol,maxit,\
54 ## FIXME: add input check!
55 if ((nargin < 5) || (nargin > 7))
56 error("nls_newton_raphson: wrong number of input parameters.");
59 if ~exist("verbosity")
71 res_y = RES(y,uptodate{:});
72 resnrm(1) = norm(res_y,inf);
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{:});
81 resnrm(ii+1) = norm(res_y,inf);
84 while ((resnrm(ii+1)>resnrm(ii))&&(jj<10))
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);
98 fprintf(1,"Converged in %d newton iterations and ",ii);
99 fprintf(1,"%d damping iterations.\n",jjtot);
104 fprintf(1,"Not converged, nrm=%g.\n",resnrm(maxit))