Jacobi method in Matlab

My jacobi method not working. I feel like I'm missing something. Currently my numbers that are result back are m = 1 and be = 1.6. I wonder if I'm breaking out of one of the loops to fast
How to run :
n=100
[a,b]=sparsesetup(100);
[m,be] = Jacobi(a,b,0.00000025)
should result with -- getting wrong values
m = 36
be = 4.5785 * 10^-7
function [m,be] = Jacobi(a,b,tol)
n = length(b)
d= diag(a); % extract diagonal of a[i,j]
r=a-diag(d);
x=zeros(n,1);
p=zeros(n,1)
c=zeros([n,1]);
e=zeros([n,1]);
n1 =0;
err =0;
relerr = 1;
while(1)
x=b(b-r*x) ./d;
err =abs(norm(x-p));
relerr = err/(norm(x) + eps);
p=x;
n1= n1+1
if(err)
break
end
end
xc=x;
m=n1;
for i=1 : n
for j=1 : n
xa(j) = 1;
c(i) = c(i) + a(i,j) * xa(j);
e(i) = e(i) + a(i,j) * xc(j);
end
end
for i =1 : n
dif(i) = abs(xa(i) - xc(i));
dif2(i) = abs(e(i) - c(i));
end
fe = max (dif, [], 2);
be = max(dif2, [],2);
end
// sparesetup method
function [a,b] = sparsesetup (n)
e = ones(n,1);
a= spdiags ([-e 3*e -e], -1:1, n,n); %Entries of a{i,j}
b = zeros (n,1); % Entries of rhs b
b(1) = 2;
b(n) = 2;
b(2:n-1)=1;
end

3 Comments

if(err)
That is a floating point comparison that will lead to a break whenever err is non-zero which will be nearly always. You should not be assuming that you can get an exact 0, because of round-off error. You should be testing for err < some threshhold.
I'm not sure what the threshold should be ?
x=b(b-r*x) ./d;
That is indexing. Is it possible that you want multiplication instead of indexing?

Sign in to comment.

Answers (0)

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Tags

Asked:

on 27 Feb 2021

Commented:

on 27 Feb 2021

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!