## Bug #1171

### RealRoots: first point is sometimes wrong?

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 Bigatti4 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 Bigatti4 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 Abbott4 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 Abbott4 months ago

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

#### #5 Updated by John Abbott4 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 Abbott4 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 Abbott2 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?)