]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | <html> |
2 | <head> | |
3 | <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII"> | |
4 | <title>Calculating a Derivative</title> | |
5 | <link rel="stylesheet" href="../../../../../../../../doc/src/boostbook.css" type="text/css"> | |
6 | <meta name="generator" content="DocBook XSL Stylesheets V1.77.1"> | |
7 | <link rel="home" href="../../../../index.html" title="Chapter 1. Boost.Multiprecision"> | |
8 | <link rel="up" href="../fp_eg.html" title="Examples"> | |
9 | <link rel="prev" href="jel.html" title="Defining a Special Function."> | |
10 | <link rel="next" href="gi.html" title="Calculating an Integral"> | |
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="jel.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../fp_eg.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="gi.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a> | |
24 | </div> | |
25 | <div class="section"> | |
26 | <div class="titlepage"><div><div><h5 class="title"> | |
27 | <a name="boost_multiprecision.tut.floats.fp_eg.nd"></a><a class="link" href="nd.html" title="Calculating a Derivative">Calculating | |
28 | a Derivative</a> | |
29 | </h5></div></div></div> | |
30 | <p> | |
31 | In this example we'll add even more power to generic numeric programming | |
32 | using not only different floating-point types but also function objects | |
33 | as template parameters. Consider some well-known central difference rules | |
34 | for numerically computing the first derivative of a function <span class="emphasis"><em>f′(x)</em></span> | |
35 | with <span class="emphasis"><em>x ∈ ℜ</em></span>: | |
36 | </p> | |
37 | <p> | |
38 | <span class="inlinemediaobject"><img src="../../../../../floating_point_eg1.png"></span> | |
39 | </p> | |
40 | <p> | |
41 | Where the difference terms <span class="emphasis"><em>m<sub>n</sub></em></span> are given by: | |
42 | </p> | |
43 | <p> | |
44 | <span class="inlinemediaobject"><img src="../../../../../floating_point_eg2.png"></span> | |
45 | </p> | |
46 | <p> | |
47 | and <span class="emphasis"><em>dx</em></span> is the step-size of the derivative. | |
48 | </p> | |
49 | <p> | |
50 | The third formula in Equation 1 is a three-point central difference rule. | |
51 | It calculates the first derivative of <span class="emphasis"><em>f′(x)</em></span> to <span class="emphasis"><em>O(dx<sup>6</sup>)</em></span>, | |
52 | where <span class="emphasis"><em>dx</em></span> is the given step-size. For example, if | |
53 | the step-size is 0.01 this derivative calculation has about 6 decimal | |
54 | digits of precision - just about right for the 7 decimal digits of single-precision | |
55 | float. Let's make a generic template subroutine using this three-point | |
56 | central difference rule. In particular: | |
57 | </p> | |
58 | <pre class="programlisting"><span class="keyword">template</span><span class="special"><</span><span class="keyword">typename</span> <span class="identifier">value_type</span><span class="special">,</span> <span class="keyword">typename</span> <span class="identifier">function_type</span><span class="special">></span> | |
59 | <span class="identifier">value_type</span> <span class="identifier">derivative</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">x</span><span class="special">,</span> <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">dx</span><span class="special">,</span> <span class="identifier">function_type</span> <span class="identifier">func</span><span class="special">)</span> | |
60 | <span class="special">{</span> | |
61 | <span class="comment">// Compute d/dx[func(*first)] using a three-point</span> | |
62 | <span class="comment">// central difference rule of O(dx^6).</span> | |
63 | ||
64 | <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">dx1</span> <span class="special">=</span> <span class="identifier">dx</span><span class="special">;</span> | |
65 | <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">dx2</span> <span class="special">=</span> <span class="identifier">dx1</span> <span class="special">*</span> <span class="number">2</span><span class="special">;</span> | |
66 | <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">dx3</span> <span class="special">=</span> <span class="identifier">dx1</span> <span class="special">*</span> <span class="number">3</span><span class="special">;</span> | |
67 | ||
68 | <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">m1</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">func</span><span class="special">(</span><span class="identifier">x</span> <span class="special">+</span> <span class="identifier">dx1</span><span class="special">)</span> <span class="special">-</span> <span class="identifier">func</span><span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="identifier">dx1</span><span class="special">))</span> <span class="special">/</span> <span class="number">2</span><span class="special">;</span> | |
69 | <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">m2</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">func</span><span class="special">(</span><span class="identifier">x</span> <span class="special">+</span> <span class="identifier">dx2</span><span class="special">)</span> <span class="special">-</span> <span class="identifier">func</span><span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="identifier">dx2</span><span class="special">))</span> <span class="special">/</span> <span class="number">4</span><span class="special">;</span> | |
70 | <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">m3</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">func</span><span class="special">(</span><span class="identifier">x</span> <span class="special">+</span> <span class="identifier">dx3</span><span class="special">)</span> <span class="special">-</span> <span class="identifier">func</span><span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="identifier">dx3</span><span class="special">))</span> <span class="special">/</span> <span class="number">6</span><span class="special">;</span> | |
71 | ||
72 | <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">fifteen_m1</span> <span class="special">=</span> <span class="number">15</span> <span class="special">*</span> <span class="identifier">m1</span><span class="special">;</span> | |
73 | <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">six_m2</span> <span class="special">=</span> <span class="number">6</span> <span class="special">*</span> <span class="identifier">m2</span><span class="special">;</span> | |
74 | <span class="keyword">const</span> <span class="identifier">value_type</span> <span class="identifier">ten_dx1</span> <span class="special">=</span> <span class="number">10</span> <span class="special">*</span> <span class="identifier">dx1</span><span class="special">;</span> | |
75 | ||
76 | <span class="keyword">return</span> <span class="special">((</span><span class="identifier">fifteen_m1</span> <span class="special">-</span> <span class="identifier">six_m2</span><span class="special">)</span> <span class="special">+</span> <span class="identifier">m3</span><span class="special">)</span> <span class="special">/</span> <span class="identifier">ten_dx1</span><span class="special">;</span> | |
77 | <span class="special">}</span> | |
78 | </pre> | |
79 | <p> | |
80 | The <code class="computeroutput"><span class="identifier">derivative</span><span class="special">()</span></code> | |
81 | template function can be used to compute the first derivative of any | |
82 | function to <span class="emphasis"><em>O(dx<sup>6</sup>)</em></span>. For example, consider the first | |
83 | derivative of <span class="emphasis"><em>sin(x)</em></span> evaluated at <span class="emphasis"><em>x = | |
84 | π/3</em></span>. In other words, | |
85 | </p> | |
86 | <p> | |
87 | <span class="inlinemediaobject"><img src="../../../../../floating_point_eg3.png"></span> | |
88 | </p> | |
89 | <p> | |
90 | The code below computes the derivative in Equation 3 for float, double | |
91 | and boost's multiple-precision type cpp_dec_float_50. | |
92 | </p> | |
93 | <pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">iostream</span><span class="special">></span> | |
94 | <span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">iomanip</span><span class="special">></span> | |
95 | <span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">multiprecision</span><span class="special">/</span><span class="identifier">cpp_dec_float</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> | |
96 | <span class="preprocessor">#include</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">constants</span><span class="special">/</span><span class="identifier">constants</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> | |
97 | ||
98 | ||
99 | <span class="keyword">int</span> <span class="identifier">main</span><span class="special">(</span><span class="keyword">int</span><span class="special">,</span> <span class="keyword">char</span><span class="special">**)</span> | |
100 | <span class="special">{</span> | |
101 | <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">constants</span><span class="special">::</span><span class="identifier">pi</span><span class="special">;</span> | |
102 | <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float_50</span><span class="special">;</span> | |
103 | <span class="comment">//</span> | |
104 | <span class="comment">// We'll pass a function pointer for the function object passed to derivative,</span> | |
105 | <span class="comment">// the typecast is needed to select the correct overload of std::sin:</span> | |
106 | <span class="comment">//</span> | |
107 | <span class="keyword">const</span> <span class="keyword">float</span> <span class="identifier">d_f</span> <span class="special">=</span> <span class="identifier">derivative</span><span class="special">(</span> | |
108 | <span class="identifier">pi</span><span class="special"><</span><span class="keyword">float</span><span class="special">>()</span> <span class="special">/</span> <span class="number">3</span><span class="special">,</span> | |
109 | <span class="number">0.01F</span><span class="special">,</span> | |
110 | <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">float</span><span class="special">(*)(</span><span class="keyword">float</span><span class="special">)>(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">sin</span><span class="special">)</span> | |
111 | <span class="special">);</span> | |
112 | ||
113 | <span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">d_d</span> <span class="special">=</span> <span class="identifier">derivative</span><span class="special">(</span> | |
114 | <span class="identifier">pi</span><span class="special"><</span><span class="keyword">double</span><span class="special">>()</span> <span class="special">/</span> <span class="number">3</span><span class="special">,</span> | |
115 | <span class="number">0.001</span><span class="special">,</span> | |
116 | <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">(*)(</span><span class="keyword">double</span><span class="special">)>(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">sin</span><span class="special">)</span> | |
117 | <span class="special">);</span> | |
118 | <span class="comment">//</span> | |
119 | <span class="comment">// In the cpp_dec_float_50 case, the sin function is multiply overloaded</span> | |
120 | <span class="comment">// to handle expression templates etc. As a result it's hard to take its</span> | |
121 | <span class="comment">// address without knowing about its implementation details. We'll use a </span> | |
122 | <span class="comment">// C++11 lambda expression to capture the call.</span> | |
123 | <span class="comment">// We also need a typecast on the first argument so we don't accidentally pass</span> | |
124 | <span class="comment">// an expression template to a template function:</span> | |
125 | <span class="comment">//</span> | |
126 | <span class="keyword">const</span> <span class="identifier">cpp_dec_float_50</span> <span class="identifier">d_mp</span> <span class="special">=</span> <span class="identifier">derivative</span><span class="special">(</span> | |
127 | <span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="identifier">pi</span><span class="special"><</span><span class="identifier">cpp_dec_float_50</span><span class="special">>()</span> <span class="special">/</span> <span class="number">3</span><span class="special">),</span> | |
128 | <span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="number">1.0E-9</span><span class="special">),</span> | |
129 | <span class="special">[](</span><span class="keyword">const</span> <span class="identifier">cpp_dec_float_50</span><span class="special">&</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">-></span> <span class="identifier">cpp_dec_float_50</span> | |
130 | <span class="special">{</span> | |
131 | <span class="keyword">return</span> <span class="identifier">sin</span><span class="special">(</span><span class="identifier">x</span><span class="special">);</span> | |
132 | <span class="special">}</span> | |
133 | <span class="special">);</span> | |
134 | ||
135 | <span class="comment">// 5.000029e-001</span> | |
136 | <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> | |
137 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">float</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">)</span> | |
138 | <span class="special"><<</span> <span class="identifier">d_f</span> | |
139 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
140 | ||
141 | <span class="comment">// 4.999999999998876e-001</span> | |
142 | <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> | |
143 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">)</span> | |
144 | <span class="special"><<</span> <span class="identifier">d_d</span> | |
145 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
146 | ||
147 | <span class="comment">// 4.99999999999999999999999999999999999999999999999999e-01</span> | |
148 | <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> | |
149 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">cpp_dec_float_50</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">)</span> | |
150 | <span class="special"><<</span> <span class="identifier">d_mp</span> | |
151 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
152 | <span class="special">}</span> | |
153 | </pre> | |
154 | <p> | |
155 | The expected value of the derivative is 0.5. This central difference | |
156 | rule in this example is ill-conditioned, meaning it suffers from slight | |
157 | loss of precision. With that in mind, the results agree with the expected | |
158 | value of 0.5. | |
159 | </p> | |
160 | <p> | |
161 | We can take this a step further and use our derivative function to compute | |
162 | a partial derivative. For example if we take the incomplete gamma function | |
163 | <span class="emphasis"><em>P(a, z)</em></span>, and take the derivative with respect to | |
164 | <span class="emphasis"><em>z</em></span> at <span class="emphasis"><em>(2,2)</em></span> then we can calculate | |
165 | the result as shown below, for good measure we'll compare with the "correct" | |
166 | result obtained from a call to <span class="emphasis"><em>gamma_p_derivative</em></span>, | |
167 | the results agree to approximately 44 digits: | |
168 | </p> | |
169 | <pre class="programlisting"><span class="identifier">cpp_dec_float_50</span> <span class="identifier">gd</span> <span class="special">=</span> <span class="identifier">derivative</span><span class="special">(</span> | |
170 | <span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="number">2</span><span class="special">),</span> | |
171 | <span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="number">1.0E-9</span><span class="special">),</span> | |
172 | <span class="special">[](</span><span class="keyword">const</span> <span class="identifier">cpp_dec_float_50</span><span class="special">&</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">-></span><span class="identifier">cpp_dec_float_50</span> | |
173 | <span class="special">{</span> | |
174 | <span class="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">gamma_p</span><span class="special">(</span><span class="number">2</span><span class="special">,</span> <span class="identifier">x</span><span class="special">);</span> | |
175 | <span class="special">}</span> | |
176 | <span class="special">);</span> | |
177 | <span class="comment">// 2.70670566473225383787998989944968806815263091819151e-01</span> | |
178 | <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> | |
179 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">cpp_dec_float_50</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">)</span> | |
180 | <span class="special"><<</span> <span class="identifier">gd</span> | |
181 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
182 | <span class="comment">// 2.70670566473225383787998989944968806815253190143120e-01</span> | |
183 | <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</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">gamma_p_derivative</span><span class="special">(</span><span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="number">2</span><span class="special">),</span> <span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="number">2</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
184 | </pre> | |
185 | </div> | |
186 | <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> | |
187 | <td align="left"></td> | |
188 | <td align="right"><div class="copyright-footer">Copyright © 2002-2013 John Maddock and Christopher Kormanyos<p> | |
189 | Distributed under the Boost Software License, Version 1.0. (See accompanying | |
190 | 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>) | |
191 | </p> | |
192 | </div></td> | |
193 | </tr></table> | |
194 | <hr> | |
195 | <div class="spirit-nav"> | |
196 | <a accesskey="p" href="jel.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../fp_eg.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="gi.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a> | |
197 | </div> | |
198 | </body> | |
199 | </html> |