Project

General

Profile

Bug #1171

RealRoots: first point is sometimes wrong?

Added by Anna Maria Bigatti 22 days ago. Updated 21 days 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 22 days 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 22 days 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 22 days 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 22 days ago

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

#5 Updated by John Abbott 21 days 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 21 days 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?

Also available in: Atom PDF