newton copy.m

MDPSAS工具箱是马里兰大学开发的

function A = newton(A,paramlist)

% newton.m A Newton-Raphson procedure for computing solutions
%            to naemodel-derived objects. The method uses a
%            finite-difference approximation for the Jacobian
%            array (see the jacobian.m method) and updates both
%            the 'var' and 'resid' fields of the neqmodel-derived
%            object. Called as
%
%            A = newton(A,paramlist)
%
% INPUT/OUTPUT PARAMETERS
%            A : a naemodel-derived object
%    paramlist : a cell-array of parameter names from A.param;
%                if specified, sensitivity.m is used instead of
%                jacobian.m

% Copyright (c) by Raymond A. Adomaitis, 1998-2004

unpack(A,'solverset');

if isempty(vmin) & isempty(vmax)
    checklim = 0;
else
    checklim = 1;
    vmin = columnnae(A,'solver_vmin');
    vmax = columnnae(A,'solver_vmax');
end

if lsmax < 2
    warning('Minimum value for lsmax = 2; reset to 2 in this object')
    A = set(A,'solverset',2,'lsmax');
    lsmax = get(A,'solverset','lsmax');
end

for inewt = 1:nmax

    % note that the jacobian method updates the residual in object A

    if inewt == 1 | ~mod(inewt,jup) % Newton (if true), modified Newt (false)

        if nargin == 2
            [J, R, param, template] = sensitivity(A,paramlist);
            sjflag = 1;
        else
            [J, R, var, template] = jacobian(A);
            sjflag = 0;
        end

        % 'normalize' each column of J

        for i = 1:size(J,2)
            colnorm(i) = mean(abs(J(:,i)));
            if colnorm(i) < eps
                error(['Column ',int2str(i),' of ',int2str(size(J,2)), ...
                    ' in the Jacobian has vanished'])
            else
                J(:,i) = J(:,i)/colnorm(i);
            end
        end

    else % do not update Jacobian at this step

        A   = residual(A);
        R   = columnnae(A,'resid');
        var = columnnae(A,'var');

    end
    
    if size(J,1) > size(J,2)
        update = -J\R;
    else

        if rcond(J) > sqrt(eps)
            update = -J\R;
        else
            update = linearsystemsolver(-J,R);
        end
    end
    

    if ~(abs(update) >= 0), error('Exiting the newton method'), end
    update = real( update./colnorm' );

    % Perform a line search in the direction of the update vector
    % to determine approximate minimum residual value

    Rnorm(1) = norm(R);
    uplength(1) = 0;
    uplength(2) = 1.0;

    for i = 2:lsmax

        if sjflag
            paramup = col2cell(param + uplength(i)*update,template);
            B = set(A,'param',paramup);
        else

            if checklim

                % Check is current uplength will generate an update that exceeds
                % any specified limits on variable values

                % find indices of updated variables that exceed lower bound
                exlb = find(var + uplength(i)*update <= vmin);
                if ~isempty(exlb)
                    frac = (vmin(exlb)-var(exlb))./(uplength(i)*update(exlb));
                    [minfrac ifrac] = min(abs(frac));
                    uplengthlimit = sign(frac(ifrac))*minfrac*uplength(i);
                    uplength(i) = (uplengthlimit + uplength(i-1))/2;
                end
                % now check upper bound
                exub = find(var + uplength(i)*update >= vmax);
                if ~isempty(exub)
                    frac = (vmax(exub)-var(exub))./(uplength(i)*update(exub));
                    [minfrac ifrac] = min(abs(frac));
                    uplengthlimit = sign(frac(ifrac))*minfrac*uplength(i);
                    uplength(i) = (uplengthlimit + uplength(i-1))/2;
                end

            end % ends "if checklim"

            varup = col2cell(var + uplength(i)*update,template);
            B = set(A,'var',varup);

        end % ends "if sjflag"

        B = residual(B);
        R = columnnae(B,'resid');
        Rnorm(i) = norm(R);

        if Rnorm(i) <= Rnorm(i-1)
            uplength(i+1) = uplength(i) + 2/3*(uplength(i)-uplength(i-1));
        elseif Rnorm(i) > Rnorm(i-1)
            uplength(i+1) = uplength(i) - 2/3*(uplength(i)-uplength(i-1));
        else
            error('Rnorm is NaN - quitting Newton')
        end

    end % ends i loop

    A = B;

    if get(A,'solverset','plot')
        plot(A), drawnow
    end

    disp([ 'Update/residual norm = ', ...
        num2str(norm(uplength(lsmax)*update)),' / ', ...
        num2str(Rnorm(lsmax)) ])
    if inewt == 1, normup1 = norm(uplength(lsmax)*update); end

    if norm(uplength(lsmax)*update)/normup1 < 1e-12 | ...
            norm(uplength(lsmax)*update) < 1e-12
        break
    end

end