moved to https://github.com/mikofski/NewtonRaphson

README.md

NewtonRaphson

Yet another solver that uses the backslash function to solve a set of non-linear equations.

Description

Although this is the most basic non-linear solver, it is surprisingly powerful. It is based on the Newton-Raphson method in chapter 9.6-7 of Numerical Recipes in C. In general for well behaved functions and decent initial guesses, its convergence is at least quadratic. However it may fail if the there are local minimums, the condition of the Jacobian is poor or the initial guess is relatively far from the solution. When convergence is negative, it will attempt to backtrack and line-search. It was validated with fsolve from the MATLAB Optimization Toolbox and IPOPT. Please see the help comments and the example.

Note: LSQ curve-fit type problems can also be solved using newtonraphson. These are problems where there are many data for a single function, but the coefficients of the function are unknown. Since there is more data than unknowns, and the residual is minimized for each data point, the Jacobian is not square. These problems usually exit with flag 2: "May have converged." and the solution is the best fit in the "least-squares" sense.

Credit: Moody diagram is from Wikipedia

Moody_diagram.jpg

newtonraphson.m

function [x, resnorm, F, exitflag, output, jacob] = newtonraphson(fun, x0, options)
% NEWTONRAPHSON Solve set of non-linear equations using Newton-Raphson method.
%
% [X, RESNORM, F, EXITFLAG, OUTPUT, JACOB] = NEWTONRAPHSON(FUN, X0, OPTIONS)
% FUN is a function handle that returns a vector of residuals equations, F,
% and takes a vector, x, as its only argument. When the equations are
% solved by x, then F(x) == zeros(size(F(:), 1)).
%
% Optionally FUN may return the Jacobian, Jij = dFi/dxj, as an additional
% output. The Jacobian must have the same number of rows as F and the same
% number of columns as x. The columns of the Jacobians correspond to d/dxj and
% the rows correspond to dFi/d.
%
%   EG:  J23 = dF2/dx3 is the 2nd row ad 3rd column.
%
% If FUN only returns one output, then J is estimated using a center
% difference approximation,
%
%   Jij = dFi/dxj = (Fi(xj + dx) - Fi(xj - dx))/2/dx.
%
% NOTE: If the Jacobian is not square the system is either over or under
% constrained.
%
% X0 is a vector of initial guesses.
%
% OPTIONS is a structure of solver options created using OPTIMSET.
% EG: options = optimset('TolX', 0.001).
%
% The following options can be set:
% * OPTIONS.TOLFUN is the maximum tolerance of the norm of the residuals.
%   [1e-6]
% * OPTIONS.TOLX is the minimum tolerance of the relative maximum stepsize.
%   [1e-6]
% * OPTIONS.MAXITER is the maximum number of iterations before giving up.
%   [100]
% * OPTIONS.DISPLAY sets the level of display: {'off', 'iter'}.
%   ['iter']
%
% X is the solution that solves the set of equations within the given tolerance.
% RESNORM is norm(F) and F is F(X). EXITFLAG is an integer that corresponds to
% the output conditions, OUTPUT is a structure containing the number of
% iterations, the final stepsize and exitflag message and JACOB is the J(X).
%
% See also OPTIMSET, OPTIMGET, FMINSEARCH, FZERO, FMINBND, FSOLVE, LSQNONLIN
%
% References:
% * http://en.wikipedia.org/wiki/Newton's_method
% * http://en.wikipedia.org/wiki/Newton's_method_in_optimization
% * 9.7 Globally Convergent Methods for Nonlinear Systems of Equations 383,
%   Numerical Recipes in C, Second Edition (1992),
%   http://www.nrbook.com/a/bookcpdf.php

% Version 0.5
% * allow sparse matrices, replace cond() with condest()
% * check if Jstar has NaN or Inf, return NaN or Inf for cond() and return
%   exitflag: -1, matrix is singular.
% * fix bug: max iteration detection and exitflag reporting typos
% Version 0.4
% * allow lsq curve fitting type problems, IE non-square matrices
% * exit if J is singular or if dx is NaN or Inf
% Version 0.3
% * Display RCOND each step.
% * Replace nargout checking in funwrapper with ducktypin.
% * Remove Ftyp and F scaling b/c F(typx)->0 & F/Ftyp->Inf!
% * User Numerical Recipies minimum Newton step, backtracking line search
%   with alpha = 1e-4, min_lambda = 0.1 and max_lambda = 0.5.
% * Output messages, exitflag and min relative step.
% Version 0.2
% * Remove `options.FinDiffRelStep` and `options.TypicalX` since not in MATLAB.
% * Set `dx = eps^(1/3)` in `jacobian` function.
% * Remove `options` argument from `funwrapper` & `jacobian` functions
%   since no longer needed.
% * Set typx = x0; typx(x0==0) = 1; % use initial guess as typx, if 0 use 1.
% * Replace `feval` with `evalf` since `feval` is builtin.

