%% 1. Exact differentiation
clear
N = 200;
step = 1/(N-1);
T = [0:step:1]';
g = 3*T.^3 - 2*T.^2;

figure
subplot(2,3,1)
plot(T,g); axis([0 1 -.2 1.1])
title('g(t) = 3t^3 - 2t^2')
subplot(2,3,4)
plot(T(1:end-1),diff(g)/step,'r'); axis([0 1 -2 6])
title('g`(t)')

%% 2. Noisy case

err = randn(N,1)/300; % White Gaussian noise
g_err = g + err;

subplot(2,3,2)
plot(T,g_err); axis([0 1 -.2 1.1])
title(['g_e = g + noise. Rel. error = ',sprintf('%.3f',norm(err)/norm(g)),'.'])
subplot(2,3,5)
plot(T(1:end-1),diff(g_err)/step); axis([0 1 -2 6])
title('g_e` by direct differentiation')

%% 3. Inversion of the matrix and SVD

A = zeros(N,N);
for i = 1:N
    for l = 1:i
        A(i,l) = 1/N;
    end
end
freal = A\g;

subplot(2,3,3)
plot(T,A\g_err,T,freal); axis([0 1 -2 6])
title('g_e` by inverting A')

[U,S,V] = svd(A);
f = zeros(N,1);
for l = 1:N
    f = f + dot(g_err,U(:,l))*V(:,l)/S(l,l);
end

subplot(2,3,6)
plot(T,f,T,freal); axis([0 1 -2 6])
title('g_e` by inverting A with SVD')

%% 4. Regularization by truncation of the SVD

f = zeros(N,1); error = f;

figure
for l = 1:N
    f = f + dot(g_err,U(:,l))*V(:,l)/S(l,l);
    subplot(1,2,1)
    plot(T,f,T,freal); axis([0 1 -2 6])
    error(l) = norm(f-freal)/norm(freal);
    title(['l <= ',num2str(l),'. Rel. error = ',num2str(error(l))])
    subplot(1,2,2)
    plot(error(1:l)); axis([1 N 0 1])
    title('Relative error as a function of the number of \sigma_l')
    pause(1/sqrt(l))
end