NLEQS_CORE Core Nonlinear Equation Solver with arbitrary update function.
[X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_CORE(SP, FCN, X0, OPT)
A function providing the core code for a standardized interface for
various methods of solving the nonlinear equation f(x) = 0, beginning
from a starting point x0. Allows for an arbitrary update function.
Inputs:
SP : solver parameters, struct with the following fields (all required)
alg : algorithm code, e.g. 'NEWTON', 'GS', etc.
name : user-visible name to be followed by 'method' in output
'Newton''s', 'Gauss-Seidel', etc.
default_max_it : default value for max_it option, if not provided
need_jac : 1 - compute Jacobian, 0 - do not compute Jacobian
update_fcn : handle to function used to update x, with syntax:
x = update_fcn(x, f)
x = update_fcn(x, f, J)
FCN : handle to function that evaluates the function f(x) to
be solved and (optionally) 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 (SP.default_max_it) - maximum number of iterations
tol (1e-8) - tolerance on Inf-norm of f(x)
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 (SP.alg)
iterations - number of iterations performed
hist - struct array with trajectories of the following:
normf
normdx
message - exit message
JAC : final Jacobian matrix, J
Calling syntax options:
[x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0);
[x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0, opt);
x = nleqs_core(...);
[x, f] = nleqs_core(...);
[x, f, exitflag] = nleqs_core(...);
[x, f, exitflag, output] = nleqs_core(...);
[x, f, exitflag, output, jac] = nleqs_core(...);
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
fcn = @(x)f1(x);
x0 = [0; 0];
opt = struct('verbose', 2);
sp = struct( ...
'alg', 'NEWTON', ...
'name', 'Newton''s', ...
'default_max_it', 10, ...
'need_jac', 1, ...
'update_fcn', @newton_update_fcn, ...
)
[x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0, opt);
See also NLEQS_MASTER, NLEQS_NEWTON, NLEQS_GAUSS_SEIDEL.0001 function [x, f, eflag, output, J] = nleqs_core(sp, fcn, x0, opt) 0002 %NLEQS_CORE Core Nonlinear Equation Solver with arbitrary update function. 0003 % [X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_CORE(SP, FCN, X0, OPT) 0004 % A function providing the core code for a standardized interface for 0005 % various methods of solving the nonlinear equation f(x) = 0, beginning 0006 % from a starting point x0. Allows for an arbitrary update function. 0007 % 0008 % Inputs: 0009 % SP : solver parameters, struct with the following fields (all required) 0010 % alg : algorithm code, e.g. 'NEWTON', 'GS', etc. 0011 % name : user-visible name to be followed by 'method' in output 0012 % 'Newton''s', 'Gauss-Seidel', etc. 0013 % default_max_it : default value for max_it option, if not provided 0014 % need_jac : 1 - compute Jacobian, 0 - do not compute Jacobian 0015 % update_fcn : handle to function used to update x, with syntax: 0016 % x = update_fcn(x, f) 0017 % x = update_fcn(x, f, J) 0018 % FCN : handle to function that evaluates the function f(x) to 0019 % be solved and (optionally) its Jacobian, J(x). Calling syntax 0020 % for this function is: 0021 % f = FCN(x) 0022 % [f, J] = FCN(x) 0023 % If f and x are n x 1, then J is the n x n matrix of partial 0024 % derivatives of f (rows) w.r.t. x (cols). 0025 % X0 : starting value, x0, of vector x 0026 % OPT : optional options structure with the following fields, 0027 % all of which are also optional (default values shown in 0028 % parentheses) 0029 % verbose (0) - controls level of progress output displayed 0030 % 0 = no progress output 0031 % 1 = some progress output 0032 % 2 = verbose progress output 0033 % max_it (SP.default_max_it) - maximum number of iterations 0034 % tol (1e-8) - tolerance on Inf-norm of f(x) 0035 % 0036 % Outputs (all optional, except X): 0037 % X : solution vector x 0038 % F : final function value, f(x) 0039 % EXITFLAG : exit flag 0040 % 1 = converged 0041 % 0 or negative values = solver specific failure codes 0042 % OUTPUT : output struct with the following fields: 0043 % alg - algorithm code of solver used (SP.alg) 0044 % iterations - number of iterations performed 0045 % hist - struct array with trajectories of the following: 0046 % normf 0047 % normdx 0048 % message - exit message 0049 % JAC : final Jacobian matrix, J 0050 % 0051 % Calling syntax options: 0052 % [x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0); 0053 % [x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0, opt); 0054 % x = nleqs_core(...); 0055 % [x, f] = nleqs_core(...); 0056 % [x, f, exitflag] = nleqs_core(...); 0057 % [x, f, exitflag, output] = nleqs_core(...); 0058 % [x, f, exitflag, output, jac] = nleqs_core(...); 0059 % 0060 % Example: (problem from https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/) 0061 % function [f, J] = f1(x) 0062 % f = [ x(1) + x(2) - 1; 0063 % -x(1)^2 + x(2) + 5 ]; 0064 % if nargout > 1 0065 % J = [1 1; -2*x(1) 1]; 0066 % end 0067 % 0068 % fcn = @(x)f1(x); 0069 % x0 = [0; 0]; 0070 % opt = struct('verbose', 2); 0071 % sp = struct( ... 0072 % 'alg', 'NEWTON', ... 0073 % 'name', 'Newton''s', ... 0074 % 'default_max_it', 10, ... 0075 % 'need_jac', 1, ... 0076 % 'update_fcn', @newton_update_fcn, ... 0077 % ) 0078 % [x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0, opt); 0079 % 0080 % See also NLEQS_MASTER, NLEQS_NEWTON, NLEQS_GAUSS_SEIDEL. 0081 0082 % MP-Opt-Model 0083 % Copyright (c) 1996-2020, Power Systems Engineering Research Center (PSERC) 0084 % by Ray Zimmerman, PSERC Cornell 0085 % 0086 % This file is part of MP-Opt-Model. 0087 % Covered by the 3-clause BSD License (see LICENSE file for details). 0088 % See https://github.com/MATPOWER/mp-opt-model for more info. 0089 0090 %% set default options 0091 opt0 = struct( 'verbose', 0, ... 0092 'max_it', sp.default_max_it, ... 0093 'tol', 1e-8 ); 0094 if isempty(opt) 0095 opt = opt0; 0096 end 0097 if isfield(opt, 'verbose') && ~isempty(opt.verbose) 0098 verbose = opt.verbose; 0099 else 0100 verbose = opt0.verbose; 0101 end 0102 if isfield(opt, 'max_it') && opt.max_it %% not empty or zero 0103 max_it = opt.max_it; 0104 else 0105 max_it = opt0.max_it; 0106 end 0107 if isfield(opt, 'tol') && opt.tol %% not empty or zero 0108 tol = opt.tol; 0109 else 0110 tol = opt0.tol; 0111 end 0112 0113 %% initialize 0114 eflag = 0; 0115 i = 0; 0116 x = x0; 0117 hist(max_it+1) = struct('normf', 0, 'normdx', 0); 0118 0119 %% evaluate f(x0) 0120 if sp.need_jac 0121 [f, J] = fcn(x); 0122 else 0123 f = fcn(x); 0124 end 0125 0126 %% check tolerance 0127 normf = norm(f, inf); 0128 if verbose > 1 0129 fprintf('\n it max residual max ∆x'); 0130 fprintf('\n---- -------------- --------------'); 0131 fprintf('\n%3d %10.3e - ', i, normf); 0132 end 0133 if normf < tol 0134 eflag = 1; 0135 msg = sprintf('%s method converged in %d iterations.', sp.name, i); 0136 if verbose > 1 0137 fprintf('\nConverged!\n'); 0138 end 0139 end 0140 0141 %% save history 0142 hist(i+1).normf = normf; 0143 0144 %% do solver iterations 0145 while (~eflag && i < max_it) 0146 %% update iteration counter 0147 i = i + 1; 0148 0149 %% save previous x 0150 xp = x; 0151 0152 if sp.need_jac 0153 %% update x 0154 x = sp.update_fcn(x, f, J); 0155 0156 %% evalute f(x) and J(x) 0157 [f, J] = fcn(x); 0158 else 0159 %% update x 0160 x = sp.update_fcn(x, f); 0161 0162 %% evaluate f(x) 0163 f = fcn(x); 0164 end 0165 0166 %% check for convergence 0167 normf = norm(f, inf); 0168 normdx = norm(x-xp, inf); 0169 if verbose > 1 0170 fprintf('\n%3d %10.3e %10.3e', i, normf, normdx); 0171 end 0172 if normf < tol 0173 eflag = 1; 0174 msg = sprintf('%s method converged in %d iterations.', sp.name, i); 0175 end 0176 0177 %% save history 0178 hist(i+1).normf = normf; 0179 hist(i+1).normfx = normdx; 0180 end 0181 if eflag ~= 1 0182 msg = sprintf('%s method did not converge in %d iterations.', sp.name, i); 0183 end 0184 if verbose 0185 fprintf('\n%s\n', msg); 0186 end 0187 if nargout > 3 0188 output = struct('alg', sp.alg, ... 0189 'iterations', i, ... 0190 'hist', hist(1:i+1), ... 0191 'message', msg ); 0192 if nargout > 4 && ~sp.need_jac 0193 J = []; 0194 end 0195 end