Maple 15 (X86 64 LINUX) numeric integration

[Marko Riedel]

Greetings to all.

I am seeing some possible bugs in MAPLE's numeric integration routines.

Consider the following numeric integrals:

> evalf(int(cos(t)/(1+t^10), t=-infinity..infinity));
2.061888084

> evalf(int(cos(t)/(1+t^14), t=-infinity..infinity));
1.954506433

> evalf(int(cos(t)/(1+t^18), t=-infinity..infinity));
1.894657021

These values are wrong. Here is what Mathematica says:

In[1]:= NumberForm[N[Integrate[Cos[t]/(1+t^10), {t, -Infinity, Infinity}]], 15]
Out[1]//NumberForm= 1.67010286508619+(2.42861286636753*10^(-16)) I

In[2]:= NumberForm[N[Integrate[Cos[t]/(1+t^14), {t, -Infinity, Infinity}]], 15]
Out[2]//NumberForm= 1.67711957523059+(4.35983562251079*10^(-17)) I

In[3]:= NumberForm[N[Integrate[Cos[t]/(1+t^18), {t, -Infinity, Infinity}]], 15]
Out[3]//NumberForm= 1.67960327297086+(1.0899589056277*10^(-16)) I

These are correct.

Your feedback is appreciated. The explicit answer is
I(k) = 2 pi i sum_{q=0..k-1} Res(exp(it)/(1+t^{2k}); t=v(q))
where
v(q) = exp(i pi/(2k) + 2 pi i q/(2k)))

for k = 5, 6, 7.

The entire calculation can be found here:
https://groups.google.com/group/es.ciencia.matematicas/browse_thread/thread/689b080f711f2922?hl=es#

Best regards,

Marko Riedel

[Marko Riedel]

That would be k = 5, 7, 9 for I(k). The numerical results for even k by
MAPLE are correct. It's the odd k, k >= 5, that are not working.

Marko Riedel

[RGVi...@shaw.ca]

I do not have access to Maple 15, but in Maple 11 all the answers are correct (but without an error analysis). However, there is a BIG difference between using evalf(int(cos(t)/(1+t^n), t=-infinity..infinity)) and
evalf(Int(cos(t)/(1+t^n), t=-infinity..infinity)). The latter form evaluates almost instantly, while the former first computes the integral symbolically, then evaluates it numerically. Even for n = 10 it is nasty: here is the symbolic integral for n = 10:

