]> git.proxmox.com Git - ceph.git/blame - ceph/src/boost/libs/math/doc/html/math_toolkit/roots/root_finding_examples/elliptic_eg.html
projects / ceph.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bump version to 12.2.2-pve1
[ceph.git] / ceph / src / boost / libs / math / doc / html / math_toolkit / roots / root_finding_examples / elliptic_eg.html
CommitLineData
7c673cae
FG		 1<html>
2<head>
3<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
4<title>A More complex example - Inverting the Elliptic Integrals</title>
5<link rel="stylesheet" href="../../../math.css" type="text/css">
6<meta name="generator" content="DocBook XSL Stylesheets V1.77.1">
7<link rel="home" href="../../../index.html" title="Math Toolkit 2.5.1">
8<link rel="up" href="../root_finding_examples.html" title="Examples of Root-Finding (with and without derivatives)">
9<link rel="prev" href="nth_root.html" title="Generalizing to Compute the nth root">
10<link rel="next" href="../bad_guess.html" title="The Effect of a Poor Initial Guess">
11</head>
12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
13<table cellpadding="2" width="100%"><tr>
14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../boost.png"></td>
15<td align="center"><a href="../../../../../../../index.html">Home</a></td>
16<td align="center"><a href="../../../../../../../libs/libraries.htm">Libraries</a></td>
17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
19<td align="center"><a href="../../../../../../../more/index.htm">More</a></td>
20</tr></table>
21<hr>
22<div class="spirit-nav">
23<a accesskey="p" href="nth_root.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding_examples.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="../bad_guess.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a>
24</div>
25<div class="section">
26<div class="titlepage"><div><div><h4 class="title">
27<a name="math_toolkit.roots.root_finding_examples.elliptic_eg"></a><a class="link" href="elliptic_eg.html" title="A More complex example - Inverting the Elliptic Integrals">A
28 More complex example - Inverting the Elliptic Integrals</a>
29</h4></div></div></div>
30<p>
31 The arc length of an ellipse with radii <span class="emphasis"><em>a</em></span> and <span class="emphasis"><em>b</em></span>
32 is given by:
33 </p>
34<pre class="programlisting">L(a, b) = 4aE(k)</pre>
35<p>
36 with:
37 </p>
38<pre class="programlisting">k = &#8730;(1 - b<sup>2</sup>/a<sup>2</sup>)</pre>
39<p>
40 where <span class="emphasis"><em>E(k)</em></span> is the complete elliptic integral of the
41 second kind - see <a class="link" href="../../ellint/ellint_2.html" title="Elliptic Integrals of the Second Kind - Legendre Form">ellint_2</a>.
42 </p>
43<p>
44 Let's suppose we know the arc length and one radii, we can then calculate
45 the other radius by inverting the formula above. We'll begin by encoding
46 the above formula into a functor that our root-finding algorithms can call.
47 </p>
48<p>
49 Note that while not completely obvious from the formula above, the function
50 is completely symmetrical in the two radii - which can be interchanged
51 at will - in this case we need to make sure that <code class="computeroutput"><span class="identifier">a</span>
52 <span class="special">&gt;=</span> <span class="identifier">b</span></code>
53 so that we don't accidentally take the square root of a negative number:
54 </p>
55<pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">typename</span> <span class="identifier">T</span> <span class="special">=</span> <span class="keyword">double</span><span class="special">&gt;</span>
56<span class="keyword">struct</span> <span class="identifier">elliptic_root_functor_noderiv</span>
57<span class="special">{</span> <span class="comment">// Nth root of x using only function - no derivatives.</span>
58 <span class="identifier">elliptic_root_functor_noderiv</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</span> <span class="identifier">arc</span><span class="special">,</span> <span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</span> <span class="identifier">radius</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">m_arc</span><span class="special">(</span><span class="identifier">arc</span><span class="special">),</span> <span class="identifier">m_radius</span><span class="special">(</span><span class="identifier">radius</span><span class="special">)</span>
59 <span class="special">{</span> <span class="comment">// Constructor just stores value a to find root of.</span>
60 <span class="special">}</span>
61 <span class="identifier">T</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</span> <span class="identifier">x</span><span class="special">)</span>
62 <span class="special">{</span>
63 <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sqrt</span><span class="special">;</span>
64 <span class="comment">// return the difference between required arc-length, and the calculated arc-length for an</span>
65 <span class="comment">// ellipse with radii m_radius and x:</span>
66 <span class="identifier">T</span> <span class="identifier">a</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">max</span><span class="special">)(</span><span class="identifier">m_radius</span><span class="special">,</span> <span class="identifier">x</span><span class="special">);</span>
67 <span class="identifier">T</span> <span class="identifier">b</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">min</span><span class="special">)(</span><span class="identifier">m_radius</span><span class="special">,</span> <span class="identifier">x</span><span class="special">);</span>
68 <span class="identifier">T</span> <span class="identifier">k</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="number">1</span> <span class="special">-</span> <span class="identifier">b</span> <span class="special">*</span> <span class="identifier">b</span> <span class="special">/</span> <span class="special">(</span><span class="identifier">a</span> <span class="special">*</span> <span class="identifier">a</span><span class="special">));</span>
69 <span class="keyword">return</span> <span class="number">4</span> <span class="special">*</span> <span class="identifier">a</span> <span class="special">*</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">ellint_2</span><span class="special">(</span><span class="identifier">k</span><span class="special">)</span> <span class="special">-</span> <span class="identifier">m_arc</span><span class="special">;</span>
70 <span class="special">}</span>
71<span class="keyword">private</span><span class="special">:</span>
72 <span class="identifier">T</span> <span class="identifier">m_arc</span><span class="special">;</span> <span class="comment">// length of arc.</span>
73 <span class="identifier">T</span> <span class="identifier">m_radius</span><span class="special">;</span> <span class="comment">// one of the two radii of the ellipse</span>
74<span class="special">};</span> <span class="comment">// template &lt;class T&gt; struct elliptic_root_functor_noderiv</span>
75</pre>
76<p>
77 We'll also need a decent estimate to start searching from, the approximation:
78 </p>
79<pre class="programlisting">L(a, b) &#8776; 4&#8730;(a<sup>2</sup> + b<sup>2</sup>)</pre>
80<p>
81 Is easily inverted to give us what we need, which using derivative-free
82 root finding leads to the algorithm:
83 </p>
84<pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span> <span class="special">=</span> <span class="keyword">double</span><span class="special">&gt;</span>
85<span class="identifier">T</span> <span class="identifier">elliptic_root_noderiv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">radius</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">arc</span><span class="special">)</span>
86<span class="special">{</span> <span class="comment">// return the other radius of an ellipse, given one radii and the arc-length</span>
87 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">std</span><span class="special">;</span> <span class="comment">// Help ADL of std functions.</span>
88 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">;</span> <span class="comment">// For bracket_and_solve_root.</span>
89
90 <span class="identifier">T</span> <span class="identifier">guess</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">arc</span> <span class="special">*</span> <span class="identifier">arc</span> <span class="special">/</span> <span class="number">16</span> <span class="special">-</span> <span class="identifier">radius</span> <span class="special">*</span> <span class="identifier">radius</span><span class="special">);</span>
91 <span class="identifier">T</span> <span class="identifier">factor</span> <span class="special">=</span> <span class="number">1.2</span><span class="special">;</span> <span class="comment">// How big steps to take when searching.</span>
92
93 <span class="keyword">const</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">50</span><span class="special">;</span> <span class="comment">// Limit to maximum iterations.</span>
94 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span> <span class="comment">// Initally our chosen max iterations, but updated with actual.</span>
95 <span class="keyword">bool</span> <span class="identifier">is_rising</span> <span class="special">=</span> <span class="keyword">true</span><span class="special">;</span> <span class="comment">// arc-length increases if one radii increases, so function is rising</span>
96 <span class="comment">// Define a termination condition, stop when nearly all digits are correct, but allow for</span>
97 <span class="comment">// the fact that we are returning a range, and must have some inaccuracy in the elliptic integral:</span>
98 <span class="identifier">eps_tolerance</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;</span> <span class="identifier">tol</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</span><span class="identifier">digits</span> <span class="special">-</span> <span class="number">2</span><span class="special">);</span>
99 <span class="comment">// Call bracket_and_solve_root to find the solution, note that this is a rising function:</span>
100 <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">,</span> <span class="identifier">T</span><span class="special">&gt;</span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">bracket_and_solve_root</span><span class="special">(</span><span class="identifier">elliptic_root_functor_noderiv</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;(</span><span class="identifier">arc</span><span class="special">,</span> <span class="identifier">radius</span><span class="special">),</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">factor</span><span class="special">,</span> <span class="identifier">is_rising</span><span class="special">,</span> <span class="identifier">tol</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span>
101 <span class="comment">// Result is midway between the endpoints of the range:</span>
102 <span class="keyword">return</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special">+</span> <span class="special">(</span><span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span> <span class="special">-</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span><span class="special">)</span> <span class="special">/</span> <span class="number">2</span><span class="special">;</span>
103<span class="special">}</span> <span class="comment">// template &lt;class T&gt; T elliptic_root_noderiv(T x)</span>
104</pre>
105<p>
106 This function generally finds the root within 8-10 iterations, so given
107 that the runtime is completely dominated by the cost of calling the ellliptic
108 integral it would be nice to reduce that count somewhat. We'll try to do
109 that by using a derivative-based method; the derivatives of this function
110 are rather hard to work out by hand, but fortunately <a href="http://www.wolframalpha.com/input/?i=d%2Fda+%5b4+*+a+*+EllipticE%281+-+b%5e2%2Fa%5e2%29%5d" target="_top">Wolfram
111 Alpha</a> can do the grunt work for us to give:
112 </p>
113<pre class="programlisting">d/da L(a, b) = 4(a<sup>2</sup>E(k) - b<sup>2</sup>K(k)) / (a<sup>2</sup> - b<sup>2</sup>)</pre>
114<p>
115 Note that now we have <span class="bold"><strong>two</strong></span> elliptic integral
116 calls to get the derivative, so our functor will be at least twice as expensive
117 to call as the derivative-free one above: we'll have to reduce the iteration
118 count quite substantially to make a difference!
119 </p>
120<p>
121 Here's the revised functor:
122 </p>
123<pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span> <span class="special">=</span> <span class="keyword">double</span><span class="special">&gt;</span>
124<span class="keyword">struct</span> <span class="identifier">elliptic_root_functor_1deriv</span>
125<span class="special">{</span> <span class="comment">// Functor also returning 1st derviative.</span>
126 <span class="identifier">BOOST_STATIC_ASSERT_MSG</span><span class="special">(</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">is_integral</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</span><span class="identifier">value</span> <span class="special">==</span> <span class="keyword">false</span><span class="special">,</span> <span class="string">"Only floating-point type types can be used!"</span><span class="special">);</span>
127
128 <span class="identifier">elliptic_root_functor_1deriv</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</span> <span class="identifier">arc</span><span class="special">,</span> <span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</span> <span class="identifier">radius</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">m_arc</span><span class="special">(</span><span class="identifier">arc</span><span class="special">),</span> <span class="identifier">m_radius</span><span class="special">(</span><span class="identifier">radius</span><span class="special">)</span>
129 <span class="special">{</span> <span class="comment">// Constructor just stores value a to find root of.</span>
130 <span class="special">}</span>
131 <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">,</span> <span class="identifier">T</span><span class="special">&gt;</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</span> <span class="identifier">x</span><span class="special">)</span>
132 <span class="special">{</span>
133 <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sqrt</span><span class="special">;</span>
134 <span class="comment">// Return the difference between required arc-length, and the calculated arc-length for an</span>
135 <span class="comment">// ellipse with radii m_radius and x, plus it's derivative.</span>
136 <span class="comment">// See http://www.wolframalpha.com/input/?i=d%2Fda+[4+*+a+*+EllipticE%281+-+b^2%2Fa^2%29]</span>
137 <span class="comment">// We require two elliptic integral calls, but from these we can calculate both</span>
138 <span class="comment">// the function and it's derivative:</span>
139 <span class="identifier">T</span> <span class="identifier">a</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">max</span><span class="special">)(</span><span class="identifier">m_radius</span><span class="special">,</span> <span class="identifier">x</span><span class="special">);</span>
140 <span class="identifier">T</span> <span class="identifier">b</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">min</span><span class="special">)(</span><span class="identifier">m_radius</span><span class="special">,</span> <span class="identifier">x</span><span class="special">);</span>
141 <span class="identifier">T</span> <span class="identifier">a2</span> <span class="special">=</span> <span class="identifier">a</span> <span class="special">*</span> <span class="identifier">a</span><span class="special">;</span>
142 <span class="identifier">T</span> <span class="identifier">b2</span> <span class="special">=</span> <span class="identifier">b</span> <span class="special">*</span> <span class="identifier">b</span><span class="special">;</span>
143 <span class="identifier">T</span> <span class="identifier">k</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="number">1</span> <span class="special">-</span> <span class="identifier">b2</span> <span class="special">/</span> <span class="identifier">a2</span><span class="special">);</span>
144 <span class="identifier">T</span> <span class="identifier">Ek</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">ellint_2</span><span class="special">(</span><span class="identifier">k</span><span class="special">);</span>
145 <span class="identifier">T</span> <span class="identifier">Kk</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">ellint_1</span><span class="special">(</span><span class="identifier">k</span><span class="special">);</span>
146 <span class="identifier">T</span> <span class="identifier">fx</span> <span class="special">=</span> <span class="number">4</span> <span class="special">*</span> <span class="identifier">a</span> <span class="special">*</span> <span class="identifier">Ek</span> <span class="special">-</span> <span class="identifier">m_arc</span><span class="special">;</span>
147 <span class="identifier">T</span> <span class="identifier">dfx</span> <span class="special">=</span> <span class="number">4</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">a2</span> <span class="special">*</span> <span class="identifier">Ek</span> <span class="special">-</span> <span class="identifier">b2</span> <span class="special">*</span> <span class="identifier">Kk</span><span class="special">)</span> <span class="special">/</span> <span class="special">(</span><span class="identifier">a2</span> <span class="special">-</span> <span class="identifier">b2</span><span class="special">);</span>
148 <span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">make_pair</span><span class="special">(</span><span class="identifier">fx</span><span class="special">,</span> <span class="identifier">dfx</span><span class="special">);</span>
149 <span class="special">}</span>
150<span class="keyword">private</span><span class="special">:</span>
151 <span class="identifier">T</span> <span class="identifier">m_arc</span><span class="special">;</span> <span class="comment">// length of arc.</span>
152 <span class="identifier">T</span> <span class="identifier">m_radius</span><span class="special">;</span> <span class="comment">// one of the two radii of the ellipse</span>
153<span class="special">};</span> <span class="comment">// struct elliptic_root__functor_1deriv</span>
154</pre>
155<p>
156 The root-finding code is now almost the same as before, but we'll make
157 use of Newton-iteration to get the result:
158 </p>
159<pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span> <span class="special">=</span> <span class="keyword">double</span><span class="special">&gt;</span>
160<span class="identifier">T</span> <span class="identifier">elliptic_root_1deriv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">radius</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">arc</span><span class="special">)</span>
161<span class="special">{</span>
162 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">std</span><span class="special">;</span> <span class="comment">// Help ADL of std functions.</span>
163 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">;</span> <span class="comment">// For newton_raphson_iterate.</span>
164
165 <span class="identifier">BOOST_STATIC_ASSERT_MSG</span><span class="special">(</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">is_integral</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</span><span class="identifier">value</span> <span class="special">==</span> <span class="keyword">false</span><span class="special">,</span> <span class="string">"Only floating-point type types can be used!"</span><span class="special">);</span>
166
167 <span class="identifier">T</span> <span class="identifier">guess</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">arc</span> <span class="special">*</span> <span class="identifier">arc</span> <span class="special">/</span> <span class="number">16</span> <span class="special">-</span> <span class="identifier">radius</span> <span class="special">*</span> <span class="identifier">radius</span><span class="special">);</span>
168 <span class="identifier">T</span> <span class="identifier">min</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> <span class="comment">// Minimum possible value is zero.</span>
169 <span class="identifier">T</span> <span class="identifier">max</span> <span class="special">=</span> <span class="identifier">arc</span><span class="special">;</span> <span class="comment">// Maximum possible value is the arc length.</span>
170
171 <span class="comment">// Accuracy doubles at each step, so stop when just over half of the digits are</span>
172 <span class="comment">// correct, and rely on that step to polish off the remainder:</span>
173 <span class="keyword">int</span> <span class="identifier">get_digits</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special">&lt;</span><span class="keyword">int</span><span class="special">&gt;(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</span><span class="identifier">digits</span> <span class="special">*</span> <span class="number">0.6</span><span class="special">);</span>
174 <span class="keyword">const</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">20</span><span class="special">;</span>
175 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span>
176 <span class="identifier">T</span> <span class="identifier">result</span> <span class="special">=</span> <span class="identifier">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">elliptic_root_functor_1deriv</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;(</span><span class="identifier">arc</span><span class="special">,</span> <span class="identifier">radius</span><span class="special">),</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">max</span><span class="special">,</span> <span class="identifier">get_digits</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span>
177 <span class="keyword">return</span> <span class="identifier">result</span><span class="special">;</span>
178<span class="special">}</span> <span class="comment">// T elliptic_root_1_deriv Newton-Raphson</span>
179</pre>
180<p>
181 The number of iterations required for <code class="computeroutput"><span class="keyword">double</span></code>
182 precision is now usually around 4 - so we've slightly more than halved
183 the number of iterations, but made the functor twice as expensive to call!
184 </p>
185<p>
186 Interestingly though, the second derivative requires no more expensive
187 elliptic integral calls than the first does, in other words it comes essentially
188 "for free", in which case we might as well make use of it and
189 use Halley-iteration. This is quite a typical situation when inverting
190 special-functions. Here's the revised functor:
191 </p>
192<pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span> <span class="special">=</span> <span class="keyword">double</span><span class="special">&gt;</span>
193<span class="keyword">struct</span> <span class="identifier">elliptic_root_functor_2deriv</span>
194<span class="special">{</span> <span class="comment">// Functor returning both 1st and 2nd derivatives.</span>
195 <span class="identifier">BOOST_STATIC_ASSERT_MSG</span><span class="special">(</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">is_integral</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</span><span class="identifier">value</span> <span class="special">==</span> <span class="keyword">false</span><span class="special">,</span> <span class="string">"Only floating-point type types can be used!"</span><span class="special">);</span>
196
197 <span class="identifier">elliptic_root_functor_2deriv</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</span> <span class="identifier">arc</span><span class="special">,</span> <span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</span> <span class="identifier">radius</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">m_arc</span><span class="special">(</span><span class="identifier">arc</span><span class="special">),</span> <span class="identifier">m_radius</span><span class="special">(</span><span class="identifier">radius</span><span class="special">)</span> <span class="special">{}</span>
198 <span class="identifier">std</span><span class="special">::</span><span class="identifier">tuple</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">,</span> <span class="identifier">T</span><span class="special">,</span> <span class="identifier">T</span><span class="special">&gt;</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</span> <span class="identifier">x</span><span class="special">)</span>
199 <span class="special">{</span>
200 <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sqrt</span><span class="special">;</span>
201 <span class="comment">// Return the difference between required arc-length, and the calculated arc-length for an</span>
202 <span class="comment">// ellipse with radii m_radius and x, plus it's derivative.</span>
203 <span class="comment">// See http://www.wolframalpha.com/input/?i=d^2%2Fda^2+[4+*+a+*+EllipticE%281+-+b^2%2Fa^2%29]</span>
204 <span class="comment">// for the second derivative.</span>
205 <span class="identifier">T</span> <span class="identifier">a</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">max</span><span class="special">)(</span><span class="identifier">m_radius</span><span class="special">,</span> <span class="identifier">x</span><span class="special">);</span>
206 <span class="identifier">T</span> <span class="identifier">b</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">min</span><span class="special">)(</span><span class="identifier">m_radius</span><span class="special">,</span> <span class="identifier">x</span><span class="special">);</span>
207 <span class="identifier">T</span> <span class="identifier">a2</span> <span class="special">=</span> <span class="identifier">a</span> <span class="special">*</span> <span class="identifier">a</span><span class="special">;</span>
208 <span class="identifier">T</span> <span class="identifier">b2</span> <span class="special">=</span> <span class="identifier">b</span> <span class="special">*</span> <span class="identifier">b</span><span class="special">;</span>
209 <span class="identifier">T</span> <span class="identifier">k</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="number">1</span> <span class="special">-</span> <span class="identifier">b2</span> <span class="special">/</span> <span class="identifier">a2</span><span class="special">);</span>
210 <span class="identifier">T</span> <span class="identifier">Ek</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">ellint_2</span><span class="special">(</span><span class="identifier">k</span><span class="special">);</span>
211 <span class="identifier">T</span> <span class="identifier">Kk</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">ellint_1</span><span class="special">(</span><span class="identifier">k</span><span class="special">);</span>
212 <span class="identifier">T</span> <span class="identifier">fx</span> <span class="special">=</span> <span class="number">4</span> <span class="special">*</span> <span class="identifier">a</span> <span class="special">*</span> <span class="identifier">Ek</span> <span class="special">-</span> <span class="identifier">m_arc</span><span class="special">;</span>
213 <span class="identifier">T</span> <span class="identifier">dfx</span> <span class="special">=</span> <span class="number">4</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">a2</span> <span class="special">*</span> <span class="identifier">Ek</span> <span class="special">-</span> <span class="identifier">b2</span> <span class="special">*</span> <span class="identifier">Kk</span><span class="special">)</span> <span class="special">/</span> <span class="special">(</span><span class="identifier">a2</span> <span class="special">-</span> <span class="identifier">b2</span><span class="special">);</span>
214 <span class="identifier">T</span> <span class="identifier">dfx2</span> <span class="special">=</span> <span class="number">4</span> <span class="special">*</span> <span class="identifier">b2</span> <span class="special">*</span> <span class="special">((</span><span class="identifier">a2</span> <span class="special">+</span> <span class="identifier">b2</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">Kk</span> <span class="special">-</span> <span class="number">2</span> <span class="special">*</span> <span class="identifier">a2</span> <span class="special">*</span> <span class="identifier">Ek</span><span class="special">)</span> <span class="special">/</span> <span class="special">(</span><span class="identifier">a</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">a2</span> <span class="special">-</span> <span class="identifier">b2</span><span class="special">)</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">a2</span> <span class="special">-</span> <span class="identifier">b2</span><span class="special">));</span>
215 <span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span><span class="identifier">fx</span><span class="special">,</span> <span class="identifier">dfx</span><span class="special">,</span> <span class="identifier">dfx2</span><span class="special">);</span>
216 <span class="special">}</span>
217<span class="keyword">private</span><span class="special">:</span>
218 <span class="identifier">T</span> <span class="identifier">m_arc</span><span class="special">;</span> <span class="comment">// length of arc.</span>
219 <span class="identifier">T</span> <span class="identifier">m_radius</span><span class="special">;</span> <span class="comment">// one of the two radii of the ellipse</span>
220<span class="special">};</span>
221</pre>
222<p>
223 The actual root-finding code is almost the same as before, except we can
224 use Halley, rather than Newton iteration:
225 </p>
226<pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span> <span class="special">=</span> <span class="keyword">double</span><span class="special">&gt;</span>
227<span class="identifier">T</span> <span class="identifier">elliptic_root_2deriv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">radius</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">arc</span><span class="special">)</span>
228<span class="special">{</span>
229 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">std</span><span class="special">;</span> <span class="comment">// Help ADL of std functions.</span>
230 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">;</span> <span class="comment">// For halley_iterate.</span>
231
232 <span class="identifier">BOOST_STATIC_ASSERT_MSG</span><span class="special">(</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">is_integral</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</span><span class="identifier">value</span> <span class="special">==</span> <span class="keyword">false</span><span class="special">,</span> <span class="string">"Only floating-point type types can be used!"</span><span class="special">);</span>
233
234 <span class="identifier">T</span> <span class="identifier">guess</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">arc</span> <span class="special">*</span> <span class="identifier">arc</span> <span class="special">/</span> <span class="number">16</span> <span class="special">-</span> <span class="identifier">radius</span> <span class="special">*</span> <span class="identifier">radius</span><span class="special">);</span>
235 <span class="identifier">T</span> <span class="identifier">min</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> <span class="comment">// Minimum possible value is zero.</span>
236 <span class="identifier">T</span> <span class="identifier">max</span> <span class="special">=</span> <span class="identifier">arc</span><span class="special">;</span> <span class="comment">// radius can't be larger than the arc length.</span>
237
238 <span class="comment">// Accuracy triples at each step, so stop when just over one-third of the digits</span>
239 <span class="comment">// are correct, and the last iteration will polish off the remaining digits:</span>
240 <span class="keyword">int</span> <span class="identifier">get_digits</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special">&lt;</span><span class="keyword">int</span><span class="special">&gt;(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</span><span class="identifier">digits</span> <span class="special">*</span> <span class="number">0.4</span><span class="special">);</span>
241 <span class="keyword">const</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">20</span><span class="special">;</span>
242 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span>
243 <span class="identifier">T</span> <span class="identifier">result</span> <span class="special">=</span> <span class="identifier">halley_iterate</span><span class="special">(</span><span class="identifier">elliptic_root_functor_2deriv</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;(</span><span class="identifier">arc</span><span class="special">,</span> <span class="identifier">radius</span><span class="special">),</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">max</span><span class="special">,</span> <span class="identifier">get_digits</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span>
244 <span class="keyword">return</span> <span class="identifier">result</span><span class="special">;</span>
245<span class="special">}</span> <span class="comment">// nth_2deriv Halley</span>
246</pre>
247<p>
248 While this function uses only slightly fewer iterations (typically around
249 3) to find the root, compared to the original derivative-free method, we've
250 moved from 8-10 elliptic integral calls to 6.
251 </p>
252<p>
253 Full code of this example is at <a href="../../../../../example/root_elliptic_finding.cpp" target="_top">root_elliptic_finding.cpp</a>.
254 </p>
255</div>
256<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
257<td align="left"></td>
258<td align="right"><div class="copyright-footer">Copyright &#169; 2006-2010, 2012-2014 Nikhar Agrawal,
259 Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert
260 Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Johan R&#229;de, Gautam Sewani,
261 Benjamin Sobotta, Thijs van den Berg, Daryle Walker and Xiaogang Zhang<p>
262 Distributed under the Boost Software License, Version 1.0. (See accompanying
263 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
264 </p>
265</div></td>
266</tr></table>
267<hr>
268<div class="spirit-nav">
269<a accesskey="p" href="nth_root.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding_examples.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="../bad_guess.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a>
270</div>
271</body>
272</html>
Ceph