% Matlab code for the simulation of linear averaging algorithm.
xinitial = [+1 -1 +1 -1 +1]';
time_steps = 50;

% Complete graph
Acomplete = 1/5*ones(5,5);
xcompletefinal = Acomplete^(time_steps) * xinitial
% result: xcompletefinal is 0.2 [1; 1; 1; 1; 1]

% Ring topology
Aring = [1/3 1/3 0 0 1/3;
1/3 1/3 1/3 0 0;
0 1/3 1/3 1/3 0 ;
0 0 1/3 1/3 1/3;
1/3 0 0 1/3 1/3];
xringfinal = Aring^(time_steps) * xinitial
% result: xringfinal is 0.2 [1; 1; 1; 1; 1]

% Star topology. center = node 1
Astar = [1/5 1/5 1/5 1/5 1/5;
1/2 1/2 0 0 0;
1/2 0 1/2 0 0;
1/2 0 0 1/2 0
1/2 0 0 0 1/2];
xstarfinal = Astar^(time_steps) * xinitial
% result: xstarfinal is 0.3846 [1; 1; 1; 1; 1]


%% Alternative method Seen in lectures 
% Compute the left eigenvector of A associated to lambda=1
% by computing the corresponding right eigenvector of A'

[Vcomplete,Dcomplete]=eig(Acomplete');
d=diag(Dcomplete);
% find the index of the eigenvalue 1. Since A is primitive, 
% this is the largest eigenvalue. This justifies
% the following line.
[y,k]=max(d);
Vcompletenorm=Vcomplete(:,k)/sum(Vcomplete(:,k));
xcompletefinal1=Vcompletenorm'*xinitial*(ones(5,1))

[Vring,Dring]=eig(Aring');
d=diag(Dring);
[y,k]=max(d);
Vringnorm=Vring(:,k)/sum(Vring(:,k));
xringfinal1=Vringnorm'*xinitial*(ones(5,1))


[Vstar,Dstar]=eig(Astar');
d=diag(Dstar);
[y,k]=max(d);
Vstarnorm=Vstar(:,k)/sum(Vstar(:,k));
xstarfinal1=Vstarnorm'*xinitial*(ones(5,1))