%% Nowhere differentiable function

%%
% The building block is the periodic function g, defined to be |x| on
% [-1,1) and extended to be 2-periodic. 

%%
% I'll define it on the positive reals in MATLAB using frac, the fractional
% part of a number, and the Heaviside function, or unit step function,
% which is 0 on negative numbers and 1 on positive numbers.

syms x k n
g1 = 2*frac((x/2))
g2 = 2*(1-frac(x/2))
g = heaviside(g2-g1)*g1 + heaviside(g1-g2)*g2


%%
% I graph g to make sure it's right.
ezplot(g,[0,6])
%%
% The function f is defined by:

f = symsum((3/4)^k*subs(g,x,4^k*x),k, 0,n)

%%
% Here's a plot showing g(x) (blue), g(4x) (red) and g(4^2(x)) (green)

X = 0:0.001:2;
G = @(x) eval(vectorize(g))
g1 = subs(g,x,4*x)
G1 = @(x) eval(vectorize(g1))
g2 = subs(g,x,4^2*x)
G2 = @(x) eval(vectorize(g2))
plot(X,G(X)); hold on
plot(X,G1(X),'color','r')
plot(X,G2(X),'color','g')
hold off
axis equal
%%
% Here are the plots of the partial sums for n=1,2,3,4,5.

ezplot(subs(f,n,1),[0,1])

%%

ezplot(subs(f,n,2),[0,1])

%%

ezplot(subs(f,n,3),[0,1])

%%
X = 0:0.0001:1;
F = @(x) eval(vectorize(subs(f,n,4)))
plot(X,F(X))

%%
F = @(x) eval(vectorize(subs(f,n,5)))
plot(X,F(X))