1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
| % Calcul des racines d'une équation par la méthode de Brent
% Benjamin AUGUSTIN, Dumery Quentin et Tang Lichen - INSA-Lyon - 12 Mai 2009
%
function [xs,n,lint] = brent(f,a,b,tol,nmax)
%
fa=f(a); fb=f(b);
if sign(fa) == sign(fb)
disp('Attention, problème de signe !'),return
end
%
if abs(fa) < abs(fb)
tmp=fa;
fa=fb;
fb=tmp;
end
%
c=a;
mflag=1;
n=0;
while (abs(b-a)>tol+eps*abs(a)) & (n<nmax) & (fb~=0)
n=n+1;
lint(n)=abs(b-a);
fc=f(c);
if (fa~=fc) & (fb~=fc)
xs(n)=(a*fb*fc)/((fa-fb)*(fa-fc))+(b*fa*fc)/((fb-fa)*(fb-fc))+(c*fa*fb)/((fc-fa)*(fc-fb));
else
xs(n)=b-fb*(b-a)/(fb-fa);
end
if (xs(n)<(3*a+b)/4) | (xs(n)>b) | (mflag & (abs(xs(n)-b)>=(abs(b-c)/2))) | (~mflag & (abs(xs(n)-b)>=(abs(c-d)/2)));
xs(n)=(a+b)/2;
mflag=1;
else
mflag=0;
end
fs=f(xs(n));
d=c;
c=b;
if fa*fs<0
b=xs(n);
else
a=xs(n);
end
fa=f(a);
fb=f(b);
if abs(fa)<abs(fb)
tmp=fa;
fa=fb;
fb=tmp;
end
end |
Partager