supnorm fails on an rather easy case
Alexander N. Godunov reports the following problem:
> I = [-5.5* 2^(-17) ; 5.625* 2^(-17)];
> P = x-x^2/2+x^3/3-x^4/4 + 7205759403792815 * 2^(-55) * x^5 -11728124044213 * 2^(-46) * x^6 + 5146953844507055 * 2^(-55)*x^7;
> f = log(1+x);
> supnorm(P,f,I,absolute, 1b-20);
Warning: during supnorm computation, no suitable Taylor form could be found.
Warning: an error occurred during supremum norm computation. A safe enclosure of the supremum norm could not be computed.
Warning: the supremum norm on the error function between the given polynomial and the given function could not be computed.
Warning: the given expression or command could not be handled.
error
However, as demonstrate the following lines (which reproduce more or less what supnorm
is supposed to do on such an example), the computation of the supnorm should pose no problem here.
Let us start by defining two useful procedures:
-
dirtymax
takes an expression f and an interval d and returns a couple [|t,y|] where t belongs to d and |f(t)| is expected to be very close to the supnorm of f on d. The value y is a safe underestimate of |f(t)|. -
horner_interval
takes a list of interval coefficients [|c0,...,cn|] and an interval d and returns an interval R such that forall x in d and all a0 in c0, a1 in c1, ... an in cn, (sum ai*x^i) belongs to R.
> dirtymax = proc(f,d) {
var M, argM, L, t, r;
M=0;
L=dirtyfindzeros(diff(f),d);
for t in inf(d).:L:.sup(d) do {
r=evaluate(abs(f),[t]);
if inf(r)>M then {
M=inf(r);
argM = t;
};
};
return [|argM, M|];
};
> horner_interval = proc(coeffs, d) {
var n, i, res;
res = [0];
n = length(coeffs)-1;
for i from n to 0 by -1 do {
res = res*d+coeffs[i];
};
return res;
};
Then, we compute with dirtymax
a value which is surely an underestimate of ||P-f||, and probably a very good approximation of it. We use taylorform
to get a polynomial approximation f_approx
of f
, whose approximation error is safely enclosed in err
.
> R = dirtymax(P-f,I);
> xstar = R[0];
> a = R[1];
> TL = taylorform(f, 9, 0, I, absolute);
> f_approx = TL[0];
> err = horner_interval(TL[1],I)+TL[2];
If we were right about the fact that a
is a very good approximation of ||P-f||, then the following b
should be a safe upper bound of it.
> b = a*(1+2^(10-prec)) + sup(abs(err));
We now check that it is indeed a safe upper bound using numberroots
:
> numberroots(P-f_approx-b, I) == 0;
true
> sup(evaluate(P-f_approx-b, [inf(I)])) < 0;
true
> numberroots(P-f_approx+b, I) == 0;
true
> inf(evaluate(P-f_approx+b, [inf(I)])) > 0;
true
That tells us that P-f_approx-b is negative all over I and P-f_approx+b is positive all over I. In other words, |P-f_approx| < b on the whole interval.
Since xstar is a point of I for which |P-f|(xstar) is provably higher than a...
> inf(evaluate(abs(P-f), xstar)) >= a;
true
> xstar >= inf(I);
true
> xstar <= sup(I);
true
... the interval [a,b] is a safe enclosure of ||P-f||. Moreover, it is accurate to 25 bits:
> -log2((b-a)/a);
Warning: rounding has happened. The value displayed is a faithful rounding to 165 bits of the true result.
25.500342445511335229875058781202061056164278499963
> midpointmode = on!;
> [a,b];
0.10057~1/2~e-36