NLEQS_NEWTON Nonlinear Equation Solver based on Newton's method.
[X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_NEWTON(FCN, X0, OPT)
[X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_NEWTON(PROBLEM)
A function providing a standardized interface for using Newton's
method to solve the nonlinear equation f(x) = 0, beginning from a
starting point x0.
Calls NLEQS_CORE with a Newton update function.
Inputs:
FCN : handle to function that evaluates the function f(x) to
be solved and its Jacobian, J(x). Calling syntax for this
function is:
f = FCN(x)
[f, J] = FCN(x)
If f and x are n x 1, then J is the n x n matrix of partial
derivatives of f (rows) w.r.t. x (cols).
X0 : starting value, x0, of vector x
OPT : optional options structure with the following fields,
all of which are also optional (default values shown in
parentheses)
verbose (0) - controls level of progress output displayed
0 = no progress output
1 = some progress output
2 = verbose progress output
max_it (10) - maximum number of iterations for Newton's method
tol (1e-8) - tolerance on Inf-norm of f(x)
newton_opt - options struct for Newton's method, with field:
lin_solver ('') - linear solver passed to MPLINSOLVE to solve
Newton update step
[ '' - default, '\' for small probs, 'LU3' larger ]
[ '\' - built-in backslash operator ]
[ 'LU' - explicit default LU decomp & back substitutn]
[ 'LU3' - 3 output arg form of LU, Gilbert-Peierls ]
[ algorithm with approximate minimum degree ]
[ (AMD) reordering ]
[ 'LU4' - 4 output arg form of LU, UMFPACK solver ]
[ (same as 'LU') ]
[ 'LU5' - 5 output arg form of LU, UMFPACK solver, ]
[ w/row scaling ]
[ (see MPLINSOLVE for complete list of all options) ]
PROBLEM : The inputs can alternatively be supplied in a single
PROBLEM struct with fields corresponding to the input arguments
described above: fcn, x0, opt
Outputs (all optional, except X):
X : solution vector x
F : final function value, f(x)
EXITFLAG : exit flag
1 = converged
0 or negative values = solver specific failure codes
OUTPUT : output struct with the following fields:
alg - algorithm code of solver used ('NEWTON')
iterations - number of iterations performed
hist - struct array with trajectories of the following:
normf
message - exit message
JAC : final Jacobian matrix, J(x)
Note the calling syntax is almost identical to that of FSOLVE from
MathWorks' Optimization Toolbox. The function for evaluating the
nonlinear function and Jacobian is identical.
Calling syntax options:
[x, f, exitflag, output, jac] = nleqs_newton(fcn, x0);
[x, f, exitflag, output, jac] = nleqs_newton(fcn, x0, opt);
x = nleqs_newton(problem);
where problem is a struct with fields: fcn, x0, opt
and all fields except 'fcn' and 'x0' are optional
x = nleqs_newton(...);
[x, f] = nleqs_newton(...);
[x, f, exitflag] = nleqs_newton(...);
[x, f, exitflag, output] = nleqs_newton(...);
[x, f, exitflag, output, jac] = nleqs_newton(...);
Example: (problem from https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/)
function [f, J] = f1(x)
f = [ x(1) + x(2) - 1;
-x(1)^2 + x(2) + 5 ];
if nargout > 1
J = [1 1; -2*x(1) 1];
end
problem = struct( ...
'fcn', @(x)f1(x), ...
'x0', [0; 0], ...
'opt', struct('verbose', 2) ...
);
[x, f, exitflag, output, jac] = nleqs_newton(problem);
See also NLEQS_MASTER, NLEQS_CORE.0001 function varargout = nleqs_newton(varargin) 0002 %NLEQS_NEWTON Nonlinear Equation Solver based on Newton's method. 0003 % [X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_NEWTON(FCN, X0, OPT) 0004 % [X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_NEWTON(PROBLEM) 0005 % A function providing a standardized interface for using Newton's 0006 % method to solve the nonlinear equation f(x) = 0, beginning from a 0007 % starting point x0. 0008 % 0009 % Calls NLEQS_CORE with a Newton update function. 0010 % 0011 % Inputs: 0012 % FCN : handle to function that evaluates the function f(x) to 0013 % be solved and its Jacobian, J(x). Calling syntax for this 0014 % function is: 0015 % f = FCN(x) 0016 % [f, J] = FCN(x) 0017 % If f and x are n x 1, then J is the n x n matrix of partial 0018 % derivatives of f (rows) w.r.t. x (cols). 0019 % X0 : starting value, x0, of vector x 0020 % OPT : optional options structure with the following fields, 0021 % all of which are also optional (default values shown in 0022 % parentheses) 0023 % verbose (0) - controls level of progress output displayed 0024 % 0 = no progress output 0025 % 1 = some progress output 0026 % 2 = verbose progress output 0027 % max_it (10) - maximum number of iterations for Newton's method 0028 % tol (1e-8) - tolerance on Inf-norm of f(x) 0029 % newton_opt - options struct for Newton's method, with field: 0030 % lin_solver ('') - linear solver passed to MPLINSOLVE to solve 0031 % Newton update step 0032 % [ '' - default, '\' for small probs, 'LU3' larger ] 0033 % [ '\' - built-in backslash operator ] 0034 % [ 'LU' - explicit default LU decomp & back substitutn] 0035 % [ 'LU3' - 3 output arg form of LU, Gilbert-Peierls ] 0036 % [ algorithm with approximate minimum degree ] 0037 % [ (AMD) reordering ] 0038 % [ 'LU4' - 4 output arg form of LU, UMFPACK solver ] 0039 % [ (same as 'LU') ] 0040 % [ 'LU5' - 5 output arg form of LU, UMFPACK solver, ] 0041 % [ w/row scaling ] 0042 % [ (see MPLINSOLVE for complete list of all options) ] 0043 % PROBLEM : The inputs can alternatively be supplied in a single 0044 % PROBLEM struct with fields corresponding to the input arguments 0045 % described above: fcn, x0, opt 0046 % 0047 % Outputs (all optional, except X): 0048 % X : solution vector x 0049 % F : final function value, f(x) 0050 % EXITFLAG : exit flag 0051 % 1 = converged 0052 % 0 or negative values = solver specific failure codes 0053 % OUTPUT : output struct with the following fields: 0054 % alg - algorithm code of solver used ('NEWTON') 0055 % iterations - number of iterations performed 0056 % hist - struct array with trajectories of the following: 0057 % normf 0058 % message - exit message 0059 % JAC : final Jacobian matrix, J(x) 0060 % 0061 % Note the calling syntax is almost identical to that of FSOLVE from 0062 % MathWorks' Optimization Toolbox. The function for evaluating the 0063 % nonlinear function and Jacobian is identical. 0064 % 0065 % Calling syntax options: 0066 % [x, f, exitflag, output, jac] = nleqs_newton(fcn, x0); 0067 % [x, f, exitflag, output, jac] = nleqs_newton(fcn, x0, opt); 0068 % x = nleqs_newton(problem); 0069 % where problem is a struct with fields: fcn, x0, opt 0070 % and all fields except 'fcn' and 'x0' are optional 0071 % x = nleqs_newton(...); 0072 % [x, f] = nleqs_newton(...); 0073 % [x, f, exitflag] = nleqs_newton(...); 0074 % [x, f, exitflag, output] = nleqs_newton(...); 0075 % [x, f, exitflag, output, jac] = nleqs_newton(...); 0076 % 0077 % Example: (problem from https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/) 0078 % function [f, J] = f1(x) 0079 % f = [ x(1) + x(2) - 1; 0080 % -x(1)^2 + x(2) + 5 ]; 0081 % if nargout > 1 0082 % J = [1 1; -2*x(1) 1]; 0083 % end 0084 % 0085 % problem = struct( ... 0086 % 'fcn', @(x)f1(x), ... 0087 % 'x0', [0; 0], ... 0088 % 'opt', struct('verbose', 2) ... 0089 % ); 0090 % [x, f, exitflag, output, jac] = nleqs_newton(problem); 0091 % 0092 % See also NLEQS_MASTER, NLEQS_CORE. 0093 0094 % MP-Opt-Model 0095 % Copyright (c) 1996-2020, Power Systems Engineering Research Center (PSERC) 0096 % by Ray Zimmerman, PSERC Cornell 0097 % 0098 % This file is part of MP-Opt-Model. 0099 % Covered by the 3-clause BSD License (see LICENSE file for details). 0100 % See https://github.com/MATPOWER/mp-opt-model for more info. 0101 0102 %%----- input argument handling ----- 0103 %% gather inputs 0104 if nargin == 1 && isstruct(varargin{1}) %% problem struct 0105 p = varargin{1}; 0106 fcn = p.fcn; 0107 x0 = p.x0; 0108 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0109 else %% individual args 0110 fcn = varargin{1}; 0111 x0 = varargin{2}; 0112 if nargin < 3 0113 opt = []; 0114 else 0115 opt = varargin{3}; 0116 end 0117 end 0118 0119 %% set default options 0120 if isempty(opt) 0121 opt = struct(); 0122 end 0123 if isfield(opt, 'newton_opt') && isfield(opt.newton_opt, 'lin_solver') 0124 lin_solver = opt.newton_opt.lin_solver; 0125 else 0126 lin_solver = ''; 0127 end 0128 0129 %% attempt to pick fastest linear solver, if not specified 0130 if isempty(lin_solver) 0131 nx = size(x0, 1); %% number of variables 0132 if nx <= 10 || have_feature('octave') 0133 lin_solver = '\'; %% default \ operator 0134 else %% MATLAB and nx > 10 0135 lin_solver = 'LU3'; %% LU decomp with 3 output args, AMD ordering 0136 end 0137 end 0138 0139 sp = struct( ... 0140 'alg', 'NEWTON', ... 0141 'name', 'Newton''s', ... 0142 'default_max_it', 10, ... 0143 'need_jac', 1, ... 0144 'update_fcn', @(x, f, J)newton_update_fcn(x, f, J, lin_solver) ); 0145 0146 [varargout{1:nargout}] = nleqs_core(sp, fcn, x0, opt); 0147 % opt.alg = 'CORE'; 0148 % opt.core_sp = sp; 0149 % [varargout{1:nargout}] = nleqs_master(fcn, x0, opt); 0150 0151 0152 function x = newton_update_fcn(x, f, J, lin_solver) 0153 dx = mplinsolve(J, -f, lin_solver); %% compute update step 0154 x = x + dx; %% update x