%% initialize
% There are no argument checks!
x0 = x0(:); % needs to be a column vector
% set default options
oldopts = optimset( ...
    'TolX', 1e-12, 'TolFun', 1e-6, 'MaxIter', 100, 'Display', 'iter');
if nargin<3
    options = oldopts; % use defaults
else
    options = optimset(oldopts, options); % update default with user options
end
FUN = @(x)funwrapper(fun, x); % wrap FUN so it always returns J
%% get options
TOLX = optimget(options, 'TolX'); % relative max step tolerance
TOLFUN = optimget(options, 'TolFun'); % function tolerance
MAXITER = optimget(options, 'MaxIter'); % max number of iterations
DISPLAY = strcmpi('iter', optimget(options, 'Display')); % display iterations
TYPX = max(abs(x0), 1); % x scaling value, remove zeros
ALPHA = 1e-4; % criteria for decrease
MIN_LAMBDA = 0.1; % min lambda
MAX_LAMBDA = 0.5; % max lambda
%% set scaling values
% TODO: let user set weights
weight = ones(numel(FUN(x0)),1);
J0 = weight*(1./TYPX'); % Jacobian scaling matrix
%% set display
if DISPLAY
    fprintf('\n%10s %10s %10s %10s %10s %12s\n', 'Niter', 'resnorm', 'stepnorm', ...
        'lambda', 'rcond', 'convergence')
    for n = 1:67,fprintf('-'),end,fprintf('\n')
    fmtstr = '%10d %10.4g %10.4g %10.4g %10.4g %12.4g\n';
    printout = @(n, r, s, l, rc, c)fprintf(fmtstr, n, r, s, l, rc, c);
end
%% check initial guess
x = x0; % initial guess
[F, J] = FUN(x); % evaluate initial guess
Jstar = J./J0; % scale Jacobian
if any(isnan(Jstar(:))) || any(isinf(Jstar(:)))
    exitflag = -1; % matrix may be singular
else
    exitflag = 1; % normal exit
end
if issparse(Jstar)
    rc = 1/condest(Jstar);
else
    if any(isnan(Jstar(:)))
        rc = NaN;
    elseif any(isinf(Jstar(:)))
        rc = Inf;
    else
        rc = 1/cond(Jstar); % reciprocal condition
    end
end
resnorm = norm(F); % calculate norm of the residuals
dx = zeros(size(x0));convergence = Inf; % dummy values
%% solver
Niter = 0; % start counter
lambda = 1; % backtracking
if DISPLAY,printout(Niter, resnorm, norm(dx), lambda, rc, convergence);end
while (resnorm>TOLFUN || lambda<1) && exitflag>=0 && Niter<=MAXITER
    if lambda==1
        %% Newton-Raphson solver
        Niter = Niter+1; % increment counter
        dx_star = -Jstar\F; % calculate Newton step
        % NOTE: use isnan(f) || isinf(f) instead of STPMAX
        dx = dx_star.*TYPX; % rescale x
        g = F'*Jstar; % gradient of resnorm
        slope = g*dx_star; % slope of gradient
        fold = F'*F; % objective function
        xold = x; % initial value
        lambda_min = TOLX/max(abs(dx)./max(abs(xold), 1));
    end
    if lambda<lambda_min
        exitflag = 2; % x is too close to XOLD
        break
    elseif any(isnan(dx)) || any(isinf(dx))
        exitflag = -1; % matrix may be singular
        break
    end
    x = xold+dx*lambda; % next guess
    [F, J] = FUN(x); % evaluate next residuals
    Jstar = J./J0; % scale next Jacobian
    f = F'*F; % new objective function
    %% check for convergence
    lambda1 = lambda; % save previous lambda
    if f>fold+ALPHA*lambda*slope
        if lambda==1
            lambda = -slope/2/(f-fold-slope); % calculate lambda
        else
            A = 1/(lambda1 - lambda2);
            B = [1/lambda1^2,-1/lambda2^2;-lambda2/lambda1^2,lambda1/lambda2^2];
            C = [f-fold-lambda1*slope;f2-fold-lambda2*slope];
            coeff = num2cell(A*B*C);
            [a,b] = coeff{:};
            if a==0
                lambda = -slope/2/b;
            else
                discriminant = b^2 - 3*a*slope;
                if discriminant<0
                    lambda = MAX_LAMBDA*lambda1;
                elseif b<=0
                    lambda = (-b+sqrt(discriminant))/3/a;
                else
                    lambda = -slope/(b+sqrt(discriminant));
                end
            end
            lambda = min(lambda,MAX_LAMBDA*lambda1); % minimum step length
        end
    elseif isnan(f) || isinf(f)
        % limit undefined evaluation or overflow
        lambda = MAX_LAMBDA*lambda1;
    else
        lambda = 1; % fraction of Newton step
    end
    if lambda<1
        lambda2 = lambda1;f2 = f; % save 2nd most previous value
        lambda = max(lambda,MIN_LAMBDA*lambda1); % minimum step length
        continue
    end
    %% display
    resnorm0 = resnorm; % old resnorm
    resnorm = norm(F); % calculate new resnorm
    convergence = log(resnorm0/resnorm); % calculate convergence rate
    stepnorm = norm(dx); % norm of the step
    if any(isnan(Jstar(:))) || any(isinf(Jstar(:)))
        exitflag = -1; % matrix may be singular
        break
    end
    if issparse(Jstar)
        rc = 1/condest(Jstar);
    else
        rc = 1/cond(Jstar); % reciprocal condition
    end
    if DISPLAY,printout(Niter, resnorm, stepnorm, lambda1, rc, convergence);end
end
%% output
output.iterations = Niter; % final number of iterations
output.stepsize = dx; % final stepsize
output.lambda = lambda; % final lambda
if Niter>=MAXITER
    exitflag = 0;
    output.message = 'Number of iterations exceeded OPTIONS.MAXITER.';
elseif exitflag==2
    output.message = 'May have converged, but X is too close to XOLD.';
elseif exitflag==-1
    output.message = 'Matrix may be singular. Step was NaN or Inf.';
else
    output.message = 'Normal exit.';
end
jacob = J;
end

function [F, J] = funwrapper(fun, x)
% if nargout<2 use finite differences to estimate J
try
    [F, J] = fun(x);
catch
    F = fun(x);
    J = jacobian(fun, x); % evaluate center diff if no Jacobian
end
F = F(:); % needs to be a column vector
end

function J = jacobian(fun, x)
% estimate J
dx = eps^(1/3); % finite difference delta
nx = numel(x); % degrees of freedom
nf = numel(fun(x)); % number of functions
J = zeros(nf,nx); % matrix of zeros
for n = 1:nx
    % create a vector of deltas, change delta_n by dx
    delta = zeros(nx, 1); delta(n) = delta(n)+dx;
    dF = fun(x+delta)-fun(x-delta); % delta F
    J(:, n) = dF(:)/dx/2; % derivatives dF/d_n
end
end

newtonraphson_example.m

%% newton raphson example
% Find the Darcy friction factor for pipe flow using the Colebrook
% equation.
fprintf('\n**************************************************\n')
fprintf('NON-LINEAR SYSTEM OF EQUATIONS - PIPE FLOW EXAMPLE\n')
fprintf('**************************************************\n')
%% References:
% * http://en.wikipedia.org/wiki/Darcy_friction_factor_formulae
% * http://en.wikipedia.org/wiki/Darcy%E2%80%93Weisbach_equation
% * http://en.wikipedia.org/wiki/Moody_chart
% * http://www.mathworks.com/matlabcentral/fileexchange/35710-iapwsif97-functional-form-with-no-slip
%% inputs:
p = 0.68; % [MPa] water pressure (100 psi)
dp = -0.068*1e6; % [Pa] pipe pressure drop (10 psi)
T = 323; % [K] water temperature
D = 0.10; % [m] pipe hydraulic diameter
L = 100; % [m] pipe length
roughness = 0.00015; % [m] cast iron pipe roughness
rho = 1./IAPWS_IF97('v_pT',p,T); % [kg/m^3] water density (988.1 kg/m^3)
mu = IAPWS_IF97('mu_pT',p,T); % [Pa*s] water viscosity (5.4790e-04 Pa*s)
Re = @(u) rho*u*D/mu; % Reynolds number
%% governing equations
% Use Colebrook and Darcy-Weisbach equation to solve for pipe flow.

% friction factor (Colebrook eqn.)
residual_friction = @(u, f) 1/sqrt(f) + 2*log10(roughness/3.7/D + 2.51/Re(u)/sqrt(f));
% pressure drop (Darcy-Weisbach eqn.)
residual_pressdrop = @(u, f) rho*u^2*f*L/2/D + dp;
% residuals
fun = @(x) [residual_friction(x(1),x(2)), residual_pressdrop(x(1),x(2))];
%% solve
x0 = [1,0.01]; % initial guess
fprintf('\ninitial guess: u = %g[m/s], f = %g\n',x0) % display initial guess 
options = optimset('TolX',1e-12); % set TolX
[x, resnorm, f, exitflag, output, jacob] = newtonraphson(fun, x0, options);
fprintf('\nexitflag: %d, %s\n',exitflag, output.message) % display output message
%% results
fprintf('\nOutputs:\n')
properties = {'Pressure','Pressure-drop','Temperature','Diameter','Length', ...
    'roughness','density','viscosity','Reynolds-number','speed','friction'};
units = {'[Pa]','[Pa]','[C]','[cm]','[m]','[mm]','[kg/m^3]','[Pa*s]','[]','[m/s]','[]'};
values = {p*1e6,dp,T-273.15,D*100,L,roughness*1000,rho,mu,Re(x(1)),x(1),x(2)};
fprintf('%15s %10s %10s\n','Property','Unit','Value')
results = [properties; units; values];
fprintf('%15s %10s %10.4g\n',results{:})
%% comparison
% solve using Haaland
Ntest = 10;
u0 = linspace(x(1)*0.1, x(1)*10, Ntest); % [m/s]
Re0 = Re(u0);
f0 = (1./(-1.8*log10((roughness/D/3.7)^1.11 + 6.9./Re0))).^2;
u0 = sqrt(-dp/rho./(f0*L/2/D));
% plot
plot(u0, f0, '-', u0, x(2)*ones(1,Ntest), '--', x(1)*ones(1,Ntest), f0, '--')
grid
title('Pipe flow solution using Haaland equation.')
xlabel('water speed, u [m/s]'),ylabel('friction factor, f')
legend('f_{Haaland}',['f = ',num2str(x(2))], ['u = ',num2str(x(1)),' [m/s]'])
%% LSQ Curve Fitting
fprintf('\n**********************************************\n')
fprintf('LEAST SQUARES CURVE FITTING WITH NEWTONRAPHSON\n')
fprintf('**********************************************\n')
% independent variables
[x,y] = meshgrid(0:10,0:10);
% bivariate distribution
bivar = @(x1,x2,sig,u1,u2)1/2/pi/sig^2*exp(-1/2*(((x1-u1).^2+(x2-u2).^2)/sig));
sigma = 3; ux = 4; uy = 5; % std dev, x & y means
z = bivar(x,y,sigma,ux,uy); % dist
% plot
figure,contour(x,y,z),hold('all')
title('lsq curve-fitting bivariate pdf with newtonraphson')
xlabel('x-coord'),ylabel('y-coord') % axes titles
grid,colorbar % show colorbar and grid
z_meas = z + (2*rand(11)-1)/1e4; % measured data
% fitting function
lsqfitfun = @(c)z_meas-bivar(x,y,c(1),c(2),c(3));
% fit coefficients to fun
c0 = [1,2,3]; % initial guess
fprintf('\ninitial guess: sigma = %g, ux = %g, uy = %g\n',c0) % display initial guess 
options = optimset('TolX',1e-12); % set TolX
[c, ~, ~, exitflag, output] = newtonraphson(lsqfitfun, c0, options);
fprintf('\nexitflag: %d, %s\n',exitflag, output.message) % display output message
fprintf('\ncurve-fit coefficients: sigma=%g, ux=%g, uy=%g\n',c)
lines = plot(repmat(c(2),1,11),0:10,'--',0:10,repmat(c(3),1,11),'--', ...
    c(2),c(3),'o');
set(lines,'LineWidth',2);
legend('bivariate distribution','u_x','u_y','center')