-
Notifications
You must be signed in to change notification settings - Fork 3
Expand file tree
/
Copy pathbrent-julia.html
More file actions
184 lines (160 loc) · 11.7 KB
/
brent-julia.html
File metadata and controls
184 lines (160 loc) · 11.7 KB
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
<!DOCTYPE html>
<html lang="en" itemscope itemtype="http://schema.org/Article">
<head>
<title>Brent's method in Julia</title>
<meta charset="utf-8">
<meta property="og:title" content="Brent's method in Julia">
<meta property="og:site_name" content="Modesto Mas | Blog">
<meta property="og:image" content="https://mmas.github.io/images/profile.jpg">
<meta property="og:image:width" content="200">
<meta property="og:image:height" content="200">
<meta property="og:url" content="https://mmas.github.io/brent-julia">
<meta property="og:locale" content="en_GB">
<meta name="twitter:image" content="https://mmas.github.io/images/profile.jpg">
<meta name="twitter:url" content="https://mmas.github.io/brent-julia">
<meta name="twitter:card" content="summary">
<meta name="twitter:domain" content="mmas.github.io">
<meta name="twitter:title" content="Brent's method in Julia">
<meta name="description" content="Brent's method or Wijngaarden-Brent-Dekker method is a root-finding algorithm which combines the bisection method, the secant method and inverse quadr...">
<meta name="twitter:description" content="Brent's method or Wijngaarden-Brent-Dekker method is a root-finding algorithm which combines the bisection method, the secant method and inverse quadr...">
<meta property="og:description" content="Brent's method or Wijngaarden-Brent-Dekker method is a root-finding algorithm which combines the bisection method, the secant method and inverse quadr...">
<meta name="keywords" content="julia,numerical-analysis,root-finding">
<meta property="og:type" content="blog">
<meta name="viewport" content="width=device-width, initial-scale=1">
<meta property="og:type" content="article">
<meta property="article:author" content="https://github.com/mmas">
<meta property="article:section" content="julia">
<meta property="article:tag" content="julia,numerical-analysis,root-finding">
<meta property="article:published_time" content="2016-06-29">
<meta property="article:modified_time" content="2016-06-29">
<link rel="stylesheet" type="text/css" href="/css/main.css">
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
CommonHTML: {
scale: 93,
showMathMenu: false
},
tex2jax: {
"inlineMath": [["$","$"], ["\\(","\\)"]]
}
});
</script>
<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_CHTML"></script>
</head>
<body class="entry-detail">
<header>
<div>
<img src="https://mmas.github.io/images/profile.jpg">
<a class="brand" href="/">Modesto Mas</a>
<span>Data/Python/DevOps Engineer</span>
<nav>
<ul>
<li><a href="/tags">Tags</a></li>
<li><a href="https://github.com/mmas/mmas.github.io/issues" target="_blank">Issues</a></li>
</ul>
</nav>
</div>
</header>
<section id="content" role="main">
<article>
<header>
<h1><a href="/brent-julia">Brent's method in Julia</a></h1>
<time datetime="2016-06-29">Jun 29, 2016</time>
<a class="tag" href="/tags?tag=julia">julia</a>
<a class="tag" href="/tags?tag=numerical-analysis">numerical-analysis</a>
<a class="tag" href="/tags?tag=root-finding">root-finding</a>
</header>
<aside id="article-nav"></aside>
<section class="body">
<p><a target="_blank" href="https://en.wikipedia.org/wiki/Brent's_method">Brent's method</a> or Wijngaarden-Brent-Dekker method is a root-finding algorithm which combines <a href="/bisection-method-julia">the bisection method</a>, <a href="/secant-julia">the secant method</a> and <a href="/inverse-quadratic-interpolation-julia">inverse quadratic interpolation</a>. This method always converges as long as the values of the function are computable within a given region containing a root.</p>
<p>Given a function <span class="math">\(f(x)\)</span> and the bracket <span class="math">\([x_0, x_1]\)</span> two new points, <span class="math">\(x_2\)</span> and <span class="math">\(x_3\)</span>, are initialized with the <span class="math">\(x_1\)</span> value. <span class="math">\(x_0\)</span> and <span class="math">\(x_1\)</span> are swapped if <span class="math">\(|f(x_0)| < |f(x_1)|\)</span>.</p>
<p>Then, in each iteration if the evaluation of the points <span class="math">\(x_0\)</span>, <span class="math">\(x_1\)</span> and <span class="math">\(x_2\)</span> are different (according to a certain tolerance) the inverse quadratic interpolation is used to get the new guess <span class="math">\(x\)</span>. Otherwise, the linear interpolation (secant method) is used to obtain the guess. After that, if any of the following conditions are satisfied <span class="math">\(x\)</span> will be redefined using the bisection method:</p>
<ul>
<li>If <span class="math">\(x\)</span> does not lie between <span class="math">\((3x_0+x_1)/4\)</span> and <span class="math">\(x_1\)</span>.</li>
<li>If in the previous iteration the bisection method was used or it is the first iteration and <span class="math">\(|x-x_1| \geq {1 \over 2} \,|x_1-x_2|\)</span>.</li>
<li>If in the previous iteration the bisection method was not used and <span class="math">\(|x-x_1| \geq {1 \over 2} \,|x_2-x_3|\)</span>.</li>
<li>If in the previous iteration the bisection method was used or it is the first iteration and <span class="math">\(|x_1-x_2| < |\delta|\)</span>, where <span class="math">\(\delta\)</span> is a tolerance factor.</li>
<li>If in the previous iteration the bisection method was not used and <span class="math">\(|x_2-x_3| < |\delta|\)</span>, where <span class="math">\(\delta\)</span> is a tolerance factor.</li>
</ul>
<p>We define <span class="math">\(\delta\)</span> as <span class="math">\(2 \epsilon x_1\)</span>, where <span class="math">\(\epsilon\)</span> is <a target="_blank" href="https://en.wikipedia.org/wiki/Machine_epsilon">the machine epsilon</a>.</p>
<p><span class="math">\(x_3\)</span> and <span class="math">\(x_2\)</span> are redefined in each iteration with <span class="math">\(x_2\)</span> and <span class="math">\(x_1\)</span> value, respectively, and the new guess <span class="math">\(x\)</span> will be set as <span class="math">\(x_1\)</span> if <span class="math">\(f(x_0)f(x)<0\)</span> or as <span class="math">\(x_2\)</span> otherwise. Again, <span class="math">\(x_0\)</span> and <span class="math">\(x_1\)</span> are swapped if <span class="math">\(|f(x_0)| < |f(x_1)|\)</span>. The algorithm converges when <span class="math">\(f(x)\)</span> or <span class="math">\(|x_1-x_0|\)</span> are small enough, both according to tolerance factors.</p>
<p>In the following implementation, the inverse quadratic interpolation is applied directly. In other cases, like the implementation in <a target="_blank" href="http://apps.nrbook.com/empanel/index.html#pg=454">Numerical recipes</a>, used for example in <a target="_blank" href="http://www.boost.org/doc/libs/1_55_0/libs/math/doc/html/math_toolkit/internals1/minima.html">Boost</a>, the Lagrange polynomial is reduced defining the variables <span class="math">\(p\)</span>, <span class="math">\(q\)</span>, <span class="math">\(r\)</span>, <span class="math">\(s\)</span> and <span class="math">\(t\)</span> as explained in <a href="http://mathworld.wolfram.com/BrentsMethod.html">MathWorld</a> and <span class="math">\(x\)</span> value is not overwritten with the bisection method, but modified. For simplicity of the code, here the inverse quadratic interpolation is applied directly as in the entry <a href="/inverse-quadratic-interpolation-julia">Inverse quadratic interpolation in Julia</a> and the new guess is overwritten if needed.</p>
<p>In Julia:</p>
<pre class="sourceCode julia"><code class="sourceCode julia"><span class="kw">function</span> brent(f::<span class="dt">Function</span>, x0::<span class="dt">Number</span>, x1::<span class="dt">Number</span>, args::<span class="dt">Tuple</span>=();
xtol::AbstractFloat=<span class="fl">1e-7</span>, ytol=<span class="fl">2</span>eps(<span class="dt">Float64</span>),
maxiter::<span class="dt">Integer</span>=<span class="fl">50</span>)
EPS = eps(<span class="dt">Float64</span>)
y0 = f(x0,args...)
y1 = f(x1,args...)
<span class="kw">if</span> abs(y0) < abs(y1)
<span class="co"># Swap lower and upper bounds.</span>
x0, x1 = x1, x0
y0, y1 = y1, y0
<span class="kw">end</span>
x2 = x0
y2 = y0
x3 = x2
bisection = true
<span class="kw">for</span> _ <span class="kw">in</span> <span class="fl">1</span>:maxiter
<span class="co"># x-tolerance.</span>
<span class="kw">if</span> abs(x1-x0) < xtol
<span class="kw">return</span> x1
<span class="kw">end</span>
<span class="co"># Use inverse quadratic interpolation if f(x0)!=f(x1)!=f(x2)</span>
<span class="co"># and linear interpolation (secant method) otherwise.</span>
<span class="kw">if</span> abs(y0-y2) > ytol && abs(y1-y2) > ytol
x = x0*y1*y2/((y0-y1)*(y0-y2)) +
x1*y0*y2/((y1-y0)*(y1-y2)) +
x2*y0*y1/((y2-y0)*(y2-y1))
<span class="kw">else</span>
x = x1 - y1 * (x1-x0)/(y1-y0)
<span class="kw">end</span>
<span class="co"># Use bisection method if satisfies the conditions.</span>
delta = abs(<span class="fl">2</span>EPS*abs(x1))
min1 = abs(x-x1)
min2 = abs(x1-x2)
min3 = abs(x2-x3)
<span class="kw">if</span> (x < (<span class="fl">3</span>x0+x1)/<span class="fl">4</span> && x > x1) ||
(bisection && min1 >= min2/<span class="fl">2</span>) ||
(!bisection && min1 >= min3/<span class="fl">2</span>) ||
(bisection && min2 < delta) ||
(!bisection && min3 < delta)
x = (x0+x1)/<span class="fl">2</span>
bisection = true
<span class="kw">else</span>
bisection = false
<span class="kw">end</span>
y = f(x,args...)
<span class="co"># y-tolerance.</span>
<span class="kw">if</span> abs(y) < ytol
<span class="kw">return</span> x
<span class="kw">end</span>
x3 = x2
x2 = x1
<span class="kw">if</span> sign(y0) != sign(y)
x1 = x
y1 = y
<span class="kw">else</span>
x0 = x
y0 = y
<span class="kw">end</span>
<span class="kw">if</span> abs(y0) < abs(y1)
<span class="co"># Swap lower and upper bounds.</span>
x0, x1 = x1, x0
y0, y1 = y1, y0
<span class="kw">end</span>
<span class="kw">end</span>
error(<span class="st">"Max iteration exceeded"</span>)
<span class="kw">end</span></code></pre>
<p>As an example, for the function <span class="math">\(f(x)=x^4-2x^2+1/4 \quad (0 \leq x \leq 1)\)</span>, the solution is <span class="math">\(\sqrt{1-\sqrt{3}/2}\)</span>:</p>
<pre class="sourceCode julia"><code class="sourceCode julia">julia> (f(x)=x^<span class="fl">4</span>-<span class="fl">2</span>x^<span class="fl">2</span>+<span class="fl">1</span>/<span class="fl">4</span>; x0=<span class="fl">0</span>; x1=<span class="fl">1</span>; brent(f,x0,x1))
<span class="fl">0.3660254037844386</span>
julia> sqrt(<span class="fl">1</span>-sqrt(<span class="fl">3</span>)/<span class="fl">2</span>)
<span class="fl">0.3660254037844387</span></code></pre>
</section>
</article>
</section>
<footer></footer>
<script type="text/javascript" src="/js/article.js"></script>
</body>
</html>