-
Notifications
You must be signed in to change notification settings - Fork 0
/
jacobi.m
63 lines (51 loc) · 1.33 KB
/
jacobi.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
clear
n = 1000;
randmat = sprand(1000,1000,.005);
A=randmat+10*speye(1000);
n = size(A,1);
b = rand(n,1);
x = zeros(n,1);
D = spdiags(diag(A),0,n,n);
R = A - D;
iter = 1;
r = b-A*x;
rnorm(iter) = norm(r);
while(rnorm(iter)/rnorm(1) > 1e-6 && iter < 100)
x=D\(b-R*x);
r = b-A*x;
iter = iter + 1;
rnorm(iter) = norm(r);
fprintf('Error at iteration %i: %f\n', iter-1, rnorm(iter))
end
A=randmat+5*speye(1000);
x = zeros(n,1);
D = spdiags(diag(A),0,n,n);
R = A - D;
iter = 1;
r = b-A*x;
rnorm2(iter) = norm(r);
while(rnorm2(iter)/rnorm2(1) > 1e-6 && iter < 100)
x=D\(b-R*x);
r = b-A*x;
iter = iter + 1;
rnorm2(iter) = norm(r);
fprintf('Error at iteration %i: %f\n', iter-1, rnorm2(iter))
end
A=randmat+2.5*speye(1000);
x = zeros(n,1);
D = spdiags(diag(A),0,n,n);
R = A - D;
iter = 1;
r = b-A*x;
rnorm3(iter) = norm(r);
while(rnorm3(iter)/rnorm3(1) > 1e-6 && iter < 100)
x=D\(b-R*x);
r = b-A*x;
iter = iter + 1;
rnorm3(iter) = norm(r);
fprintf('Error at iteration %i: %f\n', iter-1, rnorm3(iter))
end
semilogy(0:size(rnorm,2)-1,rnorm/rnorm(1), 0:size(rnorm2,2)-1,rnorm2/rnorm2(1), 0:size(rnorm3,2)-1,rnorm3/rnorm3(1))
xlabel('iteration')
ylabel('relative residual')
legend('big diag', 'medium diag', 'little diag')