function fig10_20

% modal approximation

% random generation of A values
% m=10;   % number of vertices (must be even)
% A(1:m/2,1:m/2)=10+randi(20,m/2);
% A(1+m/2:m,1+m/2:m)=10+randi(20,m/2);
% A(1+m/2,1)=1;
% A=A./sum(A);

% book example
m=4; A=zeros(m,m);
d=1/20;
A(1,1)=1/4-d; A(2,1)=3/4-d; A(3,1)=2*d;
A(1,2)=1/3; A(2,2)=2/3;
A(3,3)=1/4; A(4,3)=3/4;
A(3,4)=1/3; A(4,4)=2/3;

% sum(A)
A
eigens=eig(A)
sings=svd(A)

N=m+1;  % number of time steps
y=zeros(N,m);
y(1,1)=100;
%y(1,:)=ones(1,m);
for i=2:N
    z=A*y(i-1,:)';
    y(i,:)=z';
end
%q=0.1*(2*rand(1,m)-1);
q=0;
yd=(1+q).*y;
XD=yd(1:N-1,:);
YD=yd(2:N,:);
n=N-1;

% use the pseudo-inverse
%     X=XD;
%     Y=YD;
%     G=X'*X;
%     GI=inv(G);
%     kappa(k)=cond(G,inf)
%     P=Y'*X*GI;

% use the SVD
[U,S,V]=svd(XD);
SS=zeros(n,m);
for i=1:m
    SS(i,i)=1/S(i,i);
end
P=YD'*U*SS*V';
Ac=P;
error=norm(A-Ac,inf)/norm(A,inf)

[VP,DP]=eig(P);
[D0,I] = sort(abs(diag(DP)),'descend');
for i=1:m
    D(i)=DP(I(i),I(i));
end
D
VP = VP(:, I);
x=XD(1,:)';
alpha=VP\x
for it=1:10
    x=Ac*x;
    x1=alpha(1)*D(1)^it*VP(:,1);
    x2=x1+alpha(2)*D(2)^it*VP(:,2);
    x3=x2+alpha(3)*D(3)^it*VP(:,3);
    x4=x3+alpha(4)*D(4)^it*VP(:,4);
    err1(it)=norm(x-x1,inf)/norm(x);
    err2(it)=norm(x-x2,inf)/norm(x);
    err3(it)=norm(x-x3,inf)/norm(x);
    err4(it)=norm(x-x4,inf)/norm(x);
end
err4_error=norm(err4,inf)

figure(1)
clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.06,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
semilogy(1:10,err1,'--sk','LineWidth',1.5,'MarkerSize',8)
hold on
semilogy(1:10,err2,'--*r','LineWidth',1.5,'MarkerSize',8)
semilogy(1:10,err3,'--ob','LineWidth',1.5,'MarkerSize',8)
box on
grid on
axis([1 10 1e-10 5])
%set(gca,'ytick',[0 0.05 0.1])
xlabel('t-axis')
ylabel('Relative Error')
legend({' M = 1 ',' M = 2 ',' M = 3 '},'AutoUpdate','off','Location','SouthWest','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',16,'FontWeight','bold')

%exportgraphics(gca,'/Users/mark/Desktop/cars2.eps')





