# Numerically evaluating an integral related to Cantor's staircase

Cantor's staircase $F_C(x)$ is a well-known "pathological" function:

Plot[CantorStaircase[x], {x, 0, 1}] The MathWorld link given above claims that

$$\int_0^1 F_C(x)^{F_C(x)}\mathrm dx\approx 0.750387\dots$$

NIntegrate[], however, is not sufficiently able to deal with this integral without any help.

One recourse (apparently done in the MathWorld notebook) would be to split the unit interval into the Cantor set, integrate over those subintervals (where the staircase ought to be constant), and add up those results. This is not very efficient, since the number of subintervals grows exponentially.

Is there a more efficient/elegant way to numerically compute $\int_0^1 F_C(x)^{F_C(x)}\mathrm dx$ (say, to fifty digits)?

My own attempts yield a result $\approx 0.75037204$, which is a bit different from the posted result. I am not sure if this is only because I undersampled the function with the Cantor set splitting I described above.

• NIntegrate[CantorStaircase[x]^CantorStaircase[x], {x, 0, 1}, Method -> "DoubleExponential", MinRecursion -> 15, MaxRecursion -> 50] agrees with your solution, producing 0.750372. – shrx Dec 10 '15 at 17:02
• How long did that evaluation take, @shrx? – J. M.'s discontentment Dec 10 '15 at 17:21
• 69.8 seconds. You can drop MinRecursion to 10 and still get the same result (at least up to the digits displayed in the output) in 4.6 seconds. – shrx Dec 10 '15 at 17:26
• You're looking for NLebesgue[ ].It's coming in 10.5. – Dr. belisarius Dec 11 '15 at 2:22
• @shrx On my laptop I get the same result under 1 second with NIntegrate[CantorStaircase[x]^CantorStaircase[x], {x, 0, 1}, Method -> {"GlobalAdaptive", "SingularityHandler" -> None}, MinRecursion -> 15, MaxRecursion -> 50] – Anton Antonov Feb 21 '16 at 16:52

I don't have 50 digits. I can get close to machine precision reasonably quickly though. More importantly, there are good error bounds associated with the procedure.

Using the fact that $F_C$ assumes the values $(2k-1)/2^n$ on the intervals of length $3^n$ removed at the $n^{\text{th}}$ stage of construction of the Cantor set, we can write the integral as a double sum:

\begin{align} \int_0^1 F_C(x)^{F_C(x)} dx &= \sum_{n=1}^{\infty}\frac{1}{3^n} \sum_{k=1}^{2^{n-1}} \left(\frac{2k-1}{2^n}\right)^{(2k-1)/2^n} \\ &= \sum_{n=1}^{\infty}\frac{2^{n-1}}{3^n} \sum_{k=1}^{2^{n-1}} \left(\frac{2k-1}{2^n}\right)^{(2k-1)/2^n} \frac{1}{2^{n-1}}. \end{align}

After the little rewrite to get to the last step, the inner sum is a midpoint based Riemann sum for the integral $$I = \int_0^1 x^x dx.$$ Thus, given an $M\in\mathbb N$ we have the estimate $$\sum_{n=1}^{M}\frac{1}{3^n} \sum_{k=1}^{2^{n-1}} \left(\frac{2k-1}{2^n}\right)^{(2k-1)/2^n} + \frac{I}{2} \sum_{n=M+1}^{\infty} \left(\frac{2}{3}\right)^n.$$ Of course, the final geometric series can easily be evaluated in closed form. Using this, we can generate the following table of estimates in about 5 seconds:

i = NIntegrate[x^x, {x, 0, 1}];
results = Table[
N[Sum[Sum[((2 k - 1)/2^n)^((2 k - 1)/2^n), {k, 1, 2^(n - 1)}]/
3^n, {n, 1, m}]] + (2/3)^m*i, {m, 0, 16}];
Column[FullForm /@ results] To 12 digits, I guess the answer is $0.750372036605$.

Here is another way to look at the error estimate. First, the sum used to approximate the integral $I$ is a midpoint sum. Thus, associated error is of the order $4^{-m}$. Furthermore, that integral is multiplied by $(2/3)^m$ yielding a total absolute error of the order $1/6^m$. Note that $1/6^{16} \approx 3\times10^{-13}$, which makes the 12 digits of accuracy seem quite reasonable.

We should even be able to use Richardson extrapolation to squeeze a couple more digits out of the process:

Column[
Drop[FullForm[(6/5) #[] - (1/5) #[]] & /@
Partition[results, 2, 1], 12]] I guess a similar approach should work for $$\int_0^1 f(F_C(x)) dx,$$ whenever $f$ is a well-behaved, real function; in the case considered here, we have $f(x)=x^x$.

Finally, it might be worth considering how much improvement we get for this analysis versus, say, just truncating the sum:

$$\int_0^1 F_C(x)^{F_C(x)} dx \approx \sum_{n=1}^{M}\frac{1}{3^n} \sum_{k=1}^{2^{n-1}} \left(\frac{2k-1}{2^n}\right)^{(2k-1)/2^n}$$

or, in Mathematica, with $M=16$ as above:

N[
Sum[Sum[((2 k - 1)/2^n)^((2 k - 1)/2^n), {k, 1, 2^(n - 1)}]/3^n,
{n, 1, 16}]]
(* Out: 0.749179 *)


This doesn't even yield 2 correct digits; the second digit is almost there, since $0.783431*(2/3)^{16} \approx 0.0011927$, where $$\int_0^1 x^x dx \approx 0.78431.$$

• This looks neat! If you want to try out full Richardsonian extrapolation, here's a snippet: pts = Transpose[{6^Range[0, 1 - n, -1], results}]; Table[InterpolatingPolynomial[Take[pts, k], 0], {k, n}]. I get 0.7503720366052 with this. – J. M.'s discontentment Dec 10 '15 at 23:31
• @J.M. Glad you like it! – Mark McClure Dec 11 '15 at 1:15
• Nice answer, I tried "truncating the sum", it is nice to see how this can be improved. – Jacob Akkerboom Dec 11 '15 at 12:28