-
Notifications
You must be signed in to change notification settings - Fork 3
/
regula-falsi-julia.html
137 lines (116 loc) · 7.64 KB
/
regula-falsi-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
136
137
<!DOCTYPE html>
<html lang="en" itemscope itemtype="http://schema.org/Article">
<head>
<title>Regula falsi method in Julia</title>
<meta charset="utf-8">
<meta property="og:title" content="Regula falsi 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/regula-falsi-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/regula-falsi-julia">
<meta name="twitter:card" content="summary">
<meta name="twitter:domain" content="mmas.github.io">
<meta name="twitter:title" content="Regula falsi method in Julia">
<meta name="description" content="Regula falsi or false position method is a root-finding algorithm very similar to the secant method. In the regula falsi method, the range \([x_0,x_1]...">
<meta name="twitter:description" content="Regula falsi or false position method is a root-finding algorithm very similar to the secant method. In the regula falsi method, the range \([x_0,x_1]...">
<meta property="og:description" content="Regula falsi or false position method is a root-finding algorithm very similar to the secant method. In the regula falsi method, the range \([x_0,x_1]...">
<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-21">
<meta property="article:modified_time" content="2016-06-21">
<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="/regula-falsi-julia">Regula falsi method in Julia</a></h1>
<time datetime="2016-06-21">Jun 21, 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/False_position_method">Regula falsi or false position method</a> is a root-finding algorithm very similar to the <a href="/secant-julia">secant method</a>. In the regula falsi method, the range <span class="math">\([x_0,x_1]\)</span> where the root is found is redefined in each iteration, depending on the sign of the function evaluation in the new <span class="math">\(x\)</span>, this will be set as the minimum or maximum of the new range. Unlike the secant method, regular position method always converges and usually faster than the <a href="/bisection-method-julia">bisection method</a>.</p>
<p>In each iteration, after calculate the new <span class="math">\(x\)</span> as</p>
<p><span class="math">\[
x = x_1 - f(x_1) {x_1-x_0 \over f(x_1)-f(x_0)} \ ,
\]</span></p>
<p>this will be set as the new <span class="math">\(x_0\)</span> value in the case of <span class="math">\(f(x)f(x_0)\)</span> is positive and as the new <span class="math">\(x_1\)</span> value otherwise.</p>
<p>In Julia:</p>
<pre class="sourceCode julia"><code class="sourceCode julia"><span class="kw">function</span> falsi(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-5</span>, ytol=<span class="fl">2</span>eps(<span class="dt">Float64</span>),
maxiter::<span class="dt">Integer</span>=<span class="fl">50</span>)
y1 = f(x1,args...)
y0 = f(x0,args...)
<span class="kw">for</span> _ <span class="kw">in</span> <span class="fl">1</span>:maxiter
x = x1 - y1 * (x1-x0)/(y1-y0)
<span class="co"># x-tolerance.</span>
<span class="kw">if</span> min(abs(x-x0),abs(x-x1)) < xtol
<span class="kw">return</span> x
<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>
<span class="kw">if</span> sign(y0*y)==<span class="fl">1</span>
x0 = x
y0 = y
<span class="kw">else</span>
x1 = x
y1 = y
<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>The main differences in the structure of the regula falsi method implementation with the secant method implementation are that we reassign the variables according to the sign of <span class="math">\(f(x_0)f(x)\)</span> and that the first function evaluations are outside the loop (since we need to check the sign of the function evaluated in <span class="math">\(x\)</span> there's no point in evaluating the function again, just save the result in a variable and use it again at the beginning of the loop).</p>
<p>Also, along with the accuracy in the <span class="math">\(x\)</span> axis (<code>xtol</code>, the minimum difference between two evaluated <span class="math">\(x\)</span>) we added accuracy in the <span class="math">\(y\)</span> axis (<code>ytol</code>, the minimum to be assumed as zero), which by default is set to the machine precission times two (the same default used in <a target="_blank" href="http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#root-finding">root finding algorithms implementations in Scipy</a>).</p>
<p>As an example, for the equation <span class="math">\(f(x)=x^3-23\)</span>, the solution is <span class="math">\(x=\sqrt[3]{23}\)</span>:</p>
<pre class="sourceCode julia"><code class="sourceCode julia">julia> (f(x)=x^<span class="fl">3</span>-<span class="fl">23</span>; x0=<span class="fl">1</span>; x1=<span class="fl">5</span>; falsi(f,x0,x1))
<span class="fl">2.843859313381865</span>
julia> cbrt(<span class="fl">23</span>)
<span class="fl">2.8438669798515654</span></code></pre>
<p>Like in the other <a href="/tags?tag=root-finding">root finding algorithms implementations in this blog</a>, extra arguments and complex numbers can be used.</p>
</section>
</article>
</section>
<footer></footer>
<script type="text/javascript" src="/js/article.js"></script>
</body>
</html>