Skip to content
Snippets Groups Projects
Commit 56a20028 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

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.
parent ad64650f
No related branches found
No related tags found
No related merge requests found
...@@ -2,5 +2,6 @@ ...@@ -2,5 +2,6 @@
add_library( add_library(
PugsAnalysis PugsAnalysis
Polynomial1DRootComputer.cpp
SturmSequence.cpp SturmSequence.cpp
) )
#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}
{}
#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
...@@ -84,6 +84,7 @@ add_executable (unit_tests ...@@ -84,6 +84,7 @@ add_executable (unit_tests
test_OStreamProcessor.cpp test_OStreamProcessor.cpp
test_ParseError.cpp test_ParseError.cpp
test_Polynomial1D.cpp test_Polynomial1D.cpp
test_Polynomial1DRootComputer.cpp
test_PugsAssert.cpp test_PugsAssert.cpp
test_PugsFunctionAdapter.cpp test_PugsFunctionAdapter.cpp
test_PugsUtils.cpp test_PugsUtils.cpp
......
#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));
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment