Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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:

plot of Thomae's function

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

erroneous plot

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

Wanted to duplicate this

share|improve this question
    
Hi there, is the code you supplied complete, I didn't seem to be able to get it to run successfully ? –  image_doctor May 13 '12 at 11:38
    
@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 May 13 '12 at 12:42
    
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 May 13 '12 at 17:31
    
A tiny reminder: in floating-point arithmetic (which is what Mathematica internally uses when plotting), all numbers are rational. –  J. M. May 14 '12 at 11:15

3 Answers 3

up vote 16 down vote accepted

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}]

Thomae's function

share|improve this answer
1  
If you define pq as pq = {#, 1/Denominator[#]}& /@ fracs you could ListPlot pq directly without having to transform it first. –  Heike May 13 '12 at 12:39
    
@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 May 13 '12 at 13:19
    
btw could you do this for -1 to 1 –  The-Ever-Kid May 13 '12 at 13:20
    
@The-Ever-Kid yes that bit is straightforward, so I didn't bother implementing it –  acl May 13 '12 at 13:20
    
yeah cool thanks... –  The-Ever-Kid May 13 '12 at 13:23

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}]

Thomae's function applied to the Farey sequence

share|improve this answer

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$.

Thomae function

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.