% reduced_Newton.m

%

% This MATLAB m-file uses a reduced-Newton algorithm with a

% weak line search to solve a set of non-linear algebraic

% equations.

%

% The input parameters are :

%

% x0 = a column vector of the initial guess of the unknowns

%

% calc_f = the name of a MATLAB function that calculates

%     the function vector

%

% calc_Jac = the name of a function that calculates the Jacobian

%

% Options = a data structure containing optional flags

% .max_iter = max # of Newton's method iterations

% .max_iter_LS = max # of weak line search iterations

% .rtol = relative tolerance

% .atol = absolute tolerance

% .step_tol = abs. tolerance below which we switch to full Newton's method

% .verbose = return a trajectory matrices containing the history

%            of the Newton's method iterations

% .use_range = if non-zero, limit the maximum magitude of the full Newton

%             step so that the change in each component is not greater than

%             that in the vector .range

% .range = a vector of the ranges for each of the unknowns.  Each component

%          of the Newton step

%

% Param = a data structure containing parameters that are to be passed to

%         the calc_f and calc_Jac functions

%

% The output parameters are :

%

% x = the final estimate of the solution

%

% iflag = an integer flag that is 1 for convergence,

%      0 for no convergence, and negative for an error

%

% iter_conv = number of iterations required for convergence

%

% x_traj = a matrix where row # j is the solution estimate at iteration j-1

% f_traj = a matrix where row # j is the function vector at iteration j-1

function [x,iflag,iter_conv,x_traj,f_traj] = ...

    reduced_Newton(x0,calc_f,calc_Jac,Options,Param);

% First, signal no convergence.

iflag = 0;

% Set number of iterations required for convergence.

iter_conv = 0;

% Extract number of state variables.

Nvar = length(x0);

% Initialize solution estimate.

x = x0;

% Calculate initial function vector.

f = feval(calc_f,x,Param);

if(length(f) ~= Nvar)

    iflag = -1;

    error('reduced_Newton: calc_f returns vector of improper length');

end

% ensure f is a column vector

if(size(f,1)~=Nvar)

    f = f';

end

% Obtain initial norm of the function vector for later

% convergence tests.

f0_norm_inf = max(abs(f));

f_norm_2sq = dot(f,f);

% Record initial state and function vectors in trajectory.

count_traj = 1;

x_traj(count_traj,:) = x';

f_traj(count_traj,:) = f';

% Set the flag telling us to perform weak line searches.

i_do_LS = 1;

% Begin Newton's method iterations

for iter = 1:Options.max_iter 

    % calculate the Jacobian

    Jac = feval(calc_Jac,x,Param);

    % Solve the set of linear equations for the full line step

    try

        p = Jac\(-f);

    catch

        iflag = -2;

        error('reduced_Newton: full Newton step calculation error');

    end

    % Now, reduce the magnitude of the Newton step if the user has

    % specified a maximum change allowable for each component.

    if(Options.use_range)

        % Calculate the unit vector lying in the Newton line search

        % direction.

        p_length = norm(p,2);

        p_unit = p/p_length;

        % Calculate the maximum step in this direction allowable under

        % the condition that each state variable must not change by

        % a magnitude greater than the specified range for that variable.

        step_allow = max(abs(Options.range));

        for ivar=1:Nvar

            try

                step_ivar = abs(Options.range(ivar)/p_unit(ivar));

                if(step_ivar < step_allow)

                    step_allow = step_ivar;

                end

            end

        end

        step_allow = min(step_allow,p_length);

        p = p_unit*step_allow;

    end

    % Begin the weak line search

    if(i_do_LS)   % perform a weak line search

        for iter_LS = 0:Options.max_iter_LS           

            iconv_LS = 0;

            % Calculate fractional step length

            lambda = 2^(-iter_LS);

            % Calculate new solution estimate

            x_new = x + lambda*p;

            % Calculate function at the new solution estimate

            f_new = feval(calc_f,x_new,Param);

            % Check descent criterion

            f_new_norm_2sq = dot(f_new,f_new);

            if(f_new_norm_2sq <= f_norm_2sq)

                x = x_new;

                f = f_new;

                f_norm_2sq = f_new_norm_2sq;

                iconv_LS = 1;

                break;

            end

        end

        % If we did not satisfy descent condition, update

        % with final result.

        if(~iconv_LS)

            x = x_new;

            f = f_new;

            f_norm_2sq = f_new_norm_2sq;

        end

    else  % use full Newton step instead

        % Calculate new solution estimate

        x = x + p;

        % Calculate function at the new solution estimate

        f = feval(calc_f,x,Param);

    end

    % if in verbose mode, record state and function vectors

    if(Options.verbose)

        count_traj = count_traj + 1;

        x_traj(count_traj,:) = x';

        f_traj(count_traj,:) = f';

    end

    % check for convergence to the solution

    f_norm_inf = max(abs(f));

    i_conv_rel = 0;

    if(f_norm_inf <= Options.rtol*f0_norm_inf)

        i_conv_rel = 1;

    end

    i_conv_abs = 0;

    if(f_norm_inf <= Options.atol)

        i_conv_abs = 1;

    end

    if((i_conv_rel==1)&(i_conv_abs==1))

        iter_conv = iter;

        iflag = 1;

        break;

    end

    % Check to see whether need to perform a line search

    % at the next step.

    if(f_norm_inf <= Options.step_tol)

        i_do_LS = 0;

    else

        i_do_LS = 1;

    end

end

return;