%% Numerik D-MAVT, RW (Nipp), SS06
%
% This is the M-file solution to the fourth exercise sheet.
% To run, save this file to your account and ensure that the "Current directory"
% in MATLAB is the same as that where the file is saved.
%
% This file requires that model.m, rhs.m and sparsemodel.m are placed 
% in the same directory.
%
% Then to run the file type "s12" in the command window
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function s12


clear
disp(['Aufgabe 2: heatequation with source term']);
dx=0.01; N=1/dx; Ni=N-1;
tf=0.1;
dt=tf/100;
u=zeros(Ni,1); 
uu=[0;u;0];  		% including boundary values
xi=[dx:dx:1-dx]';	% only interior points
A=sparse(1/dx^2*(-2*diag(ones(Ni,1))+diag(ones(Ni-1,1),1)+diag(ones(Ni-1,1),-1)));
M=speye(size(A))-dt*A;
dtb=dt*(-cos(2*pi*xi));	% h*b
for t=0:dt:tf-dt
  u = M\(u+dtb);
  uu =[uu [0;u;0]];
end
% plot it
x=[0;xi;1];
t=[0:dt:tf]';
figure;
mesh(t,x,uu);
title('Heat transfer with source term');
xlabel('t'); ylabel('x'); zlabel('u(t_i,x_j)');


%%%%%%%%%%%%%%%%%%% Aufgabe 3, poisson mit CG %%%%%%%%

disp(['Aufgabe 3: CG for poisson equation']);
n=100;
ni=n-2;
tol=[1e-7];
A=-delsq(numgrid('S',n));
b=1/n^2*ones(length(A),1);
u0=ones(size(b));
[u,n,t,r]=cg(A,-b,u0,tol,1000);

%u=A\b;

u=reshape(u,ni,ni);


u(1,:)=0;  % boundary conditions
u(:,1)=0;
u(:,ni)=0;
u(ni,:)=0;
figure;
x=linspace(-1,1,length(u));y=x;
mesh(x,y,u);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%% CG routine form Serie 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Serie 4, A3
% a) implement cg-method
function [u,k,time,resid]=cg(A,b,u0,tol,maxIter)

% conjugate gradients to solve A*u+b=0 (or A*u=-b, mind the sign ;-))
% in:	A 	matrix of linear system
%	b	(negative) right hand side of system
% 	u0	initial guess
% 	tol 	tolerance: if norm(residual)<tol ==> stop
%	maxIter	maximal number of iterations
% out:	u	approximate solution
%	k	number of iterations needed
%	time	cpu time needed
%	resid	2-norm of residual for final approximativ solution
% auth: fab@math.ethz.ch

tic;		% start measuring time
% initialisation
number_iteration=0;
u=u0;
rm1=A*u+b; p=-rm1; n2m1=rm1'*rm1; % rm1=r^{k-1}, n2m1=<rm1,rm1>
for k=1:maxIter
  if k>=2
    e = n2m1/n2m2;
    p = -rm1+e*p;
  end
  z    = A*p;
  rho  = n2m1/(p'*z);
  u    = u + rho*p;
  rm2  = rm1;
  rm1  = rm1+rho*z;
  n2m2 = n2m1;
  n2m1 = rm1'*rm1;
  number_iteration = number_iteration+1;
  if (n2m1<tol*tol) break; end;	% converged!
end
time=toc;	% stop measuring time
if (k==maxIter) 
  error(['not converged after ',int2str(maxIter),' iterations']); 
end
disp(['Number of iteration for cg   ', int2str(  number_iteration)]);
resid=norm(rm1);
return;



