Broyden方法求解非線性方程組的Matlab實現(xiàn)

本文格式為Word版，下載可任意編輯——Broyden方法求解非線性方程組的Matlab實現(xiàn)

Broyden方法求解非線性方程組的Matlab實現(xiàn)

注：matlab代碼來自網(wǎng)絡，僅供學習參考。

1.把以下代碼復制在一個.m文件上

function[sol,it_hist,ierr]=brsola(x,f,tol,parms)%Broyden'sMethodsolver,globallyconvergent

%solverforf(x)=0,Armijorule,onevectorstorage%

%Thiscodecomeswithnoguaranteeorwarrantyofanykind.%

%function[sol,it_hist,ierr]=brsola(x,f,tol,parms)%

%inputs:

%initialiterate=x

%function=f

%tol=[atol,rtol]relative/absolute

%errortolerancesforthenonlineariteration%parms=[maxit,maxdim]

%maxit=maxmiumnumberofnonlineariterations%default=40

%maxdim=maximumnumberofBroydeniterations%beforerestart,somaxdim-1vectorsare%stored

%default=40

%

%output:

%sol=solution

%it_hist(maxit,3)=scaledl2normsofnonlinearresiduals%fortheiteration,numberfunctionevaluations,%andnumberofsteplengthreductions

%ierr=0uponsuccessfultermination

%ierr=1ifaftermaxititerations

%theterminationcriterionisnotsatsified.%ierr=2failureinthelinesearch.Theiteration%isterminatediftoomanysteplengthreductions%aretaken.

%

%

%internalparameter:

%debug=turnson/offiterationstatisticsdisplayas%theiterationprogresses

%

%alpha=1.d-4,parametertomeasuresufficientdecrease%

%maxarm=10,maximumnumberofsteplengthreductionsbefore%failureisreported

%

%setthedebugparameter,1turnsdisplayon,otherwiseoff%

debug=1;

%

%initializeit_hist,ierr,andsettheiterationparameters%

ierr=0;maxit=40;maxdim=39;

it_histx=zeros(maxit,3);

maxarm=10;

%

ifnargin==4

maxit=parms(1);maxdim=parms(2)-1;

end

rtol=tol(2);atol=tol(1);n=length(x);fnrm=1;itc=0;nbroy=0;%

%evaluatefattheinitialiterate

%computethestoptolerance

%

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,:)=[itcfnrm00];

%

%terminateonentry?

%

iffnrmstop_tol

sol=x;

return

end

%

%initializetheiterationhistorystoragematrices

%

stp=zeros(n,maxdim);

stp_nrm=zeros(maxdim,1);

lam_rec=ones(maxdim,1);

%

%Settheinitialstepto-F,computethestepnorm

%

lambda=1;

stp(:,1)=-fc;

stp_nrm(1)=stp(:,1)'*stp(:,1);

%

%mainiterationloop

%

while(itcmaxit)

%

nbroy=nbroy+1;

%

%keeptrackofsuccessiveresidualnormsand

%theiterationcounter(itc)

%

fnrmo=fnrm;itc=itc+1;

%

%computethenewpoint,testforterminationbefore

%addingtoiterationhistory

%

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;

%

%

%Linesearch,weassumethattheBroydendirectionisan%ineactNewtondirection.Ifthelinesearchfailsto

%findsufficientdecreaseaftermaxarmsteplengthreductions%brsolareturnswithfailure.

%

%Three-pointparaboliclinesearch

%

whilefnrm=(1-lambda*alpha)*fnrmoiarmmaxarm%lambda=lambda*lrat;

ifiarm==0

lambda=lambda*lrat;

else

lambda=parab3p(lamc,lamm,ff0,ffc,ffm);

end

lamm=lamc;ffm=ffc;lamc=lambda;

x=xold+lambda*stp(:,nbroy);

fc=feval(f,x);

fnrm=norm(fc)/sqrt(n);

ffc=fnrm*fnrm;

