%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% secbis.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [x,ier,nf]=secbis(func,x,x2,abserr,prt); % stabilized secant method for finding a zero x of feval(func,x) % if x,x2 form a sign change, the zero will be between x and x2 % % x starting point (default x=0) % x2 second starting point (default x2=0 if x~=0, x2=1 otherwise) % abserr: desired limit for absolute error in x (default 0) % prt=0: no results (default) % prt=1: print intermediate results % % x: computed zero approximation % ier= 0: computed function value f(x)=0 % ier=-1: zero located by sign change % ier= 1: no sign change found; x is point with smallest |f| found % nf: number of function values used % function[x,ier,nf]=secbis(func,x,x2,abserr,prt);
% initialization
nfmax=25; % maximal number of function values if nargin<2, x=0; end;
if nargin<3, if x==0,x2=1;else x2=0; end; end;
if nargin<4,abserr=0; end;
if nargin<5,prt=0; end;
f=feval(func,x);nf=1;
if f==0, ier=0; return; end;
f2=feval(func,x2);nf=2;
if f2==0, x=x2; ier=0; return; end;
if prt==1,
format long
disp('********************************************************') disp(' nf x f slow') disp(sprintf('%3i %20.16f %24.16f %2i',1,x,f,-1)) disp(sprintf('%3i %20.16f %24.16f %2i',2,x2,f2,-1)) end;
% find sign change while f*f2>0,
% place best point in position 2 ifabs(f)<abs(f2), x1=x2;f1=f2;x2=x;f2=f;
else x1=x;f1=f;
end;
% safeguarded root secant formula
x=x2-(x2-x1)/(1-max(sqrt(f1/f2),2));
if x==x2 | nf==nfmax,
% double zero or no sign change found
x=x2;ier=1; return;
end;
f=feval(func,x);nf=nf+1;
disp(sprintf('%3i %20.16f %24.16f %2i',nf,x,f,-1)) if f==0, ier=0; return; end;
end;
if prt>1,
disp(['first sign change at nf = ',num2str(nf)]);
end;
ier=-1; % we have a sign change f*f2<0
slow=0;
while nf<nfmax,
% compute new point xx and accuracy if slow==0,
% standard secant step ifabs(f)<abs(f2),
xx=x-(x-x2)*f/(f-f2);
acc=abs(xx-x);
else
xx=x2-(x2-x)*f2/(f2-f);
acc=abs(xx-x2);
end;
elseif slow==1,
% safeguarded secant extrapolation step if f1*f2>0,
quot=max(f1/f2,2*(x-x1)/(x-x2));
xx=x2-(x2-x1)/(1-quot);
acc=abs(xx-x2);
else% f1*f>0
quot=max(f1/f,2*(x2-x1)/(x2-x));
xx=x-(x-x1)/(1-quot);
acc=abs(xx-x);
end;
else % safeguarded geometric mean step if x*x2>0, xx=x*sqrt(x2/x); % preserves the sign! elseif x==0, xx=0.1*x2;
elseif x2==0, xx=0.1*x;
else xx=0;
end;
acc=max(abs(xx-[x,x2]));
end;
% stopping tests if acc<=abserr, x=xx;ier=-1; return; end;
ff=feval(func,xx);nf=nf+1;
if ff==0, x=xx;ier=0; return; end;
% compute reduction factor and update bracket if f2*ff<0, redfac=(x2-xx)/(x2-x);x1=x;f1=f;x=xx;f=ff;
else redfac=(x-xx)/(x-x2);x1=x2;f1=f2;x2=xx;f2=ff;
end;
% force two consecutive mean steps if nonlocal (slow=2) if redfac>0.7 | slow==2, slow=slow+1;
else slow=0;
end;
if prt>0,
disp(sprintf('%3i %20.16f %24.16f %2i',nf,xx,ff,slow)) end;
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS