Cantor function in matlab

How can we create a recursive sequence in matlab that converges to the cantor function?

f_0(x) will be equal to x

Then, for every integer n ≥ 0, the next function fn+1(x) will be defined in terms of fn(x) as follows:

Let fn+1(x) = 1/2 × fn(3x), when 0 ≤ x ≤ 1/3 ;

Let fn+1(x) = 1/2, when 1/3 ≤ x ≤ 2/3 ;

Let fn+1(x) = 1/2 + 1/2 × fn(3 x − 2), when 2/3 ≤ x ≤ 1.

I can define a non recursive function pretty easily(even one that is piecewise), but how to do a recursive function in matlab?

Answer 1

[Torsten]

N = 20;
f{1} = @(x) x;
for i = 1:N-1
    f{i+1} = @(x)(0.5*f{i}(3*x).*(x<=1/3) + 0.5.*(x>1/3 & x<2/3) + (0.5+0.5*f{i}(3*x-2)).*(x>=2/3 & x <=1));
end
x = linspace(0,1,1000);
plot(x,f{N}(x))

[Vivaan]

Thank you very much, very instructive!

[Vivaan]

What exactly does .* do in this context?

[Torsten]

.* multiplies two arrays for the same size elementwise. Or what exactly is your question ?

[Vivaan]

but why do we need to multiply an array?

[Torsten]

If x is an array of values, 0.5*f{i}(3*x) and (x<=1/3) are both arrays. And those arrays have to be multiplied componentwise. That's what .* is for. If you are sure your function is called only for a single value for x, you can also use * instead.

Do you understand what (x<=1/3) returns ? It returns 0 (false) for x>1/3 and 1 (true) for x<=1/3.

Thus the three different terms in the definition of f{i+1} constitute the piecewise definition of f{i+1}.

Answer 2

[Cris LaPierre]

You might find these 3 videos from the Mastering Programming in MATLAB Coursera course helpful.

• Recursion Part 1
• Recursion Part 2
• Recursion Part 3

Answer 3

[John D'Errico]

To be honest, I'd probably be lazy, and just use the simple algorithm found in wikipedia.

https://en.wikipedia.org/wiki/Cantor_function

That is, expresses the number in ternary. Easy enough to do. The largest integer power of 3 that is less than flintmax is 3^33.

flintmax

ans = 9.0072e+15

3^33

ans = 5.5591e+15

That means you can get 33 ternary digits for any decimal number between 0 and 1. We can do it like this:

ternary = @(x) dec2base(round(3^33*x),3,33);
t = ternary(0.123)

t = '010022200000020120010222011101211'

We can verify the result, as:

format long g
dot(3.^(-1:-1:-33),t - '0')

ans = 
                     0.123

Now just use the algorithm shown on the wiki page. For example, we know that f_inf(1/4) = 1/3.

T = ternary(1/4)

T = '020202020202020202020202020202021'

ind = find(T == '1',1,'first')

ind = 
    33

T(ind + 1:end) = '0' % replace the digits after the first 1, with 0.

T = '020202020202020202020202020202021'

% replace all 2's with a 1 
T(T == '2') = '1' % replace all 2's with a 1.

T = '010101010101010101010101010101011'

% finally, represent the number in base 2.
format rat
dot(T - '0',2.^(-1:-1:-33))

ans = 
       1/3     

Unfortunately, the result will have only 33 binary bits in the final representation as I did it here.

Yes, this is probably homework, since nobody is going to be computing this for any normal reason. :) And that means you were instructed to do it recursively, using the supplied set of relations.

[Vivaan]

+1 but it is not homework. I dont go to any university. just trying to learn some matlab :)