-
Notifications
You must be signed in to change notification settings - Fork 3
Expand file tree
/
Copy pathridders-julia.html
More file actions
152 lines (131 loc) · 7.64 KB
/
ridders-julia.html
File metadata and controls
152 lines (131 loc) · 7.64 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
<!DOCTYPE html>
<html lang="en" itemscope itemtype="http://schema.org/Article">
<head>
<title>Ridders' method in Julia</title>
<meta charset="utf-8">
<meta property="og:title" content="Ridders' 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/ridders-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/ridders-julia">
<meta name="twitter:card" content="summary">
<meta name="twitter:domain" content="mmas.github.io">
<meta name="twitter:title" content="Ridders' method in Julia">
<meta name="description" content="Ridders' method is a root-finding method based on the regula falsi method that uses an exponential function to fit a given function bracketed between...">
<meta name="twitter:description" content="Ridders' method is a root-finding method based on the regula falsi method that uses an exponential function to fit a given function bracketed between...">
<meta property="og:description" content="Ridders' method is a root-finding method based on the regula falsi method that uses an exponential function to fit a given function bracketed between...">
<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-24">
<meta property="article:modified_time" content="2016-06-24">
<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="/ridders-julia">Ridders' method in Julia</a></h1>
<time datetime="2016-06-24">Jun 24, 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/Ridders'_method">Ridders' method</a> is a root-finding method based on the <a href="/regula-falsi-julia">regula falsi method</a> that uses an exponential function to fit a given function bracketed between <span class="math">\(x_0\)</span> and <span class="math">\(x_1\)</span>.</p>
<p>In the algorithm, in each iteration first the function is evaluated at a third point <span class="math">\(x_2 = (x_0+x_1)/2\)</span>, then applies the unique exponential function</p>
<p><span class="math">\[
f(x_0) -2f(x_2)e^Q +f(x_1)e^{2Q} = 0
\]</span></p>
<p>to transform the function into a straight line.</p>
<p>Solving the quadratic equation for <span class="math">\(e^Q\)</span> gives</p>
<p><span class="math">\[
e^Q = {f(x_2) + \mathrm{sign}(f(x_1)) \sqrt{f(x_2)^2-f(x_0)f(x_1)} \over f(x_1)} \ ,
\]</span></p>
<p>and applying the false position method to <span class="math">\(f(x_0), f(x_2)e^Q\)</span> and <span class="math">\(f(x_1)e^{2Q}\)</span> gives the new guess</p>
<p><span class="math">\[
x = x_2 + (x_2-x_0) {\mathrm{sign}(f(x_0)-x(x_1)) f(x_2) \over \sqrt{f(x_2)^2-f(x_0)f(x_1)}} \ .
\]</span></p>
<p>This equation converges quadratically, therefore the number of significant digits approximately doubles with each iteration (two function evaluations).</p>
<p>After each iteration, <span class="math">\(x_0\)</span> and <span class="math">\(x_1\)</span> will be redefined, being the new guess one of them and the other one defined according to the sign of the evaluations of <span class="math">\(f\)</span> at <span class="math">\(x_0\)</span>, <span class="math">\(x_1\)</span> and <span class="math">\(x_2\)</span>.</p>
<p>In Julia:</p>
<pre class="sourceCode julia"><code class="sourceCode julia"><span class="kw">function</span> ridders(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
x2 = mean([x0,x1])
y2 = f(x2,args...)
x = x2 + (x2-x0) * sign(y0-y1)*y2/sqrt(y2^<span class="fl">2</span>-y0*y1)
<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(y2) != sign(y)
x0 = x2
y0 = y2
x1 = x
y1 = y
<span class="kw">elseif</span> sign(y1) != sign(y)
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 code follows the same structure as the false position method implementation, supports complex numbers and defined accuracy in <span class="math">\(x\)</span> and <span class="math">\(y\)</span> dimensions.</p>
<p>As an example, for the equation <span class="math">\(f(x)=x^2/12+x-4 \quad (1 \leq x \leq 5)\)</span>, the solution is <span class="math">\(x = \sqrt{84}-6\)</span>:</p>
<pre class="sourceCode julia"><code class="sourceCode julia">julia> (f(x)=x^<span class="fl">2</span>/<span class="fl">12</span>+x-<span class="fl">4</span>; x0=<span class="fl">1</span>; x1=<span class="fl">5</span>; ridders(f,x0,x1))
<span class="fl">3.1651513899117836</span>
julia> sqrt(<span class="fl">84</span>)-<span class="fl">6</span>
<span class="fl">3.1651513899116797</span></code></pre>
</section>
</article>
</section>
<footer></footer>
<script type="text/javascript" src="/js/article.js"></script>
</body>
</html>