Problem 1:

What is limit(int(x^(-1)*cos(x^(-1)*log(x)),x = epsilon .. 1),epsilon = 0)  ?

Solution:

Use a change of variables y = -x^(-1)*log(x)  to make the cosine linear in the variable of integration.  

>    restart;

>    ti:= time(): Digits:= 20:

>    J:= Int(x^(-1)*cos(x^(-1)*ln(x)),x=epsilon..1);

J := Int(1/x*cos(1/x*ln(x)),x = epsilon .. 1)

>    alias(W=LambertW): xs:= solve(y=-x^(-1)*ln(x),x);

xs := 1/y*W(y)

>    simplify(PDEtools[dchange](x=xs,J));

Int(1/y*cos(y)*W(y)/(1+W(y)),y = 0 .. -1/epsilon*ln(epsilon))

Of course, limit(ln(1/epsilon)/epsilon,epsilon = 0,right) = infinity .

>    integrand:= cos(y)*LambertW(y)/(1+LambertW(y))/y;

integrand := 1/y*cos(y)*W(y)/(1+W(y))

>    answer:= Int(integrand,y=0..infinity);

answer := Int(1/y*cos(y)*W(y)/(1+W(y)),y = 0 .. infinity)

This is still not a very pleasant oscillatory integral, although the singularity at 0 is removable.  We break it up into 3 pieces.  For y from 0 to 199*Pi , use ordinary numerical integration, (but on the original form of the integral, to avoid the apparent singularity at y=0 as well as not having to evaluate W numerically; and using the _Dexp method, as otherwise Maple will have trouble because it's close to the singularity at x = 0 ).  

>    Int(1/x * cos(ln(x)/x) , x=subs( y=199*Pi, xs ) ... 1 );
r1:= evalf(Int(op(%),_Dexp));

Int(1/x*cos(1/x*ln(x)),x = 1/199/Pi*W(199*Pi) .. 1)

r1 := .32336949155993421570

For each interval (2*n-1)*Pi < y   `` < (2*n+1)*Pi  for 100 <= n ,  write the integral as   int(cos(t)*W(t+2*n*Pi)/(t+2*n*Pi)/(1+W(t+2*n*Pi)),t = -Pi .. Pi) .  Expand this in an asymptotic series in n and integrate term-by-term.   The terms for even n vanish by symmetry.

>    map(factor,asympt(cos(t)*W(t+2*n*Pi)/(t+2*n*Pi)/(1+W(t+2*n*Pi)), n, 8));

1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...
1/2*cos(t)*W(2*n*Pi)/Pi/(W(2*n*Pi)+1)/n-1/4*cos(t)*t*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/16*cos(t)*t^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/Pi^3/(W(2*n*Pi)+1)^5/n^3-1/96*cos(t)*t^...

>    J2:= map(factor,map(int, eval(%,O=0), t=-Pi..Pi));

J2 := -1/4/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3-1/96*W(2*n*Pi)^5*(24*W(2*n*Pi)^4+192*W(2*n*Pi)^3+622*W(2*n*Pi)^2+974*W(2*n*Pi)+625)*(Pi^2-6)/Pi^4/(W(2*n*Pi)+1)^9/n^5-1/768...
J2 := -1/4/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3-1/96*W(2*n*Pi)^5*(24*W(2*n*Pi)^4+192*W(2*n*Pi)^3+622*W(2*n*Pi)^2+974*W(2*n*Pi)+625)*(Pi^2-6)/Pi^4/(W(2*n*Pi)+1)^9/n^5-1/768...
J2 := -1/4/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3-1/96*W(2*n*Pi)^5*(24*W(2*n*Pi)^4+192*W(2*n*Pi)^3+622*W(2*n*Pi)^2+974*W(2*n*Pi)+625)*(Pi^2-6)/Pi^4/(W(2*n*Pi)+1)^9/n^5-1/768...
J2 := -1/4/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3-1/96*W(2*n*Pi)^5*(24*W(2*n*Pi)^4+192*W(2*n*Pi)^3+622*W(2*n*Pi)^2+974*W(2*n*Pi)+625)*(Pi^2-6)/Pi^4/(W(2*n*Pi)+1)^9/n^5-1/768...
J2 := -1/4/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3-1/96*W(2*n*Pi)^5*(24*W(2*n*Pi)^4+192*W(2*n*Pi)^3+622*W(2*n*Pi)^2+974*W(2*n*Pi)+625)*(Pi^2-6)/Pi^4/(W(2*n*Pi)+1)^9/n^5-1/768...
J2 := -1/4/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3-1/96*W(2*n*Pi)^5*(24*W(2*n*Pi)^4+192*W(2*n*Pi)^3+622*W(2*n*Pi)^2+974*W(2*n*Pi)+625)*(Pi^2-6)/Pi^4/(W(2*n*Pi)+1)^9/n^5-1/768...

>    L:= convert(%,list): evalf(eval(L,n=100));

[-.40028302126504243120e-7, -.76396717344542427096e-12, -.14718534949661344200e-16]
[-.40028302126504243120e-7, -.76396717344542427096e-12, -.14718534949661344200e-16]

Apparently by n = 100  we only need the n^(-3)  and n^(-5)  terms.  We'll use these for n from 100 to 1000.

>    J3:= op(1,J2)+op(2,J2);
r2:= add(evalf(J3),n=100..1000);

J3 := -1/4/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3-1/96*W(2*n*Pi)^5*(24*W(2*n*Pi)^4+192*W(2*n*Pi)^3+622*W(2*n*Pi)^2+974*W(2*n*Pi)+625)*(Pi^2-6)/Pi^4/(W(2*n*Pi)+1)^9/n^5
J3 := -1/4/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3-1/96*W(2*n*Pi)^5*(24*W(2*n*Pi)^4+192*W(2*n*Pi)^3+622*W(2*n*Pi)^2+974*W(2*n*Pi)+625)*(Pi^2-6)/Pi^4/(W(2*n*Pi)+1)^9/n^5
J3 := -1/4/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3-1/96*W(2*n*Pi)^5*(24*W(2*n*Pi)^4+192*W(2*n*Pi)^3+622*W(2*n*Pi)^2+974*W(2*n*Pi)+625)*(Pi^2-6)/Pi^4/(W(2*n*Pi)+1)^9/n^5

r2 := -.20381712617415828093e-5

For n > 1000, use the asymptotics of the Euler-Maclaurin form for the summation of the n^(-3)  and n^(-5)  terms.

>    J4:= map(factor,asympt(eulermac(J3,n),n));

J4 := 1/4*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/8/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3+1/32*W(2*n*Pi)^4*(6*W(2*n*Pi)^3+36*W(2*n*Pi)^2+79*W(2*n*Pi)+64)*(Pi^2...
J4 := 1/4*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/8/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3+1/32*W(2*n*Pi)^4*(6*W(2*n*Pi)^3+36*W(2*n*Pi)^2+79*W(2*n*Pi)+64)*(Pi^2...
J4 := 1/4*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/8/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3+1/32*W(2*n*Pi)^4*(6*W(2*n*Pi)^3+36*W(2*n*Pi)^2+79*W(2*n*Pi)+64)*(Pi^2...
J4 := 1/4*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/8/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3+1/32*W(2*n*Pi)^4*(6*W(2*n*Pi)^3+36*W(2*n*Pi)^2+79*W(2*n*Pi)+64)*(Pi^2...
J4 := 1/4*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/8/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3+1/32*W(2*n*Pi)^4*(6*W(2*n*Pi)^3+36*W(2*n*Pi)^2+79*W(2*n*Pi)+64)*(Pi^2...
J4 := 1/4*W(2*n*Pi)^2*(W(2*n*Pi)+2)/Pi^2/(W(2*n*Pi)+1)^3/n^2+1/8/Pi^2*W(2*n*Pi)^3*(2*W(2*n*Pi)^2+8*W(2*n*Pi)+9)/(W(2*n*Pi)+1)^5/n^3+1/32*W(2*n*Pi)^4*(6*W(2*n*Pi)^3+36*W(2*n*Pi)^2+79*W(2*n*Pi)+64)*(Pi^2...

>    r3:= -evalf(eval(J4,{O=0, n=1001}));

r3 := -.21710893457640155616e-7

>    r:= r1+r2+r3;

r := .32336743167777901648

Actual answer, as given by Trefethen to 40 digits:

             0.3233674316 7777876139 9370087952 1704466510...

>    (time()-ti)*seconds;

48.979*seconds

Back to Solutions to the SIAM $100 Challenge