推荐-Broyden方法求解非线性方程组的Matlab实现 精品

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);