clear all;
fs = 25;

% problem #2a
n = 10;
kmax = 10;
A = [4 3 0;3 4 -2;0 -2 4];
b = [24 30 -24]';
x = A\b;
x0 = zeros(3,1);
L = tril(-A,-1);
U = triu(-A,1);
D = diag(diag(A));
[xjac,iter] = jacobi(L,U,D,b,x0,kmax);
xgs = gauss_seidel(L,U,D,b,x0,kmax);
Tj = D\(L+U);
p = sort(abs(eig(Tj)));
w = 2/(1+sqrt(1-p(1)^2));
xsor = sor(L,U,D,b,x0,kmax,w);
[rescg,xcg] = cg(A,b,x0,kmax,10^-2);
irange = [1:kmax];
clear ejac egs esor ecg
for i=irange,
    try;rjac(i) = norm(A*(x-xjac(:,i)));ejac(i) = norm(x-xjac(:,i),inf);
    catch;ejac(i)=norm(x);end
    try;rgs(i) = norm(A*(x-xgs(:,i)));egs(i) = norm(x-xgs(:,i),inf);
    catch;egs(i)=norm(x);end
    try;rsor(i) = norm(A*(x-xsor(:,i)));esor(i) = norm(x-xsor(:,i),inf);
    catch;esor(i)=norm(x);end
    try;rcg(i) = norm(A*(x-xcg(:,i)));ecg(i) = norm(x-xcg(:,i),inf);
    catch;ecg(i)=0;end
end
p1 = semilogy(irange,egs,'-*', irange,esor,'-o', ...
	irange,ecg,'-s',irange,norm(x)*cond(A)*rgs/norm(b),':', ...
	irange,norm(x)*cond(A)*rsor/norm(b),'--');
l1 = legend('Gauss-Seidel','SOR','CG','GS-bound','SOR-bound',4);
set(l1,'FontSize',fs);
axis([0 11 1e-15 1e3]);
xlabel('Iteration','FontSize',fs);
ylabel('Error','FontSize',fs);
title('Bound for Thm 7.27','FontSize',fs);
set(gca,'FontSize',fs);
print -dtiff '2a.tiff'

% problem 2b
kmax = 10;
A = [4 3 0;3 4 -4;0 -2 4];
b = [24 30 -24]';
x = A\b;
x0 = zeros(3,1);
L = tril(-A,-1);
U = triu(-A,1);
D = diag(diag(A));
[xjac,Tj,iter] = jacobi(L,U,D,b,x0,kmax);
[xgs,Tgs] = gauss_seidel(L,U,D,b,x0,kmax);
p = sort(abs(eig(Tj)));
w = 2/(1+sqrt(1-p(1)^2));
[xsor,Tw] = sor(L,U,D,b,x0,kmax,w);
[rescg,xcg] = cg(A,b,x0,kmax,10^-2);
irange = [1:kmax];
clear ejac egs esor ecg
for i=irange,
    try,ejac(i) = norm(x-xjac(:,i),inf);catch,ejac(i)=0;end
    try,egs(i) = norm(x-xgs(:,i),inf);catch,egs(i)=0;end
    try,esor(i) = norm(x-xsor(:,i),inf);catch,esor(i)=0;end
    try,ecg(i) = norm(x-xcg(:,i),inf);catch,ecg(i)=0;end
end
semilogy(irange,egs,'-*',irange,esor,'-o',irange,ecg,'-s');
l2 = legend('Gauss-Seidel','SOR','CG','GS-bound','SOR-bound',2);
set(l2,'FontSize',fs);
axis([0 11 1e1 1e13]);
xlabel('Iteration','FontSize',fs);
ylabel('Error','FontSize',fs);
set(gca,'FontSize',fs);
print -dtiff '2b.tiff'


% problem %3
for n = [10 100];
    kmax = 10;irange = [1:kmax];
    tall = linspace(0,1,n+1);
    t = tall(2:end-1);
    h = t(1)-t(2);
    A = -(diag(2*ones(n-1,1))-diag(ones(n-2,1),-1)-diag(ones(n-2,1),1));
    A = A./h^2;
    f = 10*sin(2*pi*t)';
    x = A\f;
    xall = [0;x;0];

    x0 = zeros(n-1,1);
    L = tril(-A,-1);
    U = triu(-A,1);
    D = diag(diag(A));

    [xjac,Tj] = jacobi(L,U,D,f,x0,kmax);
    [xgs,Tgs] = gauss_seidel(L,U,D,f,x0,kmax);
    p = sort(abs(eig(Tj)));
    w = 2/(1+sqrt(1-p(end)^2));
    [xsor,Tw] = sor(L,U,D,f,x0,kmax,w);
    [rescg,xcg] = cg(A,f,x0,kmax,10^-2);
    plot(t,xgs(:,end),'-*',t,xsor(:,end),'-o', ...
        t,xcg(:,end),'-s');
    l3 = legend('Gauss-Seidel','SOR','CG',2);set(l3,'FontSize',fs);
    set(gca,'FontSize',fs);
    print('-dtiff',['3a' num2str(n) '.tiff']);

    tall = linspace(0,1,n+1);
    t = tall(2:end-1)';
    x = -5*sin(2*pi*t)/2/pi^2;
    clear egs esor ecg
    for i=irange,
        try;rgs(i) = norm(A*(x-xgs(:,i)));egs(i) = norm(x-xgs(:,i),inf);
        catch;egs(i)=norm(x);end
        try;rsor(i) = norm(A*(x-xsor(:,i)));esor(i) = norm(x-xsor(:,i),inf);
        catch;esor(i)=norm(x);end
        try;rcg(i) = norm(A*(x-xcg(:,i)));ecg(i) = norm(x-xcg(:,i),inf);
        catch;ecg(i)=0;end
    end
    p3 = semilogy(irange,egs,'-*', irange,esor,'-o', ...
        irange,ecg,'-s',irange,norm(x)*cond(A)*rgs/norm(f),':', ...
        irange,norm(x)*cond(A)*rsor/norm(f),'--');
    l3 = legend('Gauss-Seidel','SOR','CG','GS-bound','SOR-bound',2);
    set(l3,'FontSize',fs);
    if n==10
        axis([0 11 7e-3 1e2]);
    else
        axis([0 11 1e-6 1e5]);
    end
    xlabel('Iteration','FontSize',fs);
    ylabel('Error','FontSize',fs);
    title('Bound from Thm 7.27','FontSize',fs)
    set(gca,'FontSize',fs);
    print('-dtiff',['3b' num2str(n) '.tiff']);
end