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 个评论

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?

请先登录,再进行评论。

回答(0 个)

类别

帮助中心File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by