-
Notifications
You must be signed in to change notification settings - Fork 3
/
bisection-method-julia.html
135 lines (112 loc) · 7.76 KB
/
bisection-method-julia.html
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
<!DOCTYPE html>
<html lang="en" itemscope itemtype="http://schema.org/Article">
<head>
<title>Bisection method in Julia</title>
<meta charset="utf-8">
<meta property="og:title" content="Bisection 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/bisection-method-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/bisection-method-julia">
<meta name="twitter:card" content="summary">
<meta name="twitter:domain" content="mmas.github.io">
<meta name="twitter:title" content="Bisection method in Julia">
<meta name="description" content="The bisection method is a simple root-finding method. Methods for finding roots are iterative and try to find an approximate root \(x\) that fulfills...">
<meta name="twitter:description" content="The bisection method is a simple root-finding method. Methods for finding roots are iterative and try to find an approximate root \(x\) that fulfills...">
<meta property="og:description" content="The bisection method is a simple root-finding method. Methods for finding roots are iterative and try to find an approximate root \(x\) that fulfills...">
<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-04-18">
<meta property="article:modified_time" content="2016-04-18">
<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="/bisection-method-julia">Bisection method in Julia</a></h1>
<time datetime="2016-04-18">Apr 18, 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>The <a href="https://en.wikipedia.org/wiki/Bisection_method">bisection method</a> is a simple root-finding method. Methods for finding roots are iterative and try to find an approximate root <span class="math">\(x\)</span> that fulfills <span class="math">\(|f(x)| \leq \epsilon\)</span>, where <span class="math">\(\epsilon\)</span> is a small number referred later as <em>tolerance</em>.</p>
<p>Being <span class="math">\(f\)</span> a nonlinear equation and <span class="math">\([a,b]\)</span> the interval where the root has to be found, in each iteration, the bisection method algorithm evaluates <span class="math">\(f\)</span> at the midpoint <span class="math">\(c\)</span> of the interval <span class="math">\([a,b]\)</span> and compares it against the evaluation at <span class="math">\(a\)</span> and <span class="math">\(b\)</span>. If <span class="math">\(f(c)=0\)</span> the problem is solved and the root is <span class="math">\(c\)</span>, otherwise the next iteration will continue with the interval <span class="math">\([a,c]\)</span> if <span class="math">\(f(x)\)</span> changes the sign in the left half interval or with <span class="math">\([c,b]\)</span> if it changes in the right half interval. This procedure will continue till find the root in the iteration or, most commonly, by reaching a minimum interval, defined by the tolerance, where the solution is found, get an approximation of the root as the last midpoint <span class="math">\(c\)</span>.</p>
<p>The Julia implementation may be written as follows:</p>
<pre class="sourceCode julia"><code class="sourceCode julia"><span class="kw">function</span> bisection(f::<span class="dt">Function</span>, a::<span class="dt">Number</span>, b::<span class="dt">Number</span>;
tol::AbstractFloat=<span class="fl">1e-5</span>, maxiter::<span class="dt">Integer</span>=<span class="fl">100</span>)
fa = f(a)
fa*f(b) <= <span class="fl">0</span> || error(<span class="st">"No real root in [a,b]"</span>)
i = <span class="fl">0</span>
<span class="kw">local</span> c
<span class="kw">while</span> b-a > tol
i += <span class="fl">1</span>
i != maxiter || error(<span class="st">"Max iteration exceeded"</span>)
c = (a+b)/<span class="fl">2</span>
fc = f(c)
<span class="kw">if</span> fc == <span class="fl">0</span>
<span class="kw">break</span>
<span class="kw">elseif</span> fa*fc > <span class="fl">0</span>
a = c <span class="co"># Root is in the right half of [a,b].</span>
fa = fc
<span class="kw">else</span>
b = c <span class="co"># Root is in the left half of [a,b].</span>
<span class="kw">end</span>
<span class="kw">end</span>
<span class="kw">return</span> c
<span class="kw">end</span></code></pre>
<p>For the equation <span class="math">\(f(x) = x^2-2,\ (x \in [0,2])\)</span>, the solution is <span class="math">\(x = +\sqrt 2 = 1.414213 \dots\)</span></p>
<pre class="sourceCode julia"><code class="sourceCode julia">julia> bisection(x -> x^<span class="fl">2</span>-<span class="fl">2</span>,<span class="fl">0</span>,<span class="fl">2</span>)
<span class="fl">1.4142074584960938</span></code></pre>
<p>In order to avoid multiple evaluations of <code>f(a)</code> and <code>f(c)</code>, its value is cached in the variables <code>fa</code> and <code>fc</code> and used both in the conditionals and the new designation <code>fa = fc</code>. Also, the comparison of the signs of the two halves of the interval is done by <code>fa*fc > 0</code> instead of <code>sign(fa) == sign(fc)</code>, since, in terms of performance, the former uses simliar, but less memory allocation for float values:</p>
<pre class="sourceCode julia"><code class="sourceCode julia">julia> x, y = rand(), -rand()
(<span class="fl">0.48656317905335866</span>,-<span class="fl">0.0722375704738607</span>)
julia> @time sign(x) == sign(y)
<span class="fl">0.000002</span> seconds (<span class="fl">6</span> allocations: <span class="fl">192</span> bytes)
false
julia> @time x*y > <span class="fl">0</span>
<span class="fl">0.000002</span> seconds (<span class="fl">5</span> allocations: <span class="fl">176</span> bytes)
false</code></pre>
</section>
</article>
</section>
<footer></footer>
<script type="text/javascript" src="/js/article.js"></script>
</body>
</html>