How do I plot the derivative of the Cantor Staircase?

I suppose I would consider — in fact, I am — what is to be understood by the numerical approximation of the derivative at an approximate real number. First of all, since the derivative CantorStaircase' is not defined symbolically in Mathematica, numerical routines like Plot[] lead to a numerical approximation of it. The numerical approximation is built into Mathematica. Second, one way of looking at a floating-point number ${\tt x}$ is that it is the rounded approximation of a real number $x$, rounded to the nearest number in the floating-point system. As such, ${\tt x}$ may represent any number in the interval between ${\tt x}-\delta_1$ and ${\tt x}+\delta_2$, which numbers are halfway between ${\tt x}$ and the next floating-point number below and above ${\tt x}$, respectively.

When a function $f$ is analytic, that a floating-point number might represent an interval is less of a concern in practice, because it is common that $f'$ varies tamely enough that $f'$ is generally close to its mean rate of change in any small interval containing $x$. However, the discrete nature of floating-point numbers mean that it is often impossible to use a finite difference formula to accurately represent the value of the derivative without using a fairly large interval. For the central difference formula, say, one whose length is on the order of $|x| \sqrt{\varepsilon}$, where $\varepsilon$ is the machine epsilon of the floating-point system, tends to be a good choice. This length we can call the "scale" of the finite difference approximation. (See ND[] or @MariuszIwaniuk's answer. The built-in procedure uses a fairly high-order formula with a fairly large scale, close to 0.5, which does poorly on functions that are not smooth to the corresponding order.)

When the function is like CantorStaircase[], which has points in its domain (1) at which it is not differentiable and (2) that are closer to each other than any given distance, then there is not a safe scale on which to base a finite-difference approximation of the derivative — at least, not one that is safe throughout the domain. Here are some points of view from which one might calculate CantorStaircase'[x]:

• For a floating-point x in a standard floating-point system, the actual number used is a rational number with a denominator which is a power of 2. So defining CantorStaircase'[_Real] = 0 will prevent numerical approximation of the derivative and yield the desired plot.

• Thinking of a floating-point x as representing an interval of real numbers, the average rate of change of a function over this interval can be taken as an approximation for the value of the derivative (compare with the Mean Value Theorem). When the interval contains points at which CantorStaircase is not differentiable, then the average rate of change can be taken as an illustration of the difficulty in differentiating it. (These mean values tend to be rather large.)

• Since Mathematica has exact representations of real numbers other than binary fractions, one might extend the symbolic definition of CantorStaircase' in the first bullet point to other numbers and numeric expressions. I don't think this can be done completely, but I don't know enough about the Cantor function to be sure.

Method 1

Clear[Derivative];
CantorStaircase'[_Real] = 0;
Plot[CantorStaircase'[x], {x, 0, 1}]

Method 2

Show the mean rate of change over the interval $[x-\varepsilon,x+\varepsilon]$. If you want a smaller interval, you need to raise the precision (via SetPrecision[]). The Cantor function does not have the same truncation/round-off error behavior as analytic functions. Also we're not really approximating the derivative. The only "values" of the derivative are zero and undefined. If the interval contains a point where the Cantor function is not differentiable, then we get a large average rate of change. Since we plot rather few points, we get rather few spikes. The average rate appears to be zero over some intervals when in fact it is not.

meanCSPrime // ClearAll;
meanCSPrime[x_Real] := If[0 <= x <= 1,
   With[{delta = SetPrecision[10^-Accuracy[x], N@Precision[x]], 
         xx = SetPrecision[x, N@Precision[x]]},
    N[(CantorStaircase[xx + delta] - 
        CantorStaircase[xx - delta])/(2 delta), Precision[x]]
    ],
   0];
(* Alt plot to specifically hit some points in the Cantor set *) 
ListPlot[
 Table[meanCSPrime[x], {x, Subdivide[0.``16, 1.`16, 2^8*3^3]}], 
 Filling -> Axis, GridLines -> {Range[3^2]/3^2, None}, 
 DataRange -> {0, 1}, PlotStyle -> Darker[Magenta, 0.2]]

Method 3

Extend the symbolic definition to all rational numbers. An alternative shown here is to use a ternary approximation of a numeric expression x to determine whether CantorStaircase'[x] is defined (and therefore 0). The definitions also leave undefined CantorStaircase'[x] for a Real x that might belong to the Cantor set. (Truncation error can lead to a mistake here, of course, but one usually puts up with this in numerical approximation. However, the nature of the Cantor function leads people to want to calculate values at these "rare" inputs, so I accept that the definition below will be disappointing to some. Feel free to post an improvement.)

(* add defs to un-Protect[]-ed Derivative; cf. ??Derivative *)
Clear[Derivative];
CantorStaircase'[x_?NumericQ] /; ! (0 <= x <= 1) := 0;
CantorStaircase'[x_Rational] /; IntegerQ[Log[3, Denominator[x]]] := 
  Indeterminate;
CantorStaircase'[x_Rational] := 0;
(* approximate: valid for the truncated ternary expansion *)
CantorStaircase'[x_Real] /; ! FreeQ[RealDigits[x, 3][[1]], 1] := 0;
$cantorMaxPrecision = 100; (* ~47 digits *)
CantorStaircase'[x_?NumericQ] /; ! 
    FreeQ[RealDigits[x, 3, $cantorMaxPrecision][[1]], 1] := 0;
(* further extensions? *)

Plot[CantorStaircase'[x], {x, 0, 1}]
(* same as first: constant zero value shown *)

Workaround:

$Version
(*"14.1.0 for Microsoft Windows (64-bit) (July 16, 2024)"*)

f[x_] := CantorStaircase[x];
Needs["NumericalCalculus`"];
g[x_] := ND[f[t], t, x, Scale -> 10^-20, Terms -> 2, WorkingPrecision -> 20];
Plot[g[x], {x, 0, 1}]

M = 100000;
L = Rationalize[Table[g[x], {x, 0, 1, 1/M}], 0];
L1 = Count[L, 0](*Counting "zeros"*)

(*99986*)

M - L1(*Counting "artifacts"*)

(*14*)

It would appear that this simple staircase has "jaggies" where no slope can be defined. I think your function lacks defined slopes, too.

Plot[CantorStaircase[x], {x, 0, 1}]

This is not an answer, but something too long to put into comments that I thought could be useful to others with more knowledge in similar subjects. I used the Wikipedia definition of the Cantor staircase to create my own version of the function where we see explicitly what is happening:

cantorFun[x0_] := Module[{x = x0},
  
  inBase3 = 
   RealDigits[x, 3, 
    100];(*you have to hope the first 1 shows up in the first 100 \
digits*)
  
  base3Digs = inBase3[[1]];
  first1 = First@FirstPosition[base3Digs, 1];
  If[IntegerQ@first1,
   base3Digs[[first1 + 1 ;; All]] = 0;
   ];
  base3Digs = base3Digs /. (2 -> 1);
  inBase3[[1]] = base3Digs;
  FromDigits[inBase3, 2]
  ]

(Tiny note: there is probably a way to optimize this to only take the base-3 expansion until there is a 1 instead of taking the first 100 digits for every $x$, but I don't know how to tell RealDigits to do that)

What we could do from here to get the derivative is find the locations where the function changes, and put (scaled) Dirac delta functions at those jumps. I am not sure how exactly to determine where those jumps occur, but my guess is it has to do with the base-3 expansion of the number.

Add-on after reading bill s's comment on the answer above, I believe looking for the "jump" locations may be infeasible, as the function is sort of fractal in the sense that the closer you zoom in, the more jumps you see and there are probably infinite jump locations from 0 to 1.

One can plot an approximation to that derivative by

Plot[(CantorStaircase[x + 1/20000] - CantorStaircase[x])*20000,{x, 0, 1}, 
WorkingPrecision -> 10000, PlotPoints -> 1000, PlotRange -> All]

Blowing up

Plot[(CantorStaircase[x + 1/20000] - CantorStaircase[x])*20000, {x, 
0.8, 0.81}, WorkingPrecision -> 10000, PlotPoints -> 1000,  PlotRange -> {0, 10^-10}, PlotTheme -> "ThickLines"]