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