From 56a2002848a69bc1a5d634d1e6fe3301f41d7a24 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 11 Jun 2021 11:30:12 +0200
Subject: [PATCH] Add root computer for 1d polynomials
By now, use a simple bisection method.
Only compute the first root (if it exists in the given interval).
Would be easy to get the full list of roots in the interval if needed.
---
src/analysis/CMakeLists.txt | 1 +
src/analysis/Polynomial1DRootComputer.cpp | 64 +++++++++++++++++++++++
src/analysis/Polynomial1DRootComputer.hpp | 31 +++++++++++
tests/CMakeLists.txt | 1 +
tests/test_Polynomial1DRootComputer.cpp | 45 ++++++++++++++++
5 files changed, 142 insertions(+)
create mode 100644 src/analysis/Polynomial1DRootComputer.cpp
create mode 100644 src/analysis/Polynomial1DRootComputer.hpp
create mode 100644 tests/test_Polynomial1DRootComputer.cpp
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
index 669af876e..fda3d1b2c 100644
--- a/src/analysis/CMakeLists.txt
+++ b/src/analysis/CMakeLists.txt
@@ -2,5 +2,6 @@
add_library(
PugsAnalysis
+ Polynomial1DRootComputer.cpp
SturmSequence.cpp
)
diff --git a/src/analysis/Polynomial1DRootComputer.cpp b/src/analysis/Polynomial1DRootComputer.cpp
new file mode 100644
index 000000000..c86d2ad31
--- /dev/null
+++ b/src/analysis/Polynomial1DRootComputer.cpp
@@ -0,0 +1,64 @@
+#include <analysis/Polynomial1DRootComputer.hpp>
+
+double
+Polynomial1DRootComputer::_computeRootBetween(const PositionSignChanges& left_bound,
+ const PositionSignChanges& right_bound) const
+{
+ double left = left_bound.position;
+ double right = right_bound.position;
+
+ const double interval_size = right - left;
+
+ const Polynomial1D& p0 = m_sturm_sequence.p(0);
+ double left_value = p0(left);
+
+ int cpt = 0;
+ while (right - left > 1E-12 * interval_size) {
+ double middle = 0.5 * (right + left);
+ double middle_value = p0(middle);
+ if (middle_value * left_value < 0) {
+ right = middle;
+ } else {
+ left_value = middle_value;
+ left = middle;
+ }
+ cpt++;
+ }
+
+ return 0.5 * (left + right);
+}
+
+std::optional<double>
+Polynomial1DRootComputer::getFirstRoot() const
+{
+ PositionSignChanges l = {m_left, m_sturm_sequence.numberOfSignChanges(m_left)};
+ PositionSignChanges r = {m_right, m_sturm_sequence.numberOfSignChanges(m_right)};
+
+ if (l.nb_sign_changes - r.nb_sign_changes > 0) {
+ if (l.nb_sign_changes - r.nb_sign_changes > 1) {
+ PositionSignChanges a = l;
+ PositionSignChanges b = r;
+ do {
+ Assert(a.nb_sign_changes - b.nb_sign_changes > 0);
+
+ double position_k = 0.5 * (a.position + b.position);
+ PositionSignChanges rk = {position_k, m_sturm_sequence.numberOfSignChanges(position_k)};
+
+ if (l.nb_sign_changes - rk.nb_sign_changes > 0) {
+ b = rk;
+ } else if (l.nb_sign_changes - rk.nb_sign_changes == 0) {
+ a = rk;
+ }
+ } while (a.nb_sign_changes - b.nb_sign_changes != 1);
+ return _computeRootBetween(a, b);
+ } else {
+ return _computeRootBetween(l, r);
+ }
+ } else {
+ return {};
+ }
+}
+
+Polynomial1DRootComputer::Polynomial1DRootComputer(const Polynomial1D& p, const double a, const double b)
+ : m_sturm_sequence{p}, m_left{a}, m_right{b}
+{}
diff --git a/src/analysis/Polynomial1DRootComputer.hpp b/src/analysis/Polynomial1DRootComputer.hpp
new file mode 100644
index 000000000..0664216e5
--- /dev/null
+++ b/src/analysis/Polynomial1DRootComputer.hpp
@@ -0,0 +1,31 @@
+#ifndef POLYNOMIAL_1D_ROOT_COMPUTER_HPP
+#define POLYNOMIAL_1D_ROOT_COMPUTER_HPP
+
+#include <analysis/SturmSequence.hpp>
+
+#include <optional>
+
+class Polynomial1DRootComputer
+{
+ private:
+ struct PositionSignChanges
+ {
+ double position;
+ size_t nb_sign_changes;
+ };
+
+ const SturmSequence m_sturm_sequence;
+ const double m_left;
+ const double m_right;
+
+ double _computeRootBetween(const PositionSignChanges& l, const PositionSignChanges& r) const;
+
+ public:
+ std::optional<double> getFirstRoot() const;
+
+ Polynomial1DRootComputer(const Polynomial1D& p, const double a, const double b);
+
+ ~Polynomial1DRootComputer() = default;
+};
+
+#endif // POLYNOMIAL_1D_ROOT_COMPUTER_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 68f1bba36..4a602a294 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -84,6 +84,7 @@ add_executable (unit_tests
test_OStreamProcessor.cpp
test_ParseError.cpp
test_Polynomial1D.cpp
+ test_Polynomial1DRootComputer.cpp
test_PugsAssert.cpp
test_PugsFunctionAdapter.cpp
test_PugsUtils.cpp
diff --git a/tests/test_Polynomial1DRootComputer.cpp b/tests/test_Polynomial1DRootComputer.cpp
new file mode 100644
index 000000000..e17509f29
--- /dev/null
+++ b/tests/test_Polynomial1DRootComputer.cpp
@@ -0,0 +1,45 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/Polynomial1DRootComputer.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Polynomial1DRootComputer", "[analysis]")
+{
+ SECTION("No root in interval")
+ {
+ Polynomial1D p({-1, 0, 0, 1});
+
+ Polynomial1DRootComputer computer(p, -4, -3);
+
+ std::optional optional_root = computer.getFirstRoot();
+
+ REQUIRE(optional_root.has_value() == false);
+ }
+
+ SECTION("one root in interval")
+ {
+ Polynomial1D p({-1, 0, 0, 1});
+
+ Polynomial1DRootComputer computer(p, -0.3, 1.2);
+
+ std::optional optional_root = computer.getFirstRoot();
+
+ REQUIRE(optional_root.has_value() == true);
+ REQUIRE(optional_root.value() == Catch::Approx(1));
+ }
+
+ SECTION("multiple roots in interval")
+ {
+ Polynomial1D p =
+ Polynomial1D{{-0.5, 1}} * Polynomial1D{{-0.23, 1}} * Polynomial1D{{-0.27, 1}} * Polynomial1D{{-0.43, 1}};
+
+ Polynomial1DRootComputer computer(p, -0.3, 1.2);
+
+ std::optional optional_root = computer.getFirstRoot();
+
+ REQUIRE(optional_root.has_value() == true);
+ REQUIRE(optional_root.value() == Catch::Approx(0.23));
+ }
+}
--
GitLab