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">
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>
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>
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
29 </h5></div></div></div>
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>:
38 <span class=
"inlinemediaobject"><img src=
"../../../../../floating_point_eg1.png"></span>
41 Where the difference terms
<span class=
"emphasis"><em>m
<sub>n
</sub></em></span> are given by:
44 <span class=
"inlinemediaobject"><img src=
"../../../../../floating_point_eg2.png"></span>
47 and
<span class=
"emphasis"><em>dx
</em></span> is the step-size of the derivative.
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:
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>
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>
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>
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>
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>
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,
87 <span class=
"inlinemediaobject"><img src=
"../../../../../floating_point_eg3.png"></span>
90 The code below computes the derivative in Equation
3 for float, double
91 and boost's multiple-precision type cpp_dec_float_50.
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>
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>
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>
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>
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>
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>
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
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:
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>
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>)
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>