Intro

“Thank you for your reply. I really feel helpless about the problem” — I wasn’t sure whether the answer was satisfactory or unsatisfactory: Do you still feel helpless? Or you felt helpless before and now it’s okay? Here’s an explanation, in case you’re asking. Plus an improvement, for fun.
The substitutions
4
x
=πu=π
8
w
give the equivalent integrals:
∞
∫
0
|sin
4
x
|
x
+
2
x
x
=
∞
∫
0
1/8
π
|sinπu|
(4+4
3/8
π
3/8
u
)
7/8
u
u
=
∞
∫
0
2
1/8
π
|sinπ
8
w
|
1+
3/8
π
3
w
w
=
1
∫
0
2
1/8
π
|sinπ
8
w
|
1+
3/8
π
3
w
w
+
∞
∫
1
1/8
π
|sinπu|
(4+4
3/8
π
3/8
u
)
7/8
u
u
≈
1
∫
0
2
1/8
π
|sinπ
8
w
|
1+
3/8
π
3
w
w
+
n
∑
k=1
k+1
∫
k
1/8
π
|sinπu|
(4+4
3/8
π
3/8
u
)
7/8
u
u
+
∞
∫
n+1
1/8
π
2/π
(4+4
3/8
π
3/8
u
)
7/8
u
u
The integrand
1/8
π
|sinπu|
(4+4
3/8
π
3/8
u
)
7/8
u
=
k
(-1)
1/8
π
sinπu
(4+4
3/8
π
3/8
u
)
7/8
u
is analytic in a (complex) neighborhood of the interval
k≤u≤k+1
, for a positive integer
k
. For
k=0
, the
w
transformation gets rid of the singularity, and the integrand becomes analytic. This means Gaussian quadrature, including Gauss-Kronrod, will converge rapidly and more so the further from
u=0
. In fact, other than the
w
integral, only
k=1
needs a somewhat high-order Gauss-Kronrod rule. The problem, of course, is that the absolute value requires us to sum infinitely many such integrals. One can look to extrapolation to approximate the sum or another approximation. Extrapolation failed for me. So I turned to approximating the integral.
The approximation
k+1
∫
k
1/8
π
|sinπu|
(4+4
3/8
π
3/8
u
)
7/8
u
u≈
k+1
∫
k
1/8
π
2/π
(4+4
3/8
π
3/8
u
)
7/8
u
u
is quite good, with an error less than
6×
-4
10
/
13/4
k
for
k≥n
, with
n
larger than, say, 1000, based numerical estimation of the constant. You didn’t say how precise an answer you wanted. Ten digits is more than enough for most applications except for things like building math libraries, and this problem doesn’t seem such a problem. Nonetheless, one can obtain a machine precision result by computing with a higher working precision.
You didn’t say you wanted a fast way, either. The analytic, positive integrand means this integral lacks some of the worse features integration has to deal with. We can bypass NIntegrate[], predict the error, and directly use Gauss or Gauss-Kronrod. It can be done in a parallelizable way that is comparatively fast. It’s kinda sad, for me, to bypass NIntegrate[] because it’s one my favorite pieces of software. But I like speed, too, especially combined with accuracy.

Setup for the improvements

The OP’s integrand, the integrand ufn for the
u
integral, and the integrand wfn for the
w
integral. The sine function in the integrand is replaced by the variable sine. The advantage builds on the fact that the values of the sine function are identical in every integral over an interval
k≤u≤k+1
. The variable sine will be replaced by the values of the sine function at the Gauss-Kronrod nodes when computing the integral.
In[]:=
integrand=
Abs[Sin[x^4]]
x
+x^2
;​​ufn=integrand*Dt[x,u]/.x->(Pi*u)^(1/4)//Simplify[#,u>0]&​​ufn0=ufn/.f_Abs:>(sinefn=f;sine)​​sinefn​​wfn=ufn*Dt[u,w]/.u->w^8//Simplify[#,w>0]&​​avgSine=Integrate[Sin[πu],{u,0,1}]
Out[]=
1/8
π
Abs[Sin[πu]]
(4+4
3/8
π
3/8
u
)
7/8
u
Out[]=
1/8
π
sine
(4+4
3/8
π
3/8
u
)
7/8
u
Out[]=
Abs[Sin[πu]]
Out[]=
2
1/8
π
Abs[Sin[π
8
w
]]
1+
3/8
π
3
w
Out[]=
2
π

Machine precision integration

