1
2
3
4
5
6
7
8
9
10 clear all;
11
12
13 brentmethod[f_,aa_,bb_,ttolerance_] := Module[{a,b,s,tolerance,temp,c,d,term1,term2,term3,mflag},
14 a=aa;
15 b=bb;
16 s=bb;
17 tolerance=ttolerance;
18 If [Abs[f[a]]<Abs[f[b]],
19 temp=a;
20 a=b;
21 b=temp;
22 ,##&[];
23 ];
24 c = a;
25 mflag=True;
26 While [ f[b]!=0 && f[s]!=0 && Abs[b-a]>tolerance,
27 If [(f[a] != f[c]) && (f[b] != f[c]),
28 term1 = (a*f[b]*f[c])/((f[a]-f[b])*(f[a]-f[c]));
29 term2 = (b*f[a]*f[c]) / ((f[b]-f[a])*(f[b]-f[c]));
30 term3 = (c*f[a]*f[b]) / ((f[c]-f[a])*(f[c]-f[b]));
31 s = term1+term2+term3;
32 ,
33 s = b-(f[b]*(b-a)/(f[b]-f[a]));
34 ];
35 If [ (!(((3*a+b)/4<s)&&(s<b))) ||
36 (mflag && ((Abs[s-b]>=Abs[b-c]/2))) ||
37 (!mflag&&(Abs[s-b]>=Abs[c-d]/2)) ||
38 (mflag&&(Abs[b-c]<Abs[tolerance])) ||
39 (!mflag&&(Abs[c-d]<Abs[tolerance])),
40 s = (a+b)/2;
41 mflag=True;
42 ,
43 mflag=False;
44 ];
45 d = c;
46 c = b;
47 If [f[a]*f[s]<0,
48 b = s;
49 ,
50 a = s;
51 ];
52 If [Abs[f[a]] < Abs[f[b]],
53 temp=a;
54 a=b;
55 b=temp;
56 ,
57 ##&[];
58 ];
59 ];
60 Return[s];
61 ];
62
63
64 f[x_] := (x+4)*(x+1);
65 a=-5.;
66 b=-3.;
67 tolerance=0.00000002;
68 s=brentmethod[f,a,b,tolerance];
69 Print["root=",s];
70 ClearAll;
71 g[x_] := Sin[x];
72 a=3.*Pi/4.;
73 b=5.*Pi/4.;
74 tolerance=0.00000002;
75 s=brentmethod[g,a,b,tolerance];
76 Print["root=",s];
77