%% Time-Varying Systems and Computations
% Klaus Diepold, Matthias Kissel, Sven Gronauer, Nehir Güzelkaya
%
%% Diagonal expansion of Toeplitz matrices
%% 1. Shift operator
% The shift operator Z\;has the function to shift all entries of a vector downward by one position. This downward shift can be interpreted as a temporal delay, that is,
% u_k =Z\;u_{k+1} ,\;\;\;\;\;k=-\infty ,\ldotp \ldotp \ldotp ,\infty.
% 1.  Create column vectors u and v, both of length m. Take for now m=5.

m = 5;
u = randn(m,1) 
v = randn(m,1)

% Create a causal shift-matrix Z that shifts a column vector u of length m down by one position and a row-vector v^{\top } of length m one position to the left. Note that empty positions in the finite dimensional vectors u and v will be filled up with 0.

Z = diag(ones(m-1,1),-1)

% Verify your results with the provided vectors u and v.

u
Z*u
v'
v'*Z

% Change the shift operator into a cyclic shift operator Z_c \;such that vector elements at the bottom of u or at the left end of v^{\top } will be reinserted at the top or at the right end, respectively.

Z_c = diag(ones(m-1,1),-1);
Z_c(1,m)=1 

% Verify your results with the provided vectors u and v.

u
Z_c*u
v'
v'*Z_c

% Compare and discuss the properties of Z and Z_c.
% 2. Now extend the shift operator such that it shifts matrices in south-east direction along the main diagonal as showed in the figure.
%
% Hint: Perform two shift operations in different directions.

A = randn(5,5);

% Verify your results with the provided matrix A, i.e. shift the matrix A south-east along the main diagonal by an appropriate application of shift operators (from both sides).

A
B = Z*A*Z'

%
% 2. Toeplitz Matrices and linear convolutions
% We consider the convolution of two vectors u and v. To do this we can use the convolution operation in Matlab.

t = conv(u,v);

% Alternatively, we can represent the convolution by a matrix-vector multiplication using Teoplitz matrix.

T = toeplitz([u;zeros(length(v)-1,1)],[u(1);zeros(length(v)-1,1)]);
t_conv = T*v;

% Finally, we can compare our solution with the Matlab operation.

norm(t-t_conv)^2

%
%% 3. Condition number of band matrices
% Please create the described band matrices and plot the condition number of the resulting matrices for different matrix sizes n.

n = 1:50;

% Example 1: We consider a lower bi-diagonal matrix with 1's on the main diagonal and -0.5's on the subdiagonal.
%
% First we take a small example of this matrix and check the invertibility of this small example.

K2 = -0.5; 
B_6 = eye(6) + K2*diag(ones(5,1),-1)
B_inv = inv(B_6)
cond(B_6)

% What can you say about the invertibility of the small example matrix?
% Now let's increase the size of the matrix step by step and check for invertibility by looking at the condition number.

conds = zeros(size(n,2),1);
for i=1:50
    B = eye(i) + diag(K2*ones(i-1,1),-1);
    conds(i,1) = cond(B);
end
semilogy(conds)
title('Example 1')
ylabel('Condition Number')
xlabel('Dimension of lower bi-diagonal matrix')

% Example 2: We consider a lower bi-diagonal matrix with 1's on the main diagonal and -2's on the subdiagonal.
%
% We again take a small example of this matrix and check the invertibility of this small example.

K1 = -2;
B_6 = eye(6) + K1*diag(ones(5,1),-1)
B_inv = inv(B_6)
cond(B_6)

% What can you say about the invertibility of the matrix B_6?
% Now let's increase the size of the matrix step by step and check again for invertibility by looking at the condition number.

conds = zeros(size(n,2),1);
for i=1:50
    B = eye(i) + K1*diag(ones(i-1,1),-1);
    conds(i,1) = cond(B);
end
semilogy(conds)
title('Example 2')
ylabel('Condition Number')
xlabel('Dimension of lower bi-diagonal matrix')

% What do you observe?
% Example 3: Finally, we mix the two previous cases and look at a lower bi-diagonal matrix with 1's on the  main diagonal, -2's on the first half of the subdiagonal and -0.5's on the second half of the subdiagonal.
%
% Again, we first look at a small example.

subdiagonal = [K1*ones(2,1);K2*ones(3,1)];
B_6 = eye(6) + diag(subdiagonal,-1)
B_inv = inv(B_6)
cond(B_6)

% What can you say about the invertibility of the small example matrix?
% Now let's inctease the size of the matrix step by step and check again for invertibility by looking at the condition number.

conds = zeros(size(n,2),1);
for i=1:50
    subdiag = K1*ones(i-1,1);
    subdiag(round(i/2):end,1) = K2;
    B = eye(i) + diag(subdiag,-1);
    conds(i,1) = cond(B);
end
semilogy(conds)
title('Example 3')
ylabel('Condition Number')
xlabel('Dimension of lower bi-diagonal matrix')

% What do you observe? Try to find an explanation for this behavior.