libs/multiprecision/doc/html/boost_multiprecision/tut/floats/fp_eg/nd.html

   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="../../../../multiprecision.css" type="text/css">
   6 <meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
   7 <link rel="home" href="../../../../index.html" title="Chapter&#160;1.&#160;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&#8242;(x)</em></span>
  35             with <span class="emphasis"><em>x &#8712; &#8476;</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&#8242;(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">&lt;</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">&gt;</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             &#960;/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">&lt;</span><span class="identifier">iostream</span><span class="special">&gt;</span>
  94 <span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">iomanip</span><span class="special">&gt;</span>
  95 <span class="preprocessor">#include</span> <span class="special">&lt;</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">&gt;</span>
  96 <span class="preprocessor">#include</span> <span class="special">&lt;</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">&gt;</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">&lt;</span><span class="keyword">float</span><span class="special">&gt;()</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">&lt;</span><span class="keyword">float</span><span class="special">(*)(</span><span class="keyword">float</span><span class="special">)&gt;(</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;()</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">&lt;</span><span class="keyword">double</span><span class="special">(*)(</span><span class="keyword">double</span><span class="special">)&gt;(</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">&lt;</span><span class="identifier">cpp_dec_float_50</span><span class="special">&gt;()</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">&amp;</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">-&gt;</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">&lt;&lt;</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">&lt;</span><span class="keyword">float</span><span class="special">&gt;::</span><span class="identifier">digits10</span><span class="special">)</span>
 138       <span class="special">&lt;&lt;</span> <span class="identifier">d_f</span>
 139       <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</span><span class="identifier">digits10</span><span class="special">)</span>
 144       <span class="special">&lt;&lt;</span> <span class="identifier">d_d</span>
 145       <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;</span><span class="identifier">cpp_dec_float_50</span><span class="special">&gt;::</span><span class="identifier">digits10</span><span class="special">)</span>
 150       <span class="special">&lt;&lt;</span> <span class="identifier">d_mp</span>
 151       <span class="special">&lt;&lt;</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">&amp;</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">-&gt;</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">&lt;&lt;</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">&lt;</span><span class="identifier">cpp_dec_float_50</span><span class="special">&gt;::</span><span class="identifier">digits10</span><span class="special">)</span>
 180    <span class="special">&lt;&lt;</span> <span class="identifier">gd</span>
 181    <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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 &#169; 2002-2019 John Maddock
 189       and Christopher Kormanyos<p>
 190         Distributed under the Boost Software License, Version 1.0. (See accompanying
 191         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>)
 192       </p>
 193 </div></td>
 194 </tr></table>
 195 <hr>
 196 <div class="spirit-nav">
 197 <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>
 198 </div>
 199 </body>
 200 </html>