function hw1_01

% Task (b)  
% solving the problem for h = 0.05;
  h = 0.05;
  Tol = 1.0e-8;
  N = round(1.0 / h) - 2;
  w = 1.3;
  
% assembly and solution of the algebraic system (SOR with w=1.3);  
  [x,sol,it, cn] = SolveDirichlet(N, Tol, w);
  
% plotting the numerical solution;  
  figure(1); surf(x, x, sol); title('Numerical Solution')
  shading interp
  
% for comparison, plotting the analytical solution  
  y = x;
  r = exp(x') .* sin(pi*x') * sin(pi*y).^2;
  figure(2); surf(x, x, r'); title('Analytical Solution');
  shading interp
  
% Task (c) and (d)
% analyzing solution errors for various h
  hh = [ 0.25 0.1 0.05 0.02 ];
  Nh = length(hh);
% errors compared at the center point of the domain  
  u0 = exp(0.5) .* sin(pi*0.5) * sin(pi*0.5).^2;
  for k=1:Nh
      
    buf = sprintf('  Solving with h=%f \n', hh(k));
    disp(buf);
    N = round(1.0 / hh(k)) - 2;
    [x,sol,it, cn] = SolveDirichlet(N, Tol, w);
    u1 = sol(N/2+1, N/2+1);
    Err(k,1) = hh(k);
    Err(k,2) = abs( (u1 - u0) / u0 );
    Err(k,3) = cn;
  end;

  figure(3); loglog(Err(:,1), Err(:,2), 'x-'); title('Solution error');
  figure(4); loglog(Err(:,1), Err(:,3), 'x-'); title('Condition Number');
    
   
function [x,mat,it_sor, cn] = SolveDirichlet(n, tol, w)

% calculating the numerical solution of

%         -laplacian operator of u = f in [0,1]x[0,1]

% 5point-star
% Dirichlet boundary conditions
% cartesian grid with element length 1/(n+2)
% lexicographic numbering

% f = exp(x) * (-sin(pi * x) - 2 * cos(pi * x) * pi + 5 * sin(pi * x) * pi ^ 2) * cos(pi * y) ^ 2 + exp(x) * (sin(pi * x) + 2 * cos(pi * x) * pi - 3 * sin(pi * x) * pi ^ 2);
% => the analytic solution u0 is:
%
%        u0=e^x * sin(pi*x) * sin(pi*y)^2


% the values of g are deduced of the known analytic solution u0
x0=linspace(0,1,n+3);
x = x0(2:n+2);

%**************  the RHS vector ****************

s1 = - exp(x) .* (sin(pi*x) + 2*cos(pi*x)*pi-3*sin(pi*x)* pi^2);
s2 = - exp(x) .* (-sin(pi*x) - 2*cos(pi*x)*pi+5*sin(pi*x)*pi^2);

%s= exp(x).*(sin(2*pi*x) + 4*cos(2*pi*x) - pi*20*sin(2*pi*x)*pi^2);
f=[];
for j=1:n+1
   f=[ f, s2 * (cos(pi*x(j)))^2 + s1 ];
end

% p will be the right side of the finite difference equation A*u=p
p=f;
p=p';


%***************  Matrix A **********************
e=ones((n+1)^2,1);
L=spdiags([-e 4*e -e],[-1 0 1],(n+1),(n+1));
A=L;
for j=1:n
   A=[A,spalloc(j*(n+1),n+1,0)
      spalloc(n+1,j*(n+1),0),L];
end
% the nodes above and beneath
A=A+spdiags([-e -e],[-n-1 n+1],(n+1)^2,(n+1)^2);
A=(n+2)^2*A;


%***************  the condition number **************
cn = condest(A);


%*************** solution u **************************

% the initial guess;
uig(1:length(p),1) = 0.0;
[u, it_sor, err_sor] = solveSOR(A, p, uig, 10000, tol, w);

if ( err_sor(it_sor) > tol )
  disp('   Solution not converged ....')
  disp('   Try allowing more iterations ... ')
end;  
u = u';

% transform the solution vector to a matrix
mat=[]; 
for j=0:n+1:n*(n+1)
   mat=[mat;u(j+1:j+n+1)];
end




% SOR solver
function [u,m,er]=solveSOR(A, f, u_s, m_max, tol , w)
  
  m=0;
  
  er = [];
  % C_SOR = 1/w*(D+wL)
	D = diag(diag(A));
	L = tril(A,-1);
	C_SOR = 1/w*(D+w*L);
  
  % Setze u gleich dem Startvektor
  u = u_s;

	% Berechne das erste Residuum
	r = f - A*u;
  
  % Schleife: Abbruchbedingung: Norm<tol oder Anzahl der
  % maximalen Iterationen erreicht.
  while ((norm(r)/norm(f)>tol) & (m<m_max))
    
    % Iterationsschritt:
    u = u + C_SOR\r;
		r = f - A*u;
    
    % Anzahl der Iterationen um 1 erhöhen:
    m = m + 1;
    er = [er; norm(r)/norm(f)];
  end
