%% Week 6
% Matlab codes for week

%% Problem 1

% Gradient method

A=[4,2;2,3]; b=[7;10];
maxiter=20;

[n,n]=size(A);
iter=0; r=b; r0=norm(r); err=r0; x=zeros(n,1);
i=0;
relative_res_error=1;
while i<maxiter && relative_res_error>10^-8
i=i+1
s=A*r;
alpha=(r'*r)/(r'*s);
x=x+alpha*r
pause
r=r-alpha*s;
relative_res_error=norm(r)/r0
end

%%
% Conjugate gradient method

A=[4,2;2,3]; b=[7;10];
maxiter=20;

[n,n]=size(A);
iter=0; r=b; p=r; r0=norm(r); err=r0; x=zeros(n,1);
i=0;
relative_res_error=1;
while i<maxiter && relative_res_error>10^-8
i=i+1
s=A*p;
alpha=(r'*r)/(p'*s);
x=x+alpha*p
pause
rnew=r-alpha*s;
beta=rnew'*rnew/(r'*r);
p=rnew+beta*p;
r=rnew;
relative_res_error=norm(r)/r0
end

%% Problem 4

A = delsq(numgrid('S',100));
b = ones(size(A,1),1);

tic, [x,flag,final_relative_residual,iter,vector_of_res_errors]=pcg(A,b,10^-8,100); toc
flag
final_relative_residual
plot(0:iter,log(vector_of_res_errors))

%% Problem 5

x=[2;6;-3];
v=x; v(1)=v(1)+norm(v);
H=eye(3)-2*v*v'/(v'*v)
H*x

c=6/norm(x(2:3)); s=3/norm(x(2:3));
G=blkdiag(1,[c,-s;s,c])
G*x

%% Problem 6

A=[0,0;1,3;0,2]
x=A(1:2,1);
c=x(1)/norm(x), s=-x(2)/norm(x)
G1=blkdiag([c,-s;s,c],1) % The first Givens matrix
A1=G1*A
x=A1(2:3,2);
c=x(1)/norm(x), s=-x(2)/norm(x)
G2=blkdiag(1,[c,-s;s,c]) % The second Givens matrix
R=G2*A1 % R
Q=G1'*G2' % Q
Q*R % Check whether A=QR 
qr(A)

%% Problem 8

A=[0,0;1,3;0,2] % the matrix from the previous problem
b=[3;4;1];

% mldivide

A\b

% normal equation

(A'*A)\(A'*b)

% QR

[Q,R]=qr(A);
c=Q'*b;
R(1:2,1:2)\c(1:2)

%% Problem 9

A=[1,1,1;
   0,0,1;
   1,1,1;
   4,2,1]

b=[1;0;3;2]

A\b %-x^2+3x

% with polynomial fitting

z=polyfit([1,0,1,2],[1,0,3,2],2)



