Broyden方法求解非线性方程组的Matlab实现
注：matlab代码来自网络，仅供学习参考。
1.把以下代码复制在一个.m文件上
function [sol, it_hist, ierr] = brsola(x,f,tol, parms)
% Broyden's Method solver, globally convergent
% solver for f(x) = 0, Armijo rule, one vector storage
%
% This code es with no guarantee or warranty of any kind.
%
% function [sol, it_hist, ierr] = brsola(x,f,tol,parms)
%
% inputs:
% initial iterate = x
% function = f
% tol = [atol, rtol] relative/absolute
% error tolerances for the nonlinear iteration
% parms = [maxit, maxdim]
% maxit = maxmium number of nonlinear iterations
% default = 40
% maxdim = maximum number of Broyden iterations
% before restart, so maxdim-1 vectors are
% stored
% default = 40
%
% output:
% sol = solution
% it_hist(maxit,3) = scaled l2 norms of nonlinear residuals % for the iteration, number function evaluations,
% and number of steplength reductions
% ierr = 0 upon successful termination
% ierr = 1 if after maxit iterations
% the termination criterion is not satsified.
% ierr = 2 failure in the line search. The iteration
% is terminated if too many steplength reductions
% are taken.
%
%
% internal parameter:
% debug = turns on/off iteration statistics display as
% the iteration progresses
%
% alpha = 1.d-4, parameter to measure sufficient decrease %
% maxarm = 10, maximum number of steplength reductions before % failure is reported
%
% set the debug parameter, 1 turns display on, otherwise off
%
debug=1;
%
% initialize it_hist, ierr, and set the iteration parameters
%
ierr = 0; maxit=40; maxdim=39;
it_histx=zeros(maxit,3);
maxarm=10;
%
if nargin == 4
maxit=parms(1); maxdim=parms(2)-1;
end
rtol=tol(2); atol=tol(1); n = length(x); fnrm=1; itc=0; nbroy=0; %
% evaluate f at the initial iterate
% pute the stop tolerance
%
f0=feval(f,x);
fc=f0;
fnrm=norm(f0)/sqrt(n);
it_hist(itc+1)=fnrm;
it_histx(itc+1,1)=fnrm; it_histx(itc+1,2)=0;
it_histx(itc+1,3)=0;
fnrmo=1;
stop_tol=atol + rtol*fnrm;
outstat(itc+1, :) = [itc fnrm 0 0];
%
% terminate on entry?
%
if fnrm < stop_tol
sol=x;
return
end
%
% initialize the iteration history storage matrices
%
stp=zeros(n,maxdim);
stp_nrm=zeros(maxdim,1);
lam_rec=ones(maxdim,1);
%
% Set the initial step to -F, pute the step norm
%
lambda=1;
stp(:,1) = -fc;
stp_nrm(1)=stp(:,1)'*stp(:,1);
%
% main iteration loop
%
while(itc < maxit)
%
nbroy=nbroy+1;
%
% keep track of successive residual norms and
% the iteration counter (itc)
%
fnrmo=fnrm; itc=itc+1;
%
% pute the new point, test for termination before
% adding to iteration history
%
xold=x; lambda=1; iarm=0; lrat=.5; alpha=1.d-4;
x = x + stp(:,nbroy);
fc=feval(f,x);
fnrm=norm(fc)/sqrt(n);
ff0=fnrmo*fnrmo; ffc=fnrm*fnrm; lamc=lambda;
%
%
% Line search, we assume that the Broyden direction is an
% ineact Newton direction. If the line search fails to
% find sufficient decrease after maxarm steplength reductions % brsola returns with failure.
%
% Three-point parabolic line search
%
while fnrm >= (1 - lambda*alpha)*fnrmo && iarm < maxarm
% lambda=lambda*lrat;
if iarm==0
lambda=lambda*lrat;
else
lambda=parab3p(lamc, lamm, ff0, ffc, ffm);
