mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
159 lines
20 KiB
HTML
159 lines
20 KiB
HTML
<html>
|
||
<head>
|
||
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
|
||
<title>Roots of Cubic Polynomials</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 4.2.1">
|
||
<link rel="up" href="../root_finding.html" title="Chapter 10. Root Finding & Minimization Algorithms">
|
||
<link rel="prev" href="roots_deriv.html" title="Root Finding With Derivatives: Newton-Raphson, Halley & Schröder">
|
||
<link rel="next" href="quartic_roots.html" title="Roots of Quartic Polynomials">
|
||
<meta name="viewport" content="width=device-width, initial-scale=1">
|
||
</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="roots_deriv.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding.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="quartic_roots.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.cubic_roots"></a><a class="link" href="cubic_roots.html" title="Roots of Cubic Polynomials">Roots of Cubic Polynomials</a>
|
||
</h2></div></div></div>
|
||
<h4>
|
||
<a name="math_toolkit.cubic_roots.h0"></a>
|
||
<span class="phrase"><a name="math_toolkit.cubic_roots.synopsis"></a></span><a class="link" href="cubic_roots.html#math_toolkit.cubic_roots.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">tools</span><span class="special">/</span><span class="identifier">cubic_roots</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="identifier">math</span><span class="special">::</span><span class="identifier">tools</span> <span class="special">{</span>
|
||
|
||
<span class="comment">// Solves ax³ + bx² + cx + d = 0.</span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span><span class="identifier">Real</span><span class="special">,</span><span class="number">3</span><span class="special">></span> <span class="identifier">cubic_roots</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">c</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">d</span><span class="special">);</span>
|
||
|
||
<span class="comment">// Computes the empirical residual p(r) (first element) and expected residual ε|rṗ(r)| (second element) for a root.</span>
|
||
<span class="comment">// Recall that for a numerically computed root r satisfying r = r⁎(1+ε) for the exact root r⁎ of a function p, |p(r)| ≲ ε|rṗ(r)|.</span>
|
||
<span class="keyword">template</span><span class="special"><</span><span class="keyword">typename</span> <span class="identifier">Real</span><span class="special">></span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">array</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">cubic_root_residual</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">c</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">d</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">root</span><span class="special">);</span>
|
||
|
||
<span class="comment">// Computes the condition number of rootfinding. Computed via Corless, A Graduate Introduction to Numerical Methods, Section 3.2.1:</span>
|
||
<span class="keyword">template</span><span class="special"><</span><span class="keyword">typename</span> <span class="identifier">Real</span><span class="special">></span>
|
||
<span class="identifier">Real</span> <span class="identifier">cubic_root_condition_number</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">c</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">d</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">root</span><span class="special">);</span>
|
||
<span class="special">}</span>
|
||
</pre>
|
||
<h4>
|
||
<a name="math_toolkit.cubic_roots.h1"></a>
|
||
<span class="phrase"><a name="math_toolkit.cubic_roots.background"></a></span><a class="link" href="cubic_roots.html#math_toolkit.cubic_roots.background">Background</a>
|
||
</h4>
|
||
<p>
|
||
The <code class="computeroutput"><span class="identifier">cubic_roots</span></code> function extracts
|
||
all real roots of a cubic polynomial ax³ + bx² + cx + d. The result is a
|
||
<code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span><span class="identifier">Real</span><span class="special">,</span> <span class="number">3</span><span class="special">></span></code>, which has length three, irrespective of
|
||
whether there are 3 real roots. There is always 1 real root, and hence (barring
|
||
overflow or other exceptional circumstances), the first element of the <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">array</span></code>
|
||
is always populated. If there is only one real root of the polynomial, then
|
||
the second and third elements are set to <code class="computeroutput"><span class="identifier">nans</span></code>.
|
||
The roots are returned in nondecreasing order.
|
||
</p>
|
||
<p>
|
||
Be careful with double roots. First, if you have a real double root, it is
|
||
numerically indistinguishable from a complex conjugate pair of roots, where
|
||
the complex part is tiny. Second, the condition number of rootfinding is infinite
|
||
at a double root, so even changes as subtle as different instruction generation
|
||
can change the outcome. We have some heuristics in place to detect double roots,
|
||
but these should be regarded with suspicion.
|
||
</p>
|
||
<h4>
|
||
<a name="math_toolkit.cubic_roots.h2"></a>
|
||
<span class="phrase"><a name="math_toolkit.cubic_roots.example"></a></span><a class="link" href="cubic_roots.html#math_toolkit.cubic_roots.example">Example</a>
|
||
</h4>
|
||
<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">iostream</span><span class="special">></span>
|
||
<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">sstream</span><span class="special">></span>
|
||
<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">tools</span><span class="special">/</span><span class="identifier">cubic_roots</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
|
||
|
||
<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">tools</span><span class="special">::</span><span class="identifier">cubic_roots</span><span class="special">;</span>
|
||
<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">tools</span><span class="special">::</span><span class="identifier">cubic_root_residual</span><span class="special">;</span>
|
||
|
||
<span class="keyword">template</span><span class="special"><</span><span class="keyword">typename</span> <span class="identifier">Real</span><span class="special">></span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">string</span> <span class="identifier">print_roots</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span><span class="identifier">Real</span><span class="special">,</span> <span class="number">3</span><span class="special">></span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">roots</span><span class="special">)</span> <span class="special">{</span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">ostringstream</span> <span class="identifier">out</span><span class="special">;</span>
|
||
<span class="identifier">out</span> <span class="special"><<</span> <span class="string">"{"</span> <span class="special"><<</span> <span class="identifier">roots</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special"><<</span> <span class="string">", "</span> <span class="special"><<</span> <span class="identifier">roots</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special"><<</span> <span class="string">", "</span> <span class="special"><<</span> <span class="identifier">roots</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special"><<</span> <span class="string">"}"</span><span class="special">;</span>
|
||
<span class="keyword">return</span> <span class="identifier">out</span><span class="special">.</span><span class="identifier">str</span><span class="special">();</span>
|
||
<span class="special">}</span>
|
||
|
||
<span class="keyword">int</span> <span class="identifier">main</span><span class="special">()</span> <span class="special">{</span>
|
||
<span class="comment">// Solves x³ - 6x² + 11x - 6 = (x-1)(x-2)(x-3).</span>
|
||
<span class="keyword">auto</span> <span class="identifier">roots</span> <span class="special">=</span> <span class="identifier">cubic_roots</span><span class="special">(</span><span class="number">1.0</span><span class="special">,</span> <span class="special">-</span><span class="number">6.0</span><span class="special">,</span> <span class="number">11.0</span><span class="special">,</span> <span class="special">-</span><span class="number">6.0</span><span class="special">);</span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"The roots of x³ - 6x² + 11x - 6 are "</span> <span class="special"><<</span> <span class="identifier">print_roots</span><span class="special">(</span><span class="identifier">roots</span><span class="special">)</span> <span class="special"><<</span> <span class="string">".\n"</span><span class="special">;</span>
|
||
|
||
<span class="comment">// Double root; YMMV:</span>
|
||
<span class="comment">// (x+1)²(x-2) = x³ - 3x - 2:</span>
|
||
<span class="identifier">roots</span> <span class="special">=</span> <span class="identifier">cubic_roots</span><span class="special">(</span><span class="number">1.0</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="special">-</span><span class="number">3.0</span><span class="special">,</span> <span class="special">-</span><span class="number">2.0</span><span class="special">);</span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"The roots of x³ - 3x - 2 are "</span> <span class="special"><<</span> <span class="identifier">print_roots</span><span class="special">(</span><span class="identifier">roots</span><span class="special">)</span> <span class="special"><<</span> <span class="string">".\n"</span><span class="special">;</span>
|
||
|
||
<span class="comment">// Single root: (x-i)(x+i)(x-3) = x³ - 3x² + x - 3:</span>
|
||
<span class="identifier">roots</span> <span class="special">=</span> <span class="identifier">cubic_roots</span><span class="special">(</span><span class="number">1.0</span><span class="special">,</span> <span class="special">-</span><span class="number">3.0</span><span class="special">,</span> <span class="number">1.0</span><span class="special">,</span> <span class="special">-</span><span class="number">3.0</span><span class="special">);</span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"The real roots of x³ - 3x² + x - 3 are "</span> <span class="special"><<</span> <span class="identifier">print_roots</span><span class="special">(</span><span class="identifier">roots</span><span class="special">)</span> <span class="special"><<</span> <span class="string">".\n"</span><span class="special">;</span>
|
||
|
||
<span class="comment">// I don't know the roots of x³ + 6.28x² + 2.3x + 3.6;</span>
|
||
<span class="comment">// how can I see if they've been computed correctly?</span>
|
||
<span class="identifier">roots</span> <span class="special">=</span> <span class="identifier">cubic_roots</span><span class="special">(</span><span class="number">1.0</span><span class="special">,</span> <span class="number">6.28</span><span class="special">,</span> <span class="number">2.3</span><span class="special">,</span> <span class="number">3.6</span><span class="special">);</span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"The real root of x³ +6.28x² + 2.3x + 3.6 is "</span> <span class="special"><<</span> <span class="identifier">roots</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special"><<</span> <span class="string">".\n"</span><span class="special">;</span>
|
||
<span class="keyword">auto</span> <span class="identifier">res</span> <span class="special">=</span> <span class="identifier">cubic_root_residual</span><span class="special">(</span><span class="number">1.0</span><span class="special">,</span> <span class="number">6.28</span><span class="special">,</span> <span class="number">2.3</span><span class="special">,</span> <span class="number">3.6</span><span class="special">,</span> <span class="identifier">roots</span><span class="special">[</span><span class="number">0</span><span class="special">]);</span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"The residual is "</span> <span class="special"><<</span> <span class="identifier">res</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special"><<</span> <span class="string">", and the expected residual is "</span> <span class="special"><<</span> <span class="identifier">res</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special"><<</span> <span class="string">". "</span><span class="special">;</span>
|
||
<span class="keyword">if</span> <span class="special">(</span><span class="identifier">abs</span><span class="special">(</span><span class="identifier">res</span><span class="special">[</span><span class="number">0</span><span class="special">])</span> <span class="special"><=</span> <span class="identifier">res</span><span class="special">[</span><span class="number">1</span><span class="special">])</span> <span class="special">{</span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"This is an expected accuracy.\n"</span><span class="special">;</span>
|
||
<span class="special">}</span> <span class="keyword">else</span> <span class="special">{</span>
|
||
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cerr</span> <span class="special"><<</span> <span class="string">"The residual is unexpectedly large.\n"</span><span class="special">;</span>
|
||
<span class="special">}</span>
|
||
<span class="special">}</span>
|
||
</pre>
|
||
<p>
|
||
This prints:
|
||
</p>
|
||
<pre class="programlisting">The roots of x<sup>3</sup> - 6x<sup>2</sup> + 11x - 6 are {1, 2, 3}.
|
||
The roots of x<sup>3</sup> - 3x - 2 are {-1, -1, 2}.
|
||
The real roots of x<sup>3</sup> - 3x<sup>2</sup> + x - 3 are {3, nan, nan}.
|
||
The real root of x<sup>3</sup> +6.28x<sup>2</sup> + 2.3x + 3.6 is -5.99656.
|
||
The residual is -1.56586e-14, and the expected residual is 4.64155e-14. This is an expected accuracy.
|
||
</pre>
|
||
<h4>
|
||
<a name="math_toolkit.cubic_roots.h3"></a>
|
||
<span class="phrase"><a name="math_toolkit.cubic_roots.performance_and_accuracy"></a></span><a class="link" href="cubic_roots.html#math_toolkit.cubic_roots.performance_and_accuracy">Performance
|
||
and Accuracy</a>
|
||
</h4>
|
||
<p>
|
||
On an Intel laptop chip running at 2700 MHz, we observe 3 roots taking ~90ns
|
||
to compute. If the polynomial only possesses a single real root, it takes ~50ns.
|
||
See <code class="computeroutput"><span class="identifier">reporting</span><span class="special">/</span><span class="identifier">performance</span><span class="special">/</span><span class="identifier">cubic_roots_performance</span><span class="special">.</span><span class="identifier">cpp</span></code>.
|
||
</p>
|
||
<p>
|
||
The forward error cannot be effectively bounded. However, for an approximate
|
||
root r returned by this routine, the residuals should be constrained by |p(r)|
|
||
≲ ε|rṗ(r)|, where '≲' should be interpreted as an order of magnitude
|
||
estimate.
|
||
</p>
|
||
</div>
|
||
<div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
|
||
Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
|
||
Lalande, John Maddock, Evan Miller, 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>
|
||
<hr>
|
||
<div class="spirit-nav">
|
||
<a accesskey="p" href="roots_deriv.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding.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="quartic_roots.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
|
||
</div>
|
||
</body>
|
||
</html>
|