123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366 |
- <html>
- <head>
- <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
- <title>Numerical Differentiation</title>
- <link rel="stylesheet" href="../math.css" type="text/css">
- <meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
- <link rel="home" href="../index.html" title="Math Toolkit 2.11.0">
- <link rel="up" href="../quadrature.html" title="Chapter 13. Quadrature and Differentiation">
- <link rel="prev" href="naive_monte_carlo.html" title="Naive Monte Carlo Integration">
- <link rel="next" href="autodiff.html" title="Automatic Differentiation">
- </head>
- <body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
- <table cellpadding="2" width="100%"><tr>
- <td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td>
- <td align="center"><a href="../../../../../index.html">Home</a></td>
- <td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td>
- <td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
- <td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
- <td align="center"><a href="../../../../../more/index.htm">More</a></td>
- </tr></table>
- <hr>
- <div class="spirit-nav">
- <a accesskey="p" href="naive_monte_carlo.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.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="autodiff.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
- </div>
- <div class="section">
- <div class="titlepage"><div><div><h2 class="title" style="clear: both">
- <a name="math_toolkit.diff"></a><a class="link" href="diff.html" title="Numerical Differentiation">Numerical Differentiation</a>
- </h2></div></div></div>
- <h4>
- <a name="math_toolkit.diff.h0"></a>
- <span class="phrase"><a name="math_toolkit.diff.synopsis"></a></span><a class="link" href="diff.html#math_toolkit.diff.synopsis">Synopsis</a>
- </h4>
- <pre class="programlisting"><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">differentiaton</span><span class="special">/</span><span class="identifier">finite_difference</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
- <span class="keyword">namespace</span> <span class="identifier">boost</span> <span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span> <span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">differentiation</span> <span class="special">{</span>
- <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Real</span><span class="special">></span>
- <span class="identifier">Real</span> <span class="identifier">complex_step_derivative</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">x</span><span class="special">);</span>
- <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Real</span><span class="special">,</span> <span class="identifier">size_t</span> <span class="identifier">order</span> <span class="special">=</span> <span class="number">6</span><span class="special">></span>
- <span class="identifier">Real</span> <span class="identifier">finite_difference_derivative</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">x</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">*</span> <span class="identifier">error</span> <span class="special">=</span> <span class="keyword">nullptr</span><span class="special">);</span>
- <span class="special">}}}</span> <span class="comment">// namespaces</span>
- </pre>
- <h4>
- <a name="math_toolkit.diff.h1"></a>
- <span class="phrase"><a name="math_toolkit.diff.description"></a></span><a class="link" href="diff.html#math_toolkit.diff.description">Description</a>
- </h4>
- <p>
- The function <code class="computeroutput"><span class="identifier">finite_difference_derivative</span></code>
- calculates a finite-difference approximation to the derivative of a function
- <span class="emphasis"><em>f</em></span> at point <span class="emphasis"><em>x</em></span>. A basic usage is
- </p>
- <pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">exp</span><span class="special">(</span><span class="identifier">x</span><span class="special">);</span> <span class="special">};</span>
- <span class="keyword">double</span> <span class="identifier">x</span> <span class="special">=</span> <span class="number">1.7</span><span class="special">;</span>
- <span class="keyword">double</span> <span class="identifier">dfdx</span> <span class="special">=</span> <span class="identifier">finite_difference_derivative</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="identifier">x</span><span class="special">);</span>
- </pre>
- <p>
- Finite differencing is complicated, as differentiation is <span class="emphasis"><em>infinitely</em></span>
- ill-conditioned. In addition, for any function implemented in finite-precision
- arithmetic, the "true" derivative is <span class="emphasis"><em>zero</em></span> almost
- everywhere, and undefined at representables. However, some tricks allow for
- reasonable results to be obtained in many cases.
- </p>
- <p>
- There are two sources of error from finite differences: One, the truncation
- error arising from using a finite number of samples to cancel out higher order
- terms in the Taylor series. The second is the roundoff error involved in evaluating
- the function. The truncation error goes to zero as <span class="emphasis"><em>h</em></span> →
- 0, but the roundoff error becomes unbounded. By balancing these two sources
- of error, we can choose a value of <span class="emphasis"><em>h</em></span> that minimizes the
- maximum total error. For this reason boost's <code class="computeroutput"><span class="identifier">finite_difference_derivative</span></code>
- does not require the user to input a stepsize. For more details about the theoretical
- error analysis involved in finite-difference approximations to the derivative,
- see <a href="http://web.archive.org/web/20150420195907/http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf" target="_top">here</a>.
- </p>
- <p>
- Despite the effort that has went into choosing a reasonable value of <span class="emphasis"><em>h</em></span>,
- the problem is still fundamentally ill-conditioned, and hence an error estimate
- is essential. It can be queried as follows
- </p>
- <pre class="programlisting"><span class="keyword">double</span> <span class="identifier">error_estimate</span><span class="special">;</span>
- <span class="keyword">double</span> <span class="identifier">d</span> <span class="special">=</span> <span class="identifier">finite_difference_derivative</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="identifier">x</span><span class="special">,</span> <span class="special">&</span><span class="identifier">error_estimate</span><span class="special">);</span>
- </pre>
- <p>
- N.B.: Producing an error estimate requires additional function evaluations
- and as such is slower than simple evaluation of the derivative. It also expands
- the domain over which the function must be differentiable and requires the
- function to have two more continuous derivatives. The error estimate is computed
- under the assumption that <span class="emphasis"><em>f</em></span> is evaluated to 1ULP. This
- might seem an extreme assumption, but it is the only sensible one, as the routine
- cannot know the functions rounding error. If the function cannot be evaluated
- with very great accuracy, Lanczos's smoothing differentiation is recommended
- as an alternative.
- </p>
- <p>
- The default order of accuracy is 6, which reflects that fact that people tend
- to be interested in functions with many continuous derivatives. If your function
- does not have 7 continuous derivatives, is may be of interest to use a lower
- order method, which can be achieved via (say)
- </p>
- <pre class="programlisting"><span class="keyword">double</span> <span class="identifier">d</span> <span class="special">=</span> <span class="identifier">finite_difference_derivative</span><span class="special"><</span><span class="keyword">decltype</span><span class="special">(</span><span class="identifier">f</span><span class="special">),</span> <span class="identifier">Real</span><span class="special">,</span> <span class="number">2</span><span class="special">>(</span><span class="identifier">f</span><span class="special">,</span> <span class="identifier">x</span><span class="special">);</span>
- </pre>
- <p>
- This requests a second-order accurate derivative be computed.
- </p>
- <p>
- It is emphatically <span class="emphasis"><em>not</em></span> the case that higher order methods
- always give higher accuracy for smooth functions. Higher order methods require
- more addition of positive and negative terms, which can lead to catastrophic
- cancellation. A function which is very good at making a mockery of finite-difference
- differentiation is exp(x)/(cos(x)<sup>3</sup> + sin(x)<sup>3</sup>). Differentiating this function
- by <code class="computeroutput"><span class="identifier">finite_difference_derivative</span></code>
- in double precision at <span class="emphasis"><em>x=5.5</em></span> gives zero correct digits
- at order 4, 6, and 8, but recovers 5 correct digits at order 2. These are dangerous
- waters; use the error estimates to tread carefully.
- </p>
- <p>
- For a finite-difference method of order <span class="emphasis"><em>k</em></span>, the error is
- <span class="emphasis"><em>C</em></span> ε<sup>k/k+1</sup>. In the limit <span class="emphasis"><em>k</em></span> →
- ∞, we see that the error tends to ε, recovering the full precision
- for the type. However, this ignores the fact that higher-order methods require
- subtracting more nearly-equal (perhaps noisy) terms, so the constant <span class="emphasis"><em>C</em></span>
- grows with <span class="emphasis"><em>k</em></span>. Since <span class="emphasis"><em>C</em></span> grows quickly
- and ε<sup>k/k+1</sup> approaches ε slowly, we can see there is a compromise
- between high-order accuracy and conditioning of the difference quotient. In
- practice we have found that <span class="emphasis"><em>k=6</em></span> seems to be a good compromise
- between the two (and have made this the default), but users are encouraged
- to examine the error estimates to choose an optimal order of accuracy for the
- given problem.
- </p>
- <div class="table">
- <a name="math_toolkit.diff.id"></a><p class="title"><b>Table 13.1. Cost of Finite-Difference Numerical Differentiation</b></p>
- <div class="table-contents"><table class="table" summary="Cost of Finite-Difference Numerical Differentiation">
- <colgroup>
- <col>
- <col>
- <col>
- <col>
- <col>
- </colgroup>
- <thead><tr>
- <th>
- <p>
- Order of Accuracy
- </p>
- </th>
- <th>
- <p>
- Function Evaluations
- </p>
- </th>
- <th>
- <p>
- Error
- </p>
- </th>
- <th>
- <p>
- Continuous Derivatives Required for Error Estimate to Hold
- </p>
- </th>
- <th>
- <p>
- Additional Function Evaluations to Produce Error Estimates
- </p>
- </th>
- </tr></thead>
- <tbody>
- <tr>
- <td>
- <p>
- 1
- </p>
- </td>
- <td>
- <p>
- 2
- </p>
- </td>
- <td>
- <p>
- ε<sup>1/2</sup>
- </p>
- </td>
- <td>
- <p>
- 2
- </p>
- </td>
- <td>
- <p>
- 1
- </p>
- </td>
- </tr>
- <tr>
- <td>
- <p>
- 2
- </p>
- </td>
- <td>
- <p>
- 2
- </p>
- </td>
- <td>
- <p>
- ε<sup>2/3</sup>
- </p>
- </td>
- <td>
- <p>
- 3
- </p>
- </td>
- <td>
- <p>
- 2
- </p>
- </td>
- </tr>
- <tr>
- <td>
- <p>
- 4
- </p>
- </td>
- <td>
- <p>
- 4
- </p>
- </td>
- <td>
- <p>
- ε<sup>4/5</sup>
- </p>
- </td>
- <td>
- <p>
- 5
- </p>
- </td>
- <td>
- <p>
- 2
- </p>
- </td>
- </tr>
- <tr>
- <td>
- <p>
- 6
- </p>
- </td>
- <td>
- <p>
- 6
- </p>
- </td>
- <td>
- <p>
- ε<sup>6/7</sup>
- </p>
- </td>
- <td>
- <p>
- 7
- </p>
- </td>
- <td>
- <p>
- 2
- </p>
- </td>
- </tr>
- <tr>
- <td>
- <p>
- 8
- </p>
- </td>
- <td>
- <p>
- 8
- </p>
- </td>
- <td>
- <p>
- ε<sup>8/9</sup>
- </p>
- </td>
- <td>
- <p>
- 9
- </p>
- </td>
- <td>
- <p>
- 2
- </p>
- </td>
- </tr>
- </tbody>
- </table></div>
- </div>
- <br class="table-break"><p>
- Given all the caveats which must be kept in mind for successful use of finite-difference
- differentiation, it is reasonable to try to avoid it if possible. Boost provides
- two possibilities: The Chebyshev transform (see <a class="link" href="sf_poly/chebyshev.html" title="Chebyshev Polynomials">here</a>)
- and the complex step derivative. If your function is the restriction to the
- real line of a holomorphic function which takes real values at real argument,
- then the <span class="bold"><strong>complex step derivative</strong></span> can be used.
- The idea is very simple: Since <span class="emphasis"><em>f</em></span> is complex-differentiable,
- <span class="emphasis"><em>f(x+ⅈ h) = f(x) + ⅈ hf'(x) - h<sup>2</sup>f''(x) + 𝑶(h<sup>3</sup>)</em></span>.
- As long as <span class="emphasis"><em>f(x)</em></span> ∈ ℝ, then <span class="emphasis"><em>f'(x)
- = ℑ f(x+ⅈ h)/h + 𝑶(h<sup>2</sup>)</em></span>. This method requires a single
- complex function evaluation and is not subject to the catastrophic subtractive
- cancellation that plagues finite-difference calculations.
- </p>
- <p>
- An example usage:
- </p>
- <pre class="programlisting"><span class="keyword">double</span> <span class="identifier">x</span> <span class="special">=</span> <span class="number">7.2</span><span class="special">;</span>
- <span class="keyword">double</span> <span class="identifier">e_prime</span> <span class="special">=</span> <span class="identifier">complex_step_derivative</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">exp</span><span class="special"><</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="keyword">double</span><span class="special">>>,</span> <span class="identifier">x</span><span class="special">);</span>
- </pre>
- <p>
- References:
- </p>
- <div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
- <li class="listitem">
- Squire, William, and George Trapp. <span class="emphasis"><em>Using complex variables to
- estimate derivatives of real functions.</em></span> Siam Review 40.1 (1998):
- 110-112.
- </li>
- <li class="listitem">
- Fornberg, Bengt. <span class="emphasis"><em>Generation of finite difference formulas on
- arbitrarily spaced grids.</em></span> Mathematics of computation 51.184
- (1988): 699-706.
- </li>
- <li class="listitem">
- Corless, Robert M., and Nicolas Fillion. <span class="emphasis"><em>A graduate introduction
- to numerical methods.</em></span> AMC 10 (2013): 12.
- </li>
- </ul></div>
- </div>
- <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
- <td align="left"></td>
- <td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar
- Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
- Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
- Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg,
- Daryle Walker and Xiaogang Zhang<p>
- Distributed under the Boost Software License, Version 1.0. (See accompanying
- 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>)
- </p>
- </div></td>
- </tr></table>
- <hr>
- <div class="spirit-nav">
- <a accesskey="p" href="naive_monte_carlo.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.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="autodiff.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
- </div>
- </body>
- </html>
|