-
Notifications
You must be signed in to change notification settings - Fork 3
Expand file tree
/
Copy pathinverse-quadratic-interpolation-julia.html
More file actions
141 lines (120 loc) · 7.64 KB
/
inverse-quadratic-interpolation-julia.html
File metadata and controls
141 lines (120 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
<!DOCTYPE html>
<html lang="en" itemscope itemtype="http://schema.org/Article">
<head>
<title>Inverse quadratic interpolation in Julia</title>
<meta charset="utf-8">
<meta property="og:title" content="Inverse quadratic interpolation 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/inverse-quadratic-interpolation-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/inverse-quadratic-interpolation-julia">
<meta name="twitter:card" content="summary">
<meta name="twitter:domain" content="mmas.github.io">
<meta name="twitter:title" content="Inverse quadratic interpolation in Julia">
<meta name="description" content="Inverse quadratic interpolation is rarely used directly as a root-finding algorithm, but is is part of the popular Brent's method. In the inverse quad...">
<meta name="twitter:description" content="Inverse quadratic interpolation is rarely used directly as a root-finding algorithm, but is is part of the popular Brent's method. In the inverse quad...">
<meta property="og:description" content="Inverse quadratic interpolation is rarely used directly as a root-finding algorithm, but is is part of the popular Brent's method. In the inverse quad...">
<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-27">
<meta property="article:modified_time" content="2016-06-27">
<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="/inverse-quadratic-interpolation-julia">Inverse quadratic interpolation in Julia</a></h1>
<time datetime="2016-06-27">Jun 27, 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/Inverse_quadratic_interpolation">Inverse quadratic interpolation</a> is rarely used directly as a root-finding algorithm, but is is part of the popular <a target="_blank" href="https://en.wikipedia.org/wiki/Brent's_method">Brent's method</a>.</p>
<p>In the inverse quadratic interpolation, a function <span class="math">\(f(x)\)</span> is evaluated in three points (<span class="math">\(x_0\)</span>, <span class="math">\(x_1\)</span> and <span class="math">\(x_2\)</span>) within an interval where the root of the function is known to be found. Assuming that the function <span class="math">\(f(x)\)</span> has an inverse quadratic function and applying the <a target="_blank" href="https://en.wikipedia.org/wiki/Lagrange_polynomial">Lagrange polynomial</a> for the three points gives the inverse quadratic function</p>
<p><span class="math">\[
g(y) = {(y-f(x_1))(y-f(x_2)) \over (f(x_0)-f(x_1))(f(x_0)-f(x_2))} \ x_0 + {(y-f(x_0))(y-f(x_2)) \over (f(x_1)-f(x_0))(f(x_1)-f(x_2))} \ x_1 + {(y-f(x_0))(y-f(x_1)) \over (f(x_2)-f(x_0))(f(x_2)-f(x_1))} \ x_2 \ .
\]</span></p>
<p>Then, we replace <span class="math">\(y=0\)</span> to obtain the new guess</p>
<p><span class="math">\[
x = {f(x_1)f(x_2) \over (f(x_0)-f(x_1))(f(x_0)-f(x_2))} \ x_0 + {f(x_0)f(x_2) \over (f(x_1)-f(x_0))(f(x_1)-f(x_2))} \ x_1 + {f(x_0)f(x_1) \over (f(x_2)-f(x_0))(f(x_2)-f(x_1))} \ x_2 \ .
\]</span></p>
<p>After each iteration, <span class="math">\(x_1\)</span> will be set as <span class="math">\(x_0\)</span>, <span class="math">\(x_2\)</span> as <span class="math">\(x_1\)</span> and <span class="math">\(x\)</span> as <span class="math">\(x_2\)</span> until the algorithm converges.</p>
<p>In Julia:</p>
<pre class="sourceCode julia"><code class="sourceCode julia"><span class="kw">function</span> invquadinterp(f::<span class="dt">Function</span>, x0::<span class="dt">Number</span>, x1::<span class="dt">Number</span>, x2::<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>)
y0 = f(x0,args...)
y1 = f(x1,args...)
y2 = f(x2,args...)
<span class="kw">for</span> _ <span class="kw">in</span> <span class="fl">1</span>:maxiter
x = x0*y1*y2/((y0-y1)*(y0-y2)) +
x1*y0*y2/((y1-y0)*(y1-y2)) +
x2*y0*y1/((y2-y0)*(y2-y1))
<span class="co"># x-tolerance.</span>
<span class="kw">if</span> min(abs(x-x0),abs(x-x1),abs(x-x2)) < 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>
x0 = x1
y0 = y1
x1 = x2
y1 = y2
x2 = x
y2 = y
<span class="kw">end</span>
error(<span class="st">"Max iteration exceeded"</span>)
<span class="kw">end</span></code></pre>
<p>This implementation has support for extra function arguments and the accuracy can be defined in the <span class="math">\(x\)</span> and <span class="math">\(y\)</span> dimensions.</p>
<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>; invquadinterp(f,<span class="fl">0</span>,<span class="fl">.5</span>,<span class="fl">1</span>))
<span class="fl">0.3660254037449329</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>