i = 0;
for n = 35:5:60% number of inner nodes in one direction.
inew = i+1;
a_amp = 12; % amplitude for the function a(x_1,x_2)
f_amp = 50; % we can choose f=1, 50, 100
x_0=0.5;
y_0=0.5;
c_x=1;
c_y=1;
h = 1/(n+1); % define step length
% ----------------------------------------
% Computing all matrices and vectors
% ----------------------------------------
% Generate a n*n by n*n stiffness matrix
S = DiscretePoisson2D(n);
%% generate coefficient matrix of a((x_1)_i,(x_2)_j) = a(i*h,j*h)
C = zeros(n,n);
for i=1:n
for j=1:n
C(i,j) = 1 + a_amp*exp(-((i*h-x_0)^2/(2*c_x^2)...
+(j*h-y_0)^2/(2*c_y^2)));
end
end
% create diagonal matrix from C
D = zeros(n^2,n^2);
for i=1:n
for j=1:n
D(j+n*(i-1),j+n*(i-1)) = C(i,j);
end
end
% If f is constant.
% f = f_amp*ones(n^2,1);
% If f is Gaussian function.
f=zeros(n^2,1);
for i=1:n
for j=1:n
f(n*(i-1)+j)=f_amp*exp(-((i*h-x_0)^2/(2*c_x^2)...
+(j*h-y_0)^2/(2*c_y^2)));
end
end
% Compute vector of right hand side
% b = D^(-1)*f computed as b(i,j)=f(i,j)/a(i,j)
b=zeros(n^2,1);
for i=1:n
for j=1:n
b(n*(i-1)+j)=f(n*(i-1)+j)/C(i,j); % Use coefficient matrix C or
% diagonal matrix D to get a(i,j)
end
end
% ----------------------------------------
% Solution of 1/h^2 S u = b using iterative Gauss-Seidel method
% with red-black ordering, version II
% ----------------------------------------
err = 1; k=0; tol=10^(-9);
% Initial guess
uold = zeros(n+2, n+2);
unew= uold;
tic
while(err > tol)
% Red nodes
for i = 2:n+1
for j = 2:n+1
if(mod(i+j,2) == 0)
unew(i, j) = (uold(i-1, j) + uold(i+1, j) + uold(i, j-1) + uold(i, j+1)+ h^2*b(n*(i-2)+j-1))/4.0;
% for computation of residual
u(j-1 + n*(i-2)) = unew(i,j);
end
end
end
% Black nodes
for i = 2:n+1
for j = 2:n+1
if(mod(i+j,2) == 1)
unew(i,j) = 0.25*(unew(i-1,j) + unew(i+1,j) ...
+ unew(i,j-1) + unew(i,j+1) + h^2*b(n*(i-2)+j-1));
% for computation of residual
u(j-1 + n*(i-2)) = unew(i,j);
end
end
end
k = k+1;
% different stopping rules
err = norm(unew-uold);
%computation of residual
% err = norm(S*u' - h^2*b);
uold = unew;
end
toc
s = toc;
u = reshape(unew(2:end-1, 2:end-1)', n*n, 1);
disp('-- Number of iterations in the version II of Gauss-Seidel method----------')
k
hold on;
plot(s(inew), n(inew)^2,'-x');
title('Time vs. Matrix Size')
xlabel('Time')
ylabel('Matrix Size')
end
