## How do I plot Thomae's function in Mathematica?

15

5

I wanted to plot this function

$$f(x) =\begin{cases} 1 & \text{if } x= 0 \\ \tfrac1{q} & \text{if } x = \tfrac{p}{q}\\ 0 & \text{if } x \in \mathbb{R}-\mathbb{Q} \end{cases}$$

so I wrote

FuncThomae[x_] := If[ExactNumberQ[Rationalize[x]], If[x == 0, a = 1,
L = #^-1 & /@ Divisors[Numerator[Rationalize[x]]]], a = 0]


and

ListDomain[xmin_, xmax_] := Table[Outer[List, {x}, FuncThomae[x]],
{x,xmin,xmax,0.001}] // Flatten[#, 1] &


My result doesn't take all the real numbers (nor negatives) in its domain, and for $-1$ to $1$, it should have looked like so:

but my function does not cater to negatives, nor does it look like the above plot. It looks like this from 2 to 5:

Nearly similar, but not quite. Can someone help me to perfect the function?

Wanted to duplicate this

Hi there, is the code you supplied complete, I didn't seem to be able to get it to run successfully ? – image_doctor – 2012-05-13T11:38:00.187

@image_doctor it does run. what was the problem? BTW you need to plot the generated list using the ListDomain function i coded as argument for the ListPlot Function..... – The-Ever-Kid – 2012-05-13T12:42:00.160

Outer::ipnfm: "Positive machine-sized integer or Infinity expected at position 3 in Outer[List,{2.},0]." Is the error I get. a and L seem not to be defined in the code segment you have posted. But as you have an answer, perhaps it isn't important :) – image_doctor – 2012-05-13T17:31:53.813

A tiny reminder: in floating-point arithmetic (which is what Mathematica internally uses when plotting), all numbers are rational. – J. M.'s ennui – 2012-05-14T11:15:57.693

18

I'd suggest producing a list of rational numbers and then plot the function there, like so:

maxq = 100;
fracs = Table[p/q, {q, 2, maxq}, {p, 2, q}] // Flatten // DeleteDuplicates;
pq = {#, 1/Denominator @ #} & /@ fracs;

ListPlot[pq, PlotRange -> {0, 1}]


1If you define pq as pq = {#, 1/Denominator[#]}& /@ fracs you could ListPlot pq directly without having to transform it first. – Heike – 2012-05-13T12:39:57.253

@acl thanks for the lovely answer but the plot is missing the first bit it isnt plotting 1 when x is 0 but anyway thanks.... – The-Ever-Kid – 2012-05-13T13:19:40.340

btw could you do this for -1 to 1 – The-Ever-Kid – 2012-05-13T13:20:05.690

@The-Ever-Kid yes that bit is straightforward, so I didn't bother implementing it – acl – 2012-05-13T13:20:50.523

yeah cool thanks... – The-Ever-Kid – 2012-05-13T13:23:26.137

thanks to @Heike, David and The-Ever-Kid for improving this mess of an answer! – acl – 2012-05-13T17:48:24.633

@acl When maxq=10, the point ${ \frac{1}{10}, \frac{1}{10} }$ is missing from pq. Perhaps we should use Table[p/q, {q, 2, maxq}, {p, 1, q}] instead? – Michael Wijaya – 2012-05-13T18:02:13.290

@Michael Yes, I suppose you're right (I admit that I lifted this fragment from some code I wrote for a different purpose, which is why it is such a mess...) – acl – 2012-05-13T18:24:14.073

13

Why not use DiscretePlot directly?

DiscretePlot[1/Denominator[Rationalize[x]], {x, -1, 1, 1/(2*3*4*5*6*7)},
PlotRange -> {0, 1}, Filling -> None, RegionFunction -> Function[{x, y}, Abs[x] != 1]]


The RegionFunction throws out the cases where $\frac{p}{q}=1$.

12

Another possibility that avoids the generation of fractions not in lowest terms (and thus the use of Union[] or DeleteDuplicates[]) rests on generating a Farey sequence, and then applying the Dirichlet-Thomae function to that:

farey[n_Integer?Positive] := Module[{v = 0, w = 1, p = 1, q = n, t},
Join[{0, 1/n}, Flatten[Last[
Reap[While[p < q,
t = Quotient[n + w, q];
{{p, v}, {q, w}} = {{t p - v, p}, {t q - w, q}};
Sow[p/q]]]
]]]]

ListPlot[{#, 1/Denominator[#]} & /@ farey[100], Frame -> True, PlotRange -> {0, 1}]


(Nowadays, one would use the built-in FareySequence[].) – J. M.'s ennui – 2015-10-19T18:09:51.893

2

A quick look at a Table of the outcomes of your code shows pretty much what the problem is. If

FuncThomae[x_] := If[
ExactNumberQ[Rationalize[x]],
If[x == 0,
1,
L = #^-1 & /@ Divisors[Numerator[Rationalize[x]]]
]
, 0]


then Table[FuncThomae[x], {x, 0, 1, 0.1}] produces

{1, {1}, {1}, {1, 1/3}, {1, 1/2}, {1}, {1, 1/3},
{1, 1/7}, {1, 1/2, 1/4}, {1, 1/3, 1/9}, {1}}


which makes it clear that the function is not producing numbers as output like it should.

If you want a working functional version of the Thomae function, the naive try

FuncThomae[x_] := If[
ExactNumberQ[Rationalize[x]],
If[x == 0,
1,
1/Denominator[Rationalize[x]]
]
, 0
]


does work, and you do not need to worry about eliminating common factors with the numerator before asking for the Denominator because the latter does the procedure as standard. (If not, what unique number could it produce?) With that version, Table[FuncThomae[x], {x, 0, 1, 0.1}] produces

{1, 1/10, 1/5, 1/10, 1/5, 1/2, 1/5, 1/10, 1/5, 1/10, 1}


and if you graph it using, say, DiscretePlot[FuncThomae[x], {x, 0, 1, 0.001}, PlotRange -> Full, Axes -> False], you get

Note, however, that this plot is not correct, because any discrete scan is going to miss points like $f(1/3)=1/3$. To get those points, run a plot like

ListPlot[
Flatten[Table[
{p/q, FuncThomae[p/q]}
, {p, 0, 100}, {q, 1, 100}], 1]
, PlotRange -> {{0, 1}, {0, 1}}, Axes -> False
]


to produce