boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Complex Newton's method. Zero derivatives handled by Muller's Method.

Browse Source
This commit is contained in:
Nick Thompson
2018-11-25 00:19:11 -07:00
parent 51107199b5
commit 3ab69d00ee
4 changed files with 177 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
doc/roots/roots.qbk
View File

@@ -29,6 +29,9 @@
template <class F, class T>
T schroder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);
template <class F, class Complex>
Complex complex_newton(F f, Complex guess, int max_iterations = std::numeric_limits<typename Complex::value_type>::digits);
}}} // namespaces boost::math::tools.
[h4 Description]
@@ -40,9 +43,10 @@ These functions all perform iterative root-finding [*using derivatives]:
* `halley_iterate` and `schroder_iterate` perform third-order
__halley and __schroder iteration.
The functions all take the same parameters:
* `complex_newton` performs Newton's method on complex-analytic functions.
[variablelist Parameters of the root finding functions
[variablelist Parameters of the real-valued root finding functions
[[F f] [Type F must be a callable function object (or C++ lambda) that accepts one parameter and
returns a __tuple_type:
@@ -158,6 +162,15 @@ College Park, MD: University of Maryland, Institute for Advanced Computer Studie
This method guarantees at least quadratic convergence (the same as Newton's method), and is known to work well in the presence of multiple roots:
something that neither Newton nor Halley can do.
The complex Newton method works slightly differently than the rest of the methods:
Since there is no way to bracket roots in the complex plane,
the `min` and `max` arguments are not accepted.
Failure to reach a root is communicated by returning `nan`s.
Remember that if a function has many roots,
then which root the complex Newton's method converges to is essentially impossible to predict a priori; see the Newton's fractal for more information.
An example usage of `complex_newton` is given in `examples/daubechies_coefficients.cpp`.
[h4 Examples]
See __root_finding_examples.

45
example/daubechies_coefficients.cpp
View File

@@ -24,6 +24,7 @@ using boost::math::binomial_coefficient;
using boost::math::tools::schroder_iterate;
using boost::math::tools::halley_iterate;
using boost::math::tools::newton_raphson_iterate;
using boost::math::tools::complex_newton;
using boost::math::constants::half;
using boost::math::constants::root_two;
using boost::math::constants::pi;
@@ -31,37 +32,6 @@ using boost::math::quadrature::gauss_kronrod;
using boost::multiprecision::cpp_bin_float_100;
using boost::multiprecision::cpp_complex_100;
template<class Complex, class F>
Complex complex_newton(F g, Complex guess, typename Complex::value_type search_radius=100, int iterations=30)
{
typedef typename Complex::value_type Real;
Complex xn = guess;
Complex f;
Complex df;
do {
auto pair = g(xn);
f = pair.first;
df = pair.second;
--iterations;
xn = xn - (f/df);
if (abs(f) < sqrt(std::numeric_limits<Real>::epsilon()))
{
return xn;
}
} while(iterations);
if (abs(f) > std::numeric_limits<typename Complex::value_type>::epsilon())
{
return {std::numeric_limits<typename Complex::value_type>::quiet_NaN(), std::numeric_limits<typename Complex::value_type>::quiet_NaN()};
}
return xn;
}
template<class Complex>
std::vector<std::pair<Complex, Complex>> find_roots(size_t p)
{
@@ -76,6 +46,8 @@ std::vector<std::pair<Complex, Complex>> find_roots(size_t p)
polynomial<Complex> P(std::move(coeffs));
polynomial<Complex> Pcopy = P;
polynomial<Complex> Pcopy_prime = P.prime();
auto orig = [&](Complex z) { return std::make_pair<Complex, Complex>(Pcopy(z), Pcopy_prime(z)); };
polynomial<Complex> P_prime = P.prime();
@@ -94,7 +66,11 @@ std::vector<std::pair<Complex, Complex>> find_roots(size_t p)
return std::make_pair<Complex, Complex>(P(x), P_prime(x));
};
Complex r = complex_newton(f, guess, search_radius, 100);
Complex r = complex_newton(f, guess);
// Refine r with the original function.
// We only use the polynomial division to ensure we don't get the same root over and over.
// However, the division induces error which can grow quickly-or slowly! See Numerical Recipes, section 9.5.1.
r = complex_newton(orig, r);
// Test the root:
using std::sqrt;
Real tol = sqrt(sqrt(std::numeric_limits<Real>::epsilon()));
@@ -210,8 +186,7 @@ std::vector<typename Complex::value_type> daubechies_coefficients(std::vector<st
int main()
{
typedef cpp_complex_100 Complex;
// This is probably the most ill-conditioned algorithm I've ever written.
for(size_t p = 1; p < 57; ++p)
for(size_t p = 1; p < 89; ++p)
{
auto roots = find_roots<Complex>(p);
auto h = daubechies_coefficients(roots);
@@ -219,6 +194,6 @@ int main()
for (auto& x : h) {
std::cout << x << ", ";
}
std::cout << "}\n";
std::cout << "}\n\n\n\n";
}
}

81
include/boost/math/tools/roots.hpp
View File

@@ -9,7 +9,7 @@
#ifdef _MSC_VER
#pragma once
#endif
#include <iostream>
#include <utility>
#include <boost/config/no_tr1/cmath.hpp>
#include <stdexcept>
@@ -65,7 +65,7 @@ inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T)
{
using dummy::get;
// Rely on ADL to find the correct overload of get:
val = get<0>(t);
val = get<0>(t);
}
template <class T, class U, class V>
@@ -559,10 +559,85 @@ inline T schroeder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEP
return schroder_iterate(f, guess, min, max, digits, m);
}
#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS
/* Why do we set the default maximum number of iterations to the number of digits in the type?
* Because for double roots, the number of digits increases linearly with the number of iterations, so this default should recover full precision.
* For isolated roots, the problem is so rapidly convergent that this doesn't matter at all.
* For nonexistence, the problem is harder.
*/
template<class Complex, class F>
Complex complex_newton(F g, Complex guess, int max_iterations=std::numeric_limits<typename Complex::value_type>::digits)
{
typedef typename Complex::value_type Real;
using std::norm;
Complex z0 = guess/static_cast<Complex>(4);
Complex z1 = guess/static_cast<Complex>(2);
Complex z2 = guess;
do {
auto pair = g(z2);
if (abs(pair.first) < 10*std::numeric_limits<Real>::epsilon())
{
// Computed it, might as well use it:
if (abs(pair.second) > sqrt(std::numeric_limits<Real>::epsilon()))
{
return z2 - pair.first/pair.second;
}
return z2;
}
if (norm(pair.second) == 0)
{
// Muller's method. Notation follows Numerical Recipes, 9.5.2:
Complex q = (z2 - z1)/(z1 - z0);
auto P0 = g(z0);
auto P1 = g(z1);
Complex qp1 = static_cast<Complex>(1)+q;
Complex A = q*(pair.first - qp1*P1.first + q*P0.first);
Complex B = (static_cast<Complex>(2)*q+static_cast<Complex>(1))*pair.first - qp1*qp1*P1.first +q*q*P0.first;
Complex C = qp1*pair.first;
Complex rad = sqrt(B*B - static_cast<Complex>(4)*A*C);
Complex denom1 = B + rad;
Complex denom2 = B - rad;
Complex correction = (z1-z2)*static_cast<Complex>(2)*C;
if (norm(denom1) > norm(denom2))
{
correction /= denom1;
} else {
correction /= denom2;
}
z0 = z1;
z1 = z2;
z2 = z2 + correction;
}
else
{
z0 = z1;
z1 = z2;
z2 = z2 - (pair.first/pair.second);
}
} while(max_iterations--);
// At this point, maybe we should use complex trapezoidal quadrature to see if there are any roots?
// https://www.ams.org/journals/mcom/1967-21-100/S0025-5718-1967-0228165-4/S0025-5718-1967-0228165-4.pdf
// Exposition of the above: http://www.chebfun.org/examples/roots/ComplexRoots.html
// The idea is that if we can get abs(f) < eps, we should, but if we go through all these iterations
// and abs(f) < sqrt(eps), then roundoff error simply does not allow that we can evaluate f to < esp
auto pair = g(z2);
if (abs(pair.first) < sqrt(std::numeric_limits<Real>::epsilon()))
{
return z2;
}
return {std::numeric_limits<Real>::quiet_NaN(),
std::numeric_limits<Real>::quiet_NaN()};
}
#endif
} // namespace tools
} // namespace math
} // namespace boost
#endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP

78
test/test_roots.cpp
View File

@@ -11,14 +11,23 @@
#include <boost/test/results_collector.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <boost/math/distributions/skew_normal.hpp>
#include <boost/math/tools/polynomial.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/test/results_collector.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/array.hpp>
#include <boost/type_index.hpp>
#include "table_type.hpp"
#include <iostream>
#include <iomanip>
#if BOOST_HAS_FLOAT128
#include <boost/multiprecision/complex128.hpp>
#endif
#include <boost/multiprecision/cpp_complex.hpp>
#define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
{\
unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
@@ -317,13 +326,74 @@ void test_beta(T, const char* /* name */)
test_inverses<T>(ibeta_large_data);
}
#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS
template <class Complex>
void test_complex_newton()
{
typedef typename Complex::value_type Real;
std::cout << "Testing complex Newton's Method on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
using boost::math::tools::complex_newton;
using boost::math::tools::polynomial;
using boost::math::constants::half;
Real tol = std::numeric_limits<Real>::epsilon();
// p(z) = z^2 + 1, roots: \pm i.
polynomial<Complex> p{{1,0}, {0, 0}, {1,0}};
Complex guess{1,1};
polynomial<Complex> p_prime = p.prime();
auto f = [&](Complex z) { return std::make_pair<Complex, Complex>(p(z), p_prime(z)); };
Complex root = complex_newton(f, guess);
BOOST_CHECK(abs(root.real()) <= tol);
BOOST_CHECK_CLOSE(root.imag(), 1, tol);
guess = -guess;
root = complex_newton(f, guess);
BOOST_CHECK(abs(root.real()) <= tol);
BOOST_CHECK_CLOSE(root.imag(), -1, tol);
// Test that double roots are handled correctly-as correctly as possible.
// Convergence at a double root is not quadratic.
// This sets p = (z-i)^2:
p = polynomial<Complex>({{-1,0}, {0,-2}, {1,0}});
p_prime = p.prime();
guess = -guess;
auto g = [&](Complex z) { return std::make_pair<Complex, Complex>(p(z), p_prime(z)); };
root = complex_newton(g, guess);
BOOST_CHECK(abs(root.real()) < 10*sqrt(tol));
BOOST_CHECK_CLOSE(root.imag(), 1, tol);
// Test that zero derivatives are handled.
// p(z) = z^2 + iz + 1
p = polynomial<Complex>({{1,0}, {0,1}, {1,0}});
// p'(z) = 2z + i
p_prime = p.prime();
guess = Complex(0,-boost::math::constants::half<Real>());
auto g2 = [&](Complex z) { return std::make_pair<Complex, Complex>(p(z), p_prime(z)); };
root = complex_newton(g2, guess);
// Here's the other root, in case code changes cause it to be found:
//Complex expected_root1{0, half<Real>()*(sqrt(static_cast<Real>(5)) - static_cast<Real>(1))};
Complex expected_root2{0, -half<Real>()*(sqrt(static_cast<Real>(5)) + static_cast<Real>(1))};
BOOST_CHECK_CLOSE(expected_root2.imag(),root.imag(), tol);
BOOST_CHECK(abs(root.real()) < tol);
}
#endif
BOOST_AUTO_TEST_CASE( test_main )
{
test_beta(0.1, "double");
// bug reports:
boost::math::skew_normal_distribution<> dist(2.0, 1.0, -2.5);
BOOST_CHECK(boost::math::isfinite(quantile(dist, 0.075)));
#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS
test_complex_newton<std::complex<float>>();
test_complex_newton<std::complex<double>>();
test_complex_newton<std::complex<long double>>();
test_complex_newton<boost::multiprecision::cpp_complex_100>();
#if BOOST_HAS_FLOAT128
test_complex_newton<boost::multiprecision::complex128>();
#endif
#endif
}