iarm=iarm+1;

end

%

%seterrorflagandreturnonfailureofthelinesearch%

ifiarm==maxarm

disp('Linesearchfailureinbrsola')

ierr=2;

it_hist=it_histx(1:itc+1,:);

sol=xold;

return;

end

%

%Howmanyfunctionevaluationsdidthisiterationrequire?%

it_histx(itc+1,1)=fnrm;

it_histx(itc+1,2)=it_histx(itc,2)+iarm+1;

if(itc==1)it_histx(itc+1,2)=it_histx(itc+1,2)+1;end;it_histx(itc+1,3)=iarm;

%

%terminate?

%

iffnrmstop_tol

sol=x;

rat=fnrm/fnrmo;

outstat(itc+1,:)=[itcfnrmiarmrat];

it_hist=it_histx(1:itc+1,:);

%it_hist(itc+1)=fnrm;

ifdebug==1

disp(outstat(itc+1,:))

end

return

end

%

%

%modifythestepandstepnormifneededtoreflecttheline%search

%

lam_rec(nbroy)=lambda;

iflambda~=1

stp(:,nbroy)=lambda*stp(:,nbroy);

stp_nrm(nbroy)=lambda*lambda*stp_nrm(nbroy);

end

%

%

%it_hist(itc+1)=fnrm;

rat=fnrm/fnrmo;

outstat(itc+1,:)=[itcfnrmiarmrat];

ifdebug==1

disp(outstat(itc+1,:))

end

%

%

%ifthere'sroom,computethenextsearchdirectionandstepnormand

%addtotheiterationhistory

%

ifnbroymaxdim+1

z=-fc;

ifnbroy1

forkbr=1:nbroy-1

ztmp=stp(:,kbr+1)/lam_rec(kbr+1);

ztmp=ztmp+(1-1/lam_rec(kbr))*stp(:,kbr);ztmp=ztmp*lam_rec(kbr);

z=z+ztmp*((stp(:,kbr)'*z)/stp_nrm(kbr));

end

end

%

%storethenewsearchdirectionanditsnorm

%

a2=-lam_rec(nbroy)/stp_nrm(nbroy);

a1=1-lam_rec(nbroy);

zz=stp(:,nbroy)'*z;

a3=a1*zz/stp_nrm(nbroy);

a4=1+a2*zz;

stp(:,nbroy+1)=(z-a3*stp(:,nbroy))/a4;

stp_nrm(nbroy+1)=stp(:,nbroy+1)'*stp(:,nbroy+1);%

%

%

else

%

%outofroom,timetorestart

%

stp(:,1)=-fc;

stp_nrm(1)=stp(:,1)'*stp(:,1);

nbroy=0;

%

%

%

end

%

%endwhile

end

%

%We'renotsupposedtobehere,we'vetakenthemaximum%numberofiterationsandnotterminated.

%

sol=x;

it_hist=it_histx(1:itc+1,:);

ierr=1;

ifdebug==1

disp('outstat')

end

functionlambdap=parab3p(lambdac,lambdam,ff0,ffc,ffm)%Applythree-pointsafeguardedparabolicmodelforalinesearch.%

%Thiscodecomeswithnoguaranteeorwarrantyofanykind.%

%functionlambdap=parab3p(lambdac,lambdam,ff0,ffc,ffm)%

%input:

%lambdac=currentsteplength

%lambdam=previoussteplength

%ff0=valueof\|F(x_c)\|^2

%ffc=valueof\|F(x_c+\lambdacd)\|^2

%ffm=valueof\|F(x_c+\lambdamd)\|^2

%

%output:

%lambdap=newvalueoflambdagivenparabolicmodel%

%internalparameters:

%sigma0=.1,sigma1=.5,safeguardingboundsforthelinesearch

%

%

%setinternalparameters

%

sigma0=.1;sigma1=.5;

%

%computecoefficientsofinterpolationpolynomial

%

%p(lambda)=ff0+(c1lambda+c2lamb