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 comentarios

Walter Roberson
Walter Roberson el 27 de Feb. de 2021
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.
gracias claude
gracias claude el 27 de Feb. de 2021
I'm not sure what the threshold should be ?
Walter Roberson
Walter Roberson el 27 de Feb. de 2021
x=b(b-r*x) ./d;
That is indexing. Is it possible that you want multiplication instead of indexing?

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Loops and Conditional Statements en Centro de ayuda y File Exchange.

Etiquetas

Preguntada:

el 27 de Feb. de 2021

Comentada:

el 27 de Feb. de 2021

Community Treasure Hunt

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

Start Hunting!

Translated by