Project

General

Profile

Bug #1171

RealRoots: first point is sometimes wrong?

Added by Anna Maria Bigatti 4 months ago. Updated 2 months ago.

Status:
In Progress
Priority:
Normal
Assignee:
-
Category:
bug
Target version:
Start date:
03 Apr 2018
Due date:
% Done:

50%

Estimated time:
Spent time:

Description

It seems that sometimes (very rarely) the first root given by RealRoots is wrong

/**/ use QQ[x];
/**/ f :=   x^25 +1500*x^24 -3472*x^23 +75*x^22 -2888*x^21 +12984*x^20 +(-22927/4)*x^19 +11578*x^18 -22029*x^17 +(18981/2)*x^16 -21819*x^15 +29021*x^14 -12855*x^13 +25232*x^12 -25419*x^11 +4494*x^10 -20702*x^9 +14894*x^8 -1920*x^7 +12502*x^6 -92*x^5 +1700*x^4 -6628*x^3 -616*x^2 -520*x +920;

/**/ RR := RealRoots(f, 1/1000);
/**/ DecimalStr(subst(f, x, RR[1].inf));
-6710019372583766710138791136389510647991731589224733661634966329959500813.670
/**/ DecimalStr(subst(f, x, RR[2].inf));
0.034
/**/ DecimalStr(subst(f, x, RR[3].inf));
-0.418


Related issues

Related to CoCoALib - Feature #1173: Upper bound for value of poly in an intervalNew2018-04-04

History

#1 Updated by Anna Maria Bigatti 4 months ago

I found this bug while debugging ApproxSolve.
I had noticed that, when the answer was wrong, it was just on the first point.

#2 Updated by Anna Maria Bigatti 4 months ago

It seems that actually the problem is in the instability of the roots: the approximate root is probably close to the real solution, but the evaluation is not close to 0.

/**/ f :=   x^25 +1500*x^24 -3472*x^23 +75*x^22 -2888*x^21 +12984*x^20 +(-22927/4)*x^19 +11578*x^18 -22029*x^17 +(18981/2)*x^16 -21819*x^15 +29021*x^14 -12855*x^13 +25232*x^12 -25419*x^11 +4494*x^10 -20702*x^9 +14894*x^8 -1920*x^7 +12502*x^6 -92*x^5 +1700*x^4 -6628*x^3 -616*x^2 -520*x +920;

/**/ [DecimalStr(subst(f, x, rr.inf)) | rr in RealRoots(f, 10^(-50))];
["-1.095",  "0.000",  "-0.000",  "0.000",  "-0.000"]

/**/ [DecimalStr(subst(f, x, rr.inf)) | rr in RealRoots(f, 10^(-100))];
["-0.000",  "0.000",  "-0.000",  "0.000",  "-0.000"]

Another example on this pathology.

/**/ f :=   x^25  -3472*x^23 +11578*x^18 +x^7 +92*x^5 +x^4 -x^3 -x^2 -x +920;
/**/ [DecimalStr(subst(f, x, rr.inf)) | rr in RealRoots(f, 10^(-100))];
--> 0 0 0

/**/ [DecimalStr(subst(f, x, rr.inf)) | rr in RealRoots(f, 10^(-10))];
--> ["-1749714807099581247924392.282",  "0.000",  "-1882915635399905369900820.013"]

#3 Updated by John Abbott 4 months ago

  • Status changed from New to In Progress
  • % Done changed from 0 to 50

One can estimate the "stability" of the root by computing the value of the derivative at that root. By Taylor expansion we have f(alpha+eps) = f(alpha) + eps*deriv(f)(alpha) + other terms.

Probably we would want an upper bound for abs(Df(beta)) where beta is a "good approximation" to alpha.

#4 Updated by John Abbott 4 months ago

  • Related to Feature #1173: Upper bound for value of poly in an interval added

#5 Updated by John Abbott 4 months ago

It is not as simple as I first thought. While the first order approximation I mentioned in comment 3 is correct, it is only a first order approximation.

I made an implementation, and I tried the example in comment 2. That worked OK.
But with the polynomial f^2+f problems arose, presumably because the first order approximation was not good.

It would be nice to have a fn which returns intervals with the guarantee that max(abs(f(x)) | x in [a,b]) < epsilon where epsilon is the bound specified by a user.

#6 Updated by John Abbott 4 months ago

I suggest a function with the following semantics:
  • call RealRootsSmallValue(f, eps)
  • result is a list of closed intervals whose interiors are disjoint such that
    - width of interval is less than eps (less than or equal?)
    - value of f over each interval is less than eps (less than or equal?)

It is permitted to return intervals which are "too small".

What do you think? Suggestions for the fn name?

#7 Updated by John Abbott 2 months ago

I am not so happy with the suggested interface in comment 6: the parameter eps is doing double duty (a bound for the accuracy of the root, and a bound for the value of the function). Having two separate parameters for these bounds would be semantically neater (but possibly more awkward for the caller?)

Comments? Suggestions?

Also available in: Atom PDF