-
Notifications
You must be signed in to change notification settings - Fork 0
/
iterativecodes.m
112 lines (107 loc) · 2.81 KB
/
iterativecodes.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
clear; close all;
format long;
N=10; % size of matrix
% Set iteration tolerance and max num of iterations here:
%TOL = input('enter tolerance: ');
TOL = 1e-12;
Nmax = 5000;
% load matrix problem
% Ax = B
%
U = randn(N);
A=U*U'+150*eye(N); % matrix is diagonally dominant for small N
rr=cond(A);
fprintf('Condition number of matrix = %f \n \n',rr);
x=randn(N,1); % solution
B = A*x; % RHS
disp('PROGRAM WILL FIRST SOLVE USING JACOBI');
disp('THEN IT WILL SOLVE USING GAUSS-SEIDEL OR SOR');
counts=0;
xp=zeros(N,1);xo=xp;
% SOR AND GAUSS SEIDEL:
% initial guess assumed to be the null vector:
clear error;
%
% Successive Over-Relaxation Technique:
% Gauss-Seidel: same as SOR but w = 1
% set w=1 in SOR for Gauss-Seidel:
xp = zeros(N,1);
counts = 0;
w = input('Enter 1 for G-S. Enter SOR factor 1 <= w < 2: ')
% uncomment to test:
%w = 1.2;
tic;
for k=1:Nmax
counts = counts + 1;
for i=1:N
s =0;
for j=1:N
if i ~= j, s = s + A(i,j)*xp(j,1); end
end
xp(i,1) = w*(-s + B(i,1))/A(i,i) + (1-w)*xp(i,1);
end
error(k) = norm(A*xp-B,inf);
% fprintf(' It. Num. = %3.0f, error = %15.13e\n', ...
% k, error(k))
if error(k) < TOL, break; end
if k == Nmax,
fprintf('\n MAXIMUM NUMBER OF ITERATIONS EXCEEDED\n'),
end;
end
tgs = toc;
fprintf(' \n \n time=%f, It. Num. = %3.0f, error = %15.13e\n\n\n', ...
tgs, k, error(counts))
norm(x-xp)
semilogy(1:counts,error(1:counts));
grid on
xlabel('iteration number')
ylabel('error')
if w==1, title('GAUSS-SEIDEL ITERATION TECHNIQUE'), end;
if w~=1
omeg=num2str(w);
title(['SOR ITERATION WITH OMEGA =', omeg])
end;
disp('PRESS ANY KEY TO CONTINUE')
pause
clear error;
%
% Successive Over-Relaxation Technique:
% Gauss-Seidel: same as SOR but w = 1
% set w=1 in SOR for Gauss-Seidel:
xp = zeros(N,1);
counts = 0;
w = input('Enter 1 for G-S. Enter SOR factor 1 <= w < 2: ')
% uncomment to test:
%w = 1.2;
tic;
for k=1:Nmax
counts = counts + 1;
for i=1:N
s =0;
for j=1:N
if i ~= j, s = s + A(i,j)*xp(j,1); end
end
xp(i,1) = w*(-s + B(i,1))/A(i,i) + (1-w)*xp(i,1);
end
error(k) = norm(A*xp-B,inf);
% fprintf(' It. Num. = %3.0f, error = %15.13e\n', ...
% k, error(k))
if error(k) < TOL, break; end
if k == Nmax,
fprintf('\n MAXIMUM NUMBER OF ITERATIONS EXCEEDED\n'),
end;
end
tgs = toc;
fprintf(' \n \n time=%f, It. Num. = %3.0f, error = %15.13e\n\n\n', ...
tgs, k, error(counts))
norm(x-xp)
semilogy(1:counts,error(1:counts));
grid on
xlabel('iteration number')
ylabel('error')
if w==1, title('GAUSS-SEIDEL ITERATION TECHNIQUE'), end;
if w~=1
omeg=num2str(w);
title(['SOR ITERATION WITH OMEGA =', omeg])
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%5