UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-01-19 00:37:25 +01:00
Code Issues Packages Projects Releases Wiki Activity

[poincare/expression] Fix solutions to e^x=0

Browse Source 
Because of the limitations of the floating-point representation, e^x is
null for x <= 710, causing the nextRoot functions to find roots for it.
To prevent this, we looking for places where a function changes sign, we
actually require the function to take two non-null values of different
signs.
This commit is contained in:
Gabriel Ozouf
2020-12-09 15:50:01 +01:00
committed by EmilieNumworks
parent 78cb340065
commit c7b758f536
3 changed files with 30 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
apps/solver/test/equation_store.cpp
View File

@@ -82,6 +82,8 @@ QUIZ_CASE(equation_solve) {
assert_solves_numerically_to("cos(x)=0", -100, 100, {-90.0, 90.0});
assert_solves_numerically_to("cos(x)=0", -900, 1000, {-810.0, -630.0, -450.0, -270.0, -90.0, 90.0, 270.0, 450.0, 630.0, 810.0});
assert_solves_numerically_to("√(y)=0", -900, 1000, {0}, "y");
assert_solves_numerically_to("ℯ^x=0", -1000, 1000, {});
assert_solves_numerically_to("ℯ^x/1000=0", -1000, 1000, {});
// Long variable names
assert_solves_to("2abcde+3=4", "abcde=1/2");

24
poincare/src/expression.cpp
View File

@@ -1005,6 +1005,9 @@ Coordinate2D<double> Expression::nextMaximum(const char * symbol, double start,
}
double Expression::nextRoot(const char * symbol, double start, double step, double max, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
if (nullStatus(context) == ExpressionNode::NullStatus::Null) {
return start + step;
}
return nextIntersectionWithExpression(symbol, start, step, max,
[](double x, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const void * context1, const void * context2, const void * context3) {
const Expression * expression0 = reinterpret_cast<const Expression *>(context1);
@@ -1139,17 +1142,32 @@ double Expression::nextIntersectionWithExpression(const char * symbol, double st
void Expression::bracketRoot(const char * symbol, double start, double step, double max, double result[2], Solver::ValueAtAbscissa evaluation, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, const Expression expression) const {
double a = start;
double c = a;
double b = start+step;
double fa = evaluation(a, context, complexFormat, angleUnit, this, symbol, &expression);
double fc = fa;
double fb = evaluation(b, context, complexFormat, angleUnit, this, symbol, &expression);
while (step > 0.0 ? b <= max : b >= max) {
double fa = evaluation(a, context, complexFormat, angleUnit, this, symbol, &expression);
double fb = evaluation(b, context, complexFormat, angleUnit, this, symbol, &expression);
if (fa*fb <= 0) {
if (fa == 0. && ((fc < 0. && fb > 0.) || (fc > 0. && fb < 0.))) {
/* If fa is null, we still check that the function changes sign on ]c,b[,
* and that fc and fb are not null. Otherwise, it's more likely those
* zeroes are caused by approximation errors. */
result[0] = a;
result[1] = b;
return;
} else if (fb != 0. && ((fa < 0.) != (fb < 0.))) {
/* The function changes sign.
* The case fb = 0 is handled in the next pass with fa = 0. */
result[0] = a;
result[1] = b;
return;
}
c = a;
fc = fa;
a = b;
fa = fb;
b = b+step;
fb = evaluation(b, context, complexFormat, angleUnit, this, symbol, &expression);
}
result[0] = NAN;
result[1] = NAN;

8
poincare/test/function_solver.cpp
View File

@@ -236,9 +236,15 @@ QUIZ_CASE(poincare_function_root) {
{
constexpr int numberOfRoots = 1;
Coordinate2D<double> roots[numberOfRoots] = {
Coordinate2D<double>(99.8, 0.0)};
Coordinate2D<double>(99.9, 0.0)};
assert_points_of_interest_are(PointOfInterestType::Root, numberOfRoots, roots, "0", nullptr, "a", 100.0, -0.1, -1.0);
}
{
constexpr int numberOfRoots = 1;
Coordinate2D<double> roots[numberOfRoots] = {
Coordinate2D<double>(NAN, 0.0)};
assert_points_of_interest_are(PointOfInterestType::Root, numberOfRoots, roots, "ℯ^x", nullptr, "a", -1000.0, 0.1, -800);
}
}
QUIZ_CASE(poincare_function_intersection) {