Get the nodes and weights for the
u
and
w
integrals.
In[]:=
{wvalsMP,wwtsMP,werrwtsMP}=NIntegrate`GaussKronrodRuleData[24,MachinePrecision];​​{uvalsMP,uwtsMP,uerrwtsMP}=NIntegrate`GaussKronrodRuleData[11,MachinePrecision];​​sineValsMP=sinefn/.u->uvalsMP;
The following table show the error estimate of the machine precision integral over
[k,k+1]
for the selected number of Gauss points:
In[]:=
Join[{​​{k,"integral/{k,k+1}","error"},​​{0,{wwtsMP,werrwtsMP}.(wfn/.w->wvalsMP)}//Flatten,​​{1,{wwtsMP,werrwtsMP}.(ufn/.u->1+wvalsMP)}//Flatten​​},​​Table[Prepend[{uwtsMP,uerrwtsMP}.(ufn/.u->n+uvalsMP),n],{n,1+10^Range[0,8]}]​​]//Grid
Out[]=
k
integral/{k,k+1}
error
0
0.231138
7.44847×
-17
10
1
0.0473771
-4.87229×
-18
10
2
0.0262558
-3.29549×
-18
10
11
0.00448037
1.03857×
-18
10
101
0.000332761
-2.01864×
-18
10
1001
0.0000202316
1.2723×
-18
10
10001
1.17112×
-6
10
-7.34834×
-19
10
100001
6.66455×
-8
10
-4.14327×
-19
10
1000001
3.76656×
-9
10
1.01848×
-20
10
10000001
2.12257×
-10
10
2.20494×
-20
10
100000001
1.19468×
-11
10
1.27475×
-20
10
The error for the
u
does not decrease with the value of the integral (relative error is increasing). This does not change if we increase the number of Gauss points above 9. We have hit some sort of roundoff plateau below which we cannot reduce the error estimate. That does not mean that a rule with more than 9 Gauss points does not improve the integral, only that the Kronrod error estimate is not sensitive to it. So I put it a little above 9, and we can test the result later.

32-digit precision integration

This does the same thing as the machine precision code, but with 32-digit numbers and higher order integration rules for 32-digit precision.
{wvals32,wwts32,werrwts32}=NIntegrate`GaussKronrodRuleData[42,32];​​{uvals32,uwts32,uerrwts32}=NIntegrate`GaussKronrodRuleData[13,32];​​sineVals32=sinefn/.u->uvals32;
Error of the 32-digit precision integral over
[k,k+1]
:
In[]:=
Join[{​​{k,"integral/{k,k+1}","error"},​​{0,{wwts32,werrwts32}.(wfn/.w->wvals32)}//Flatten,​​{1,{wwts32,werrwts32}.(ufn/.u->1+wvals32)}//Flatten​​},​​Table[Prepend[{uwts32,uerrwts32}.(ufn/.u->n+uvals32),n],{n,1+10^Range[0,8]}]​​]//Grid
Out[]=
k
integral/{k,k+1}
error
0
0.2311384906318081903288504956633
0.×
-32
10
1
0.04737708190777577379740933114976
0.×
-33
10
2
0.02625581791586871941246381884384
-1.5×
-31
10
11
0.004480371667189347662606269254343
-2.9×
-32
10
101
0.0003327605726340999972696332424722
-3.6×
-33
10
1001
0.00002023162407976137571785423938706
-2.23×
-34
10
10001
1.1711221356237429197693654520479×
-6
10
-1.29×
-35
10
100001
6.664548141731831615030277034775×
-8
10
-7.3×
-37
10
1000001
3.766560209189189755329709501435×
-9
10
-4.1×
-38
10
10000001
2.122573166208294997449552115348×
-10
10
-2.33×
-39
10
100000001
1.1946756983828543296287669380904×
-11
10
-1.31×
-40
10

Working precision = machine precision

The main advantage of compiling the integral over
k≤u≤k+1
is parallelization of the many terms to be computed. Dot[] is fast library call, whether compiled or not. The
u
integral over
1≤u≤2
uses the higher order rule used for the
w
integral. Summing with Method->”CompensatedSummation” saves one or two bits of roundoff error now and then. Neither a great benefit nor a great cost, overall.
termSine=With[{fn=ufn/._Abs->sine},​​Compile[(*integralk<=u<=k+1,integern*)​​{{k,_Real},{u1,_Real,1},{sine,_Real,1},{wts,_Real,1}},​​Block[{u=k+u1},​​wts.fn],​​RuntimeAttributes->{Listable},Parallelization->True]​​];​​i0=wwtsMP.(wfn/.w->wvalsMP)+(*integral0<=w<=1*)​​wwtsMP.(ufn/.u->1+wvalsMP);(*integral1<=u<=2*)​​(ivalsMP=Table[​​i0+Total[termSine[Range[ncycles,2,-1],uvalsMP,sineValsMP,uwtsMP],Method->"CompensatedSummation"]+​​N[Integrate[ufn/._Abs:>avgSine,{u,ncycles+1,Infinity}],32],​​{ncycles,10^Range[3,8]}])//AbsoluteTiming
{15.190387,​​{0.6182229382064307,​​0.6182229377606658,​​0.6182229377580635,​​0.6182229377580486,​​0.6182229377580486,​​0.6182229377580486}}
In[]:=
Differences@ivalsMP
Out[]=
{-4.45765×
-10
10
,-2.60236×
-12
10
,-1.4877×
-14
10
,0.,0.}
Note that the first to converge, ncycles=10^6, takes about 0.4 sec, most of which is Integrate[]. The million terms take 0.14 sec.

Working precision = 32

i0=wwts32.(wfn/.w->wvals32)+(*integral0<=u<=1*)​​wwts32.(ufn/.u->1+wvals32);(*integral1<=u<=2*)​​(ivals32=Table[​​i0+Total@ParallelTable[uwts32.(ufn0/.{u->n+uvals32,sine->sineVals32}),{n,2,ncycles}]+​​N[Integrate[ufn/._Abs:>avgSine,{u,ncycles+1,Infinity}],32],​​{ncycles,10^Range[3,7]}])//AbsoluteTiming
Compare high precision with machine precision. They differ in the last bit.

Numerical exploration of the approximation