lprint(J10);1/5*Pi^(1/2)*5^(1/2)*(1/10*Pi^(1/2)*exp(1/1000000000*10000000000^(9/10))*csc(2/5*Pi)*cot(1/5*Pi)+1/10*Pi^(1/2)*exp(-1/1000000000*10000000000^(9/10))*csc(2/5*Pi)*cot(1/5*Pi)+1/5000000000*10000000000^(9/10)*Pi^(1/2)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(2/5*Pi)*cot(3/10*Pi)*cot(1/10*Pi)*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/50000*10000000000^(2/5)*Pi^(1/2)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cos(1/5*Pi)*tan(3/10*Pi)*tan(1/10*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/100000*10000000000^(2/5)*Pi^(1/2)*exp(-1/1000000000*10000000000^(9/10))*tan(3/10*Pi)*tan(1/10*Pi)+1/50000*10000000000^(2/5)*Pi^(1/2)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(2/5*Pi)*tan(3/10*Pi)*tan(1/10*Pi)*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/50000*10000000000^(2/5)*Pi^(1/2)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(2/5*Pi)*tan(3/10*Pi)*tan(1/10*Pi)*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))+1/50000*10000000000^(2/5)*Pi^(1/2)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cos(2/5*Pi)*tan(3/10*Pi)*tan(1/10*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))+1/100000*10000000000^(2/5)*Pi^(1/2)*exp(1/1000000000*10000000000^(9/10))*tan(3/10*Pi)*tan(1/10*Pi)+1/5*Pi^(1/2)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cot(1/5*Pi)*csc(2/5*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/5000000*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(3/5)*sin(1/5*Pi)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))-1/500000000*10000000000^(4/5)*Pi^(1/2)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cos(1/5*Pi)*cot(1/5*Pi)*csc(2/5*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/5000000000*10000000000^(9/10)*Pi^(1/2)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cos(2/5*Pi)*cot(3/10*Pi)*cot(1/10*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/50000*10000000000^(2/5)*Pi^(1/2)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*sin(1/5*Pi)*tan(3/10*Pi)*tan(1/10*Pi)*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/50000*10000000000^(2/5)*Pi^(1/2)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cos(2/5*Pi)*tan(3/10*Pi)*tan(1/10*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/500*10000000000^(1/5)*Pi^(1/2)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cot(2/5*Pi)*sec(3/10*Pi)*sin(1/10*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/5*Pi^(1/2)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cot(1/5*Pi)*csc(2/5*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))+1/1000000000*10000000000^(4/5)*Pi^(1/2)*exp(1/1000000000*10000000000^(9/10))*csc(2/5*Pi)*cot(1/5*Pi)+1/500*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(1/5)*sin(1/5*Pi)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/500000000*10000000000^(4/5)*Pi^(1/2)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cos(1/5*Pi)*cot(1/5*Pi)*csc(2/5*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))+1/5000000000*10000000000^(9/10)*Pi^(1/2)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*sin(1/5*Pi)*cot(3/10*Pi)*cot(1/10*Pi)*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))-1/500*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(1/5)*sin(1/5*Pi)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))+1/500*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(1/5)*cos(1/5*Pi)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/10000000000*10000000000^(9/10)*Pi^(1/2)*exp(1/1000000000*10000000000^(9/10))*cot(3/10*Pi)*cot(1/10*Pi)+1/500*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(1/5)*cos(1/5*Pi)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/50000*10000000000^(2/5)*Pi^(1/2)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cos(1/5*Pi)*tan(3/10*Pi)*tan(1/10*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/500000000*Pi^(1/2)*csc(2/5*Pi)*cos(1/5*Pi)*10000000000^(4/5)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/1000*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(1/5)*exp(-1/1000000000*10000000000^(9/10))-1/500000000*10000000000^(4/5)*Pi^(1/2)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cot(1/5*Pi)*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/5000000*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(3/5)*cos(1/5*Pi)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))-1/5000000*10000000000^(3/5)*Pi^(1/2)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cot(2/5*Pi)*sec(3/10*Pi)*sin(1/10*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))+1/10000000000*10000000000^(9/10)*Pi^(1/2)*exp(-1/1000000000*10000000000^(9/10))*cot(3/10*Pi)*cot(1/10*Pi)-1/5000000000*10000000000^(9/10)*Pi^(1/2)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cos(1/5*Pi)*cot(3/10*Pi)*cot(1/10*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/5000000000*10000000000^(9/10)*Pi^(1/2)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*sin(1/5*Pi)*cot(3/10*Pi)*cot(1/10*Pi)*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/50000*10000000000^(2/5)*Pi^(1/2)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*sin(1/5*Pi)*tan(3/10*Pi)*tan(1/10*Pi)*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/5*Pi^(1/2)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cot(1/5*Pi)*csc(2/5*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/5000000000*10000000000^(9/10)*Pi^(1/2)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cos(2/5*Pi)*cot(3/10*Pi)*cot(1/10*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/10000000*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(3/5)*exp(-1/1000000000*10000000000^(9/10))+1/5000000000*10000000000^(9/10)*Pi^(1/2)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(2/5*Pi)*cot(3/10*Pi)*cot(1/10*Pi)*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/1000*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(1/5)*exp(1/1000000000*10000000000^(9/10))+1/1000000000*10000000000^(4/5)*Pi^(1/2)*exp(-1/1000000000*10000000000^(9/10))*csc(2/5*Pi)*cot(1/5*Pi)-1/5000000*Pi^(1/2)*sec(3/10*Pi)*sin(1/10*Pi)*10000000000^(3/5)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/10000000*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(3/5)*exp(1/1000000000*10000000000^(9/10))+1/500000000*10000000000^(4/5)*Pi^(1/2)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cot(1/5*Pi)*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/5000000*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(3/5)*sin(1/5*Pi)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))-1/500000000*Pi^(1/2)*csc(2/5*Pi)*cos(1/5*Pi)*10000000000^(4/5)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))+1/5000000*Pi^(1/2)*sec(3/10*Pi)*csc(2/5*Pi)*sin(1/10*Pi)*10000000000^(3/5)*cos(1/5*Pi)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))-1/5000000*10000000000^(3/5)*Pi^(1/2)*exp(1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*cot(2/5*Pi)*sec(3/10*Pi)*sin(1/10*Pi)*cos(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))-1/500*10000000000^(1/5)*Pi^(1/2)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cot(2/5*Pi)*sec(3/10*Pi)*sin(1/10*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))-1/500*Pi^(1/2)*sec(3/10*Pi)*sin(1/10*Pi)*10000000000^(1/5)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/5000000*Pi^(1/2)*sec(3/10*Pi)*sin(1/10*Pi)*10000000000^(3/5)*exp(-1/1000000000*cos(2/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(2/5*Pi)*10000000000^(9/10))+1/5000000000*10000000000^(9/10)*Pi^(1/2)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cos(1/5*Pi)*cot(3/10*Pi)*cot(1/10*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/500000000*10000000000^(4/5)*Pi^(1/2)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cot(2/5*Pi)*cot(1/5*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/500000000*10000000000^(4/5)*Pi^(1/2)*exp(1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cot(2/5*Pi)*cot(1/5*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/5*Pi^(1/2)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*cot(1/5*Pi)*csc(2/5*Pi)*cos(1/1000000000*sin(1/5*Pi)*10000000000^(9/10))+1/500*Pi^(1/2)*sec(3/10*Pi)*sin(1/10*Pi)*10000000000^(1/5)*exp(-1/1000000000*cos(1/5*Pi)*10000000000^(9/10))*sin(1/1000000000*sin(1/5*Pi)*10000000000^(9/10)))


It gets much, much worse for n = 14 and n = 18. Also, I suspect that roundoff errors could be a serious issue when proceeding from the symbolic answer. So, definitely use evalf(Int...)) instead of the other form.

RGV

[Marko Riedel]

[...]

>
>
> It gets much, much worse for n = 14 and n = 18. Also, I suspect that
> roundoff errors could be a serious issue when proceeding from the
> symbolic answer. So, definitely use evalf(Int...)) instead of the
> other form.
>
> RGV

Hello there,

thank you for the hint, I just tried "Int" instead of "int" and it gives
the right answers:

> evalf(Int(cos(t)/(1+t^10), t=-infinity..infinity));
1.670102865

> evalf(Int(cos(t)/(1+t^14), t=-infinity..infinity));
1.677119575

> evalf(Int(cos(t)/(1+t^18), t=-infinity..infinity));
1.679603273

> evalf(Int(cos(t)/(1+t^22), t=-infinity..infinity));
1.680769183

So the error with "int" lies not with the numerics but with the symbolic
integrator. As I pointed out above where I give a closed formula these
integrals have simple explicit representations, namely

> print(I1);
exp(t I)
k -> 2 I Pi add(residue(----------, t = v(q, k)), q = 0 .. k - 1)
(2 k)
1 + t

> print(v);
1/2 I Pi Pi q I
(q, k) -> exp(-------- + ------)
k k



> I1(5);
/
| (3/5) (1/10) (4/5) (3/10)
2 I Pi |-1/10 exp((-1) ) (-1) - 1/10 exp((-1) ) (-1)
|
\

7/10 9/10 \
(-1) (-1) |
- 1/10 I exp(-1) - 1/10 -------------- - 1/10 --------------|
(1/5) (2/5) |
exp((-1) ) exp((-1) )/


> I1(7);
/
|
|
| (4/7) (1/14) (5/7) (3/14)
2 I Pi |-1/14 exp((-1) ) (-1) - 1/14 exp((-1) ) (-1)
|
\

9/14
(6/7) (5/14) (-1)
- 1/14 exp((-1) ) (-1) - 1/14 I exp(-1) - 1/14 --------------
(1/7)
exp((-1) )

11 13 \
-- -- |
14 14 |
(-1) (-1) |
- 1/14 -------------- - 1/14 --------------|
(2/7) (3/7) |
exp((-1) ) exp((-1) )/

> I1(9);
/
|
| 1/2
| (5/9) (1/18) 2 exp(1/2 I 3 ) exp(-1/2)
2 I Pi |-1/18 exp((-1) ) (-1) + ---------------------------
| 1/2
\ 18 I - 18 3

(7/9) (5/18) (8/9) (7/18)
- 1/18 exp((-1) ) (-1) - 1/18 exp((-1) ) (-1)

11 13
-- --
18 18
(-1) (-1)
- 1/18 I exp(-1) - 1/18 -------------- - 1/18 --------------
(1/9) (2/9)
exp((-1) ) exp((-1) )

17 \
-- |
1/2 18 |
2 exp(-1/2 I 3 ) exp(-1/2) (-1) |
- ---------------------------- - 1/18 --------------|
1/2 (4/9) |
-18 I - 18 3 exp((-1) )/

I can confirm that the symbolic integrator produces truly massive
results for these integrals. I guess the question to ask here is why
does the integrator not evaluate these integrals through residues, which
gives compact answers that are orders of magnitude smaller and
mathematically correct.

With the above formulas care must be taken to use the same branch of the
logarithm for all exponentials. In the expressions above MAPLE has
turned exp(i pi) from the original formula (see my first message) back
into -1, which, when exponentiated, is multivalued. That's why the
powers of -1 should be computed with a branch of the logarithm that has
log(-1) = i pi.

Here are some truly scary statistics:

> length(convert(int(cos(t)/(1+t^10), t=-infinity..infinity),string));
8806

> length(convert(int(cos(t)/(1+t^14), t=-infinity..infinity),string));
22476

> length(convert(int(cos(t)/(1+t^18), t=-infinity..infinity),string));
42751

> length(convert(int(cos(t)/(1+t^22), t=-infinity..infinity),string));
84282



The output of the integrator shows exponential growth! (Compare to my
formula, where the number of terms is linear.)

Best regards,

Marko Riedel

[Marko Riedel]

Greetings,

the formula as implemented in "I1" definitely does not produce
exponential size output as can be seen here:

> st := [seq(length(convert(I1(2*m+1),string)), m=2..14)];
st := [156, 224, 325, 370, 442, 537, 586, 658, 749, 802, 864, 964, 1018]

> seq(st[m+1]-st[m], m=1..nops(st)-1);
68, 101, 45, 72, 95, 49, 72, 91, 53, 62, 100, 54

The fluctuations are to be expected as for some k the roots of 1+t^{2k}
have more compact representations than others.

Marko Riedel

[RGVi...@shaw.ca]

My previous response used Maple 11 on an HP Probook. Now I am running Maple 14 on an HP Pavilion, and the results are different: the evalf(int...) results are no longer correct (so something broke in the symbolic integrator between Maple 11 and Maple 14). However, the evalf(Int...) results are the same. Possibly, one can change variables in such a way that Maple 14-15 uses a different symbolic integral---perhaps using residues, as you suggest---but doing that is inconvenient.

RGV

[Axel Vogt]

Besides what was already said (difficult to read ...)

To compute Int(cos(t)*f(t), t=-infinity..infinity)), f(t)= 1/(1+t^(2*k))
with Maple one can use with(inttrans) and then either use Cosine

2/sqrt(2/Pi)*fouriercos( 1/(1+t^14),t,1); # as an example
evalf(%);

or the Fourier transform

fourier( 1/(1+t^14),t,1);
evalf(%);

[Axel Vogt]

To do it purely numerically first set Digits:=15, which roughly means to work
in double precision. Then note symmetry w.r.t. 0, so half the axis is enough.

More or less the integrals are 1.6 and for double precision one therefore can
ignore the tail, if it is less than epsilon:=evalhf(DBL_EPSILON).

Define upper:= k -> (epsilon*(2*k-1))^(1/(-2*k+1)), which ensures that, and
now it can be reduced to a finite interval:

2*Int( cos(t)/(1+t^(2*k)), t=0..upper(k))

Now this should do even for large values of k (since upper ---> 1).

One should compare with the Cosine transform already given.

[Marko Riedel]

Dear Axel,

thank you very much for the hint, if it were only the numeric values
that would be sufficient.

However I find the output of "fouriercos" to be of very poor quality.
Compared to my formula above it invokes Ci and Ssi, which are not
elementary functions. However after using "allvalues" the output is
correct (and has a linear number of terms) and is essentially the same
as mine. (To read the formulas in my post use a news client that
supports fixed-width fonts, like Google Groups or GNUS. You might enjoy
reading the proof of my closed form expression on
es.ciencia.matematicas, see the link above.)

On the other hand the output of "fourier" is a near hit after passing it
through "allvalues".

I still believe that Maple ought not to output incorrect numerics
without warning, as per my earlier message. The exponential memory usage
does not exactly help it shine either (as with
int(cos(t)/(1+t^14),t=-infinity..infinity).) On the other hand the
results of Integrate[Cos[t]/(1 + t^14), {t, -Infinity, Infinity}] are
just as bad. There are fluctuations but overall it looks exponential,
too.

In[5]:= ByteCount[Integrate[Cos[t]/(1+t^10), {t, -Infinity, Infinity}]]
Out[5]= 16512
In[6]:= ByteCount[Integrate[Cos[t]/(1+t^14), {t, -Infinity, Infinity}]]
Out[6]= 155880
In[7]:= ByteCount[Integrate[Cos[t]/(1+t^18), {t, -Infinity, Infinity}]]
Out[7]= 188248
In[8]:= ByteCount[Integrate[Cos[t]/(1+t^22), {t, -Infinity, Infinity}]]
Out[8]= 270136
In[9]:= ByteCount[Integrate[Cos[t]/(1+t^26), {t, -Infinity, Infinity}]]
Out[9]= 388488
In[10]:= ByteCount[Integrate[Cos[t]/(1+t^30), {t, -Infinity, Infinity}]]
Out[10]= 581304

I believe there is no shame in having a lookup table be part of the
integrator to catch those cases where the algorithm heads down the wrong
path. There are many good books out there that describe integration by
residues, which is a very powerful method, especially when combined with
the considered deployment of the right branches of the logarithm.

Best regards,

Marko

[Marko Riedel]

One more thing, should anyone with access to MAPLE's internal algorithms
be reading this, please pay attention to the fact that the bug only
occurs with odd k in I(k):

> evalf(I1(5));
1.670102866 + 0. I

> evalf(int(cos(t)/(1+t^10), t=-infinity..infinity));
2.061888084

> evalf(I1(6));
-10
1.674641305 - 0.6283185308 10 I

> evalf(int(cos(t)/(1+t^12), t=-infinity..infinity));
1.674641307

> evalf(I1(7));
1.677119575 + 0. I

> evalf(int(cos(t)/(1+t^14), t=-infinity..infinity));
1.954506433

> evalf(I1(8));
-10
1.678621916 - 0.6283185308 10 I

> evalf(int(cos(t)/(1+t^16), t=-infinity..infinity));
1.678621913

> evalf(I1(9));
1.679603273 + 0. I

> evalf(int(cos(t)/(1+t^18), t=-infinity..infinity));
1.894657020


So this definitely looks like a software error rather than the
accumulated result of numeric errors.

Best regards,

Marko

[Axel Vogt]

Hm, well the Cosine transform looks a bit crazy because of the
strange 3-valued signum, but it actually says: roots in the
upper plane, the rest involves roots of unity and some exp.

For Numerics I gave some fast and accurate solution.

Concerning Maple's possible bug versus Mathematica:

That is difficult to judge. First we do not *know* whether MMA
uses some symbolics inbetween (Maple does not). And second it
is for sure, that Maple uses only 10 Digits (in Your setting),
while MMA uses much more.

Anyway: I also would expect Maple to get it (may be with some
higher Digits, issuing warnings). And it does not know the
generic solution (from Residue calculus). But you have 2 ways
which work to get the numerical value.

And as already said: do *not* use evalf(lOWERCASE int(...)).

[Marko Riedel]

Hello there,

thanks for your feedback. If you go back to the solution from my first
post yesterday and read it carefully (displayed with a fixed-width font)
you'll see that it uses precisely the singularities in the upper half
plane, so I can confirm what you are saying. The contour used in the
proof is a half-circle in the upper half plane and that fact plays a key
role in some of the growth estimates, because the integral along the
half-circle in the lower half plane does not disappear in the limit
(with R -> oo).

I tried with many more digits (like 30) but the bug remains. I am
convinced that this is a coding/software error or else we would not have
the pattern with it working for even k and failing for odd k. For
example the related integral int_{-oo}^oo cos(t)^q/(1+t^2) dt has
different closed form expressions for q odd and q even. I bet this is
precisely what is happening with I(k).

Marko

[Marko Riedel]

I didn't express myself well, the closed form for I(k) is the same
whether k is even or odd, I am refering to MAPLE's algorithm, which
seems to use a different formula for k even and k odd.

Marko

[Marko Riedel]

I used evalf(lOWERCASE int(...)) after reading your message, which I
understand and will remember for the future, in order to show that the
symbolic value of the integral is incorrect.

Marko