-
Notifications
You must be signed in to change notification settings - Fork 3
/
secant-julia.html
137 lines (116 loc) · 7.38 KB
/
secant-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>Secant method in Julia</title>
<meta charset="utf-8">
<meta property="og:title" content="Secant 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/secant-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/secant-julia">
<meta name="twitter:card" content="summary">
<meta name="twitter:domain" content="mmas.github.io">
<meta name="twitter:title" content="Secant method in Julia">
<meta name="description" content="The secant method is a root-finding algorithm that recursively calculates the roots of secant lines of points defined in \(f\). Starting with initial...">
<meta name="twitter:description" content="The secant method is a root-finding algorithm that recursively calculates the roots of secant lines of points defined in \(f\). Starting with initial...">
<meta property="og:description" content="The secant method is a root-finding algorithm that recursively calculates the roots of secant lines of points defined in \(f\). Starting with initial...">
<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-20">
<meta property="article:modified_time" content="2016-06-20">
<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="/secant-julia">Secant method in Julia</a></h1>
<time datetime="2016-06-20">Jun 20, 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/Secant_method">The secant method</a> is a root-finding algorithm that recursively calculates the roots of secant lines of points defined in <span class="math">\(f\)</span>.</p>
<p>Starting with initial values <span class="math">\(x_0\)</span> and <span class="math">\(x_1\)</span>, the equation to find the root of the line between <span class="math">\((x_0, f(x_0))\)</span> and <span class="math">\((x_1, f(x_1))\)</span> is</p>
<p><span class="math">\[
0 = {f(x_1)-f(x_0) \over x_1-x_0}(x - x_1) +f(x_1).
\]</span></p>
<p>Solving the equation for <span class="math">\(x\)</span> gives</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>The new <span class="math">\(x\)</span> obtained will be the next <span class="math">\(x_1\)</span> in the recurrence process, which can be expressed as</p>
<p><span class="math">\[
x_n = x_{n-1} - f(x_{n-1}){x_{n-1}-x_{n-2} \over f(x_{n-1})-f(x_{n-2})} \ .
\]</span></p>
<p>The process of solving <span class="math">\(x\)</span> will continue until it reaches a level of tolerance is reached, that is a defined minimum value of the difference between <span class="math">\(x_1\)</span> and <span class="math">\(x\)</span>.</p>
<p>In Julia:</p>
<pre class="sourceCode julia"><code class="sourceCode julia"><span class="kw">function</span> secant(f::<span class="dt">Function</span>, x0::<span class="dt">Number</span>, x1::<span class="dt">Number</span>, args::<span class="dt">Tuple</span>=();
tol::AbstractFloat=<span class="fl">1e-5</span>, maxiter::<span class="dt">Integer</span>=<span class="fl">50</span>)
<span class="kw">for</span> _ <span class="kw">in</span> <span class="fl">1</span>:maxiter
y1 = f(x1, args...)
y0 = f(x0, args...)
x = x1 - y1 * (x1-x0)/(y1-y0)
<span class="kw">if</span> abs(x-x1) < tol
<span class="kw">return</span> x
<span class="kw">end</span>
x0 = x1
x1 = x
<span class="kw">end</span>
error(<span class="st">"Max iteration exceeded"</span>)
<span class="kw">end</span></code></pre>
<p>As an example, the function <span class="math">\(f(x)=x^2-10\)</span> has the roots <span class="math">\(\pm \sqrt{10}\)</span> and with initial guesses <span class="math">\(x_0=1\)</span> and <span class="math">\(x_1=\pm 2\)</span>:</p>
<pre class="sourceCode julia"><code class="sourceCode julia">julia> (f(x)=x^<span class="fl">2</span>-<span class="fl">10</span>; x0=<span class="fl">1</span>; x1=<span class="fl">2</span>; secant(f,x0,x1))
<span class="fl">3.162277660040216</span>
julia> (f(x)=x^<span class="fl">2</span>-<span class="fl">10</span>; x0=<span class="fl">1</span>; x1=-<span class="fl">2</span>; secant(f,x0,x1))
-<span class="fl">3.1622776609633</span>
julia> sqrt(<span class="fl">10</span>)
<span class="fl">3.1622776601683795</span></code></pre>
<p>In the same way of the implementation of the <a href="/newton-julia">Newton-Raphson method</a>, extra arguments <code>args</code> can be passed to the function <code>f</code> and result may be complex numbers if complex numbers are passed as initial guesses. The function <span class="math">\(f(x,c)=x^2+10+c\)</span> with <span class="math">\(c=5\)</span> has the roots <span class="math">\(\pm i \sqrt{15}\)</span>, with initial guesses <span class="math">\(x_0=1\)</span> and <span class="math">\(x_1=2i\)</span></p>
<pre class="sourceCode julia"><code class="sourceCode julia">julia> (f(x,c)=x^<span class="fl">2</span>+<span class="fl">10</span>+c; x0=<span class="fl">1</span>; x1=<span class="fl">2im</span>; secant(f,x0,x1,(<span class="fl">5</span>,)))
-<span class="fl">8.268421911988619e-11</span> + <span class="fl">3.8729833464880765im</span>
julia> sqrt(-<span class="fl">15</span>+<span class="fl">0im</span>)
<span class="fl">0.0</span> + <span class="fl">3.872983346207417im</span></code></pre>
</section>
</article>
</section>
<footer></footer>
<script type="text/javascript" src="/js/article.js"></script>
</body>
</html>