-
Notifications
You must be signed in to change notification settings - Fork 0
/
modifiedbrent3032247544.m
135 lines (127 loc) · 3.78 KB
/
modifiedbrent3032247544.m
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
% The Modified Brent's Method
% Based on "A modified Brent's method for finding zeros of functions" by Gautam Wilkins · Ming Gu
% Implemented by Moonwon Lee for Math 128A, Fall 2019
% University of California, Berkeley / Student ID : 303 224 7544
function [root,info] = modifiedbrent3032247544(func,Int,params)
% On input:
% func : is a function handle
% Int : is the initial interval containing a root
% params = an object that contains at least three fields
% params.root_tol : should terminate once the interval containg the root is at most params.root_tol in length
% params.func_tol : or the function value at the current iterate is at most params.func_tol in absolute value.
% params.maxit
% On output :
% root = the computed root
% info = should have at least one field info.flag
% info.flag = 0 is 0 for a successful execution, and 1 otherwise.
% info.num_func_call = the number of the given function gets called.
% info.num_it = the number of iteration.
% info.debug is one of 0, 1, 2 and is for debugging purposes.
% Setup
a = Int.a;
b = Int.b;
fa = func(a); info.num_func_call = 1;
fb = func(b); info.num_func_call = 2;
info.flag = 1;
% Optimization
if abs(fa) < abs(fb)
[b, a] = deal(a,b);
[fb, fa] = deal(fa,fb);
end
% The Brent's method uses bisection to guarantee convergence.
% So it only works perfectly when f(a)f(b) <= 0. Otherwise, make an exception
if sign(fa) == sign(fb)
fprintf('Warning!! The given function must change sign on the given interval. So ouput root would be NaN.\n');
info.flag = 1;
info.num_it = 0;
root = NaN;
return;
end
% setup
% c is b's previous value
c = a;
fc = fa;
countFive = 0;
sizeOfIntByBisec = abs(a - b);
root = b;
fs = fb;
for k = 1 : params.maxit
if(k == params.maxit)
info.num_it = k;
info.debug = 0;
return;
end
% Divided Tol cases for debugging purposes.
if (fs == 0)
info.flag = 0;
info.num_it = k;
info.debug = 0;
return;
elseif (abs(b - a) < params.root_tol)
info.flag = 0;
info.num_it = k;
info.debug = 1;
return;
elseif(abs(fs) < params.func_tol)
info.flag = 0;
info.num_it = k;
info.debug = 2;
return;
end
% Inverse quadratic interpolation
if fa ~= fc && fb ~= fc
s = (a * fb * fc)/ ((fa - fb) * (fa - fc)) + (b * fa * fc)/ ((fb - fa) * (fb - fc)) + (c * fa * fb)/ ((fc - fa) * (fc - fb));
set = 0;
% The Secant's Method
else
s = b - (fb * (b-a) / (fb - fa));
set = 1;
end
fs = func(s); info.num_func_call = info.num_func_call + 1;
% Testing two conditions in the paper to decide whether to use Bisection Method or not.
% Calculating new interval.
if fa * fs < 0
newIntSize = abs(s - a);
else
newIntSize = abs(b - s);
end
% The Bisection Method
% Divided two condition cases for debugging purposes.
% The 1st Condition
if ((sizeOfIntByBisec / 2) < newIntSize)
countFive = countFive + 1;
else
countFive = 1;
end
% The 2nd Condition
if countFive == 5 && ((sizeOfIntByBisec / 2) < newIntSize)
s = (a + b) / 2;
set = 2;
elseif (abs(fs) > (abs(fc)/2))
s = (a + b) / 2;
set = 2;
end
% Updating
% Updating fs if the Bisection Method is applied.
if set == 2
fs = func(s); info.num_func_call = info.num_func_call + 1;
sizeOfIntByBisec = abs(a - b);
countFive = 0;
end
% Saving/updating root with s
root = s;
% c is the previous b
c = b;
fc = fb;
if fa * fs < 0
b = s ;
fb = fs;
else
a = s ;
fa = fs;
end
if abs(fa) < abs(fb)
[b, a] = deal(a, b);
[fb, fa]= deal(fa, fb);
end
end