i was calculate using gauss jacobi method with 4x4 matrice, why does NaN appear? what should i do?

12 次查看（过去 30 天）
Beatrice
Beatrice 2023-10-29
编辑： Beatrice 2023-11-2
function[Xi,g,H]=jacobi(A,B,X,N)
clc;
A=[1 3 1 -1; 2 0 1 1; 0 -1 4 1; 0 1 1 -5]
B=[3;1;6;16];
X=[0;0;0;0];
Xi=inv(A)*B;
N=10;
abcd=X';
n=4;
X1=X;
for k=1:N,
for i=1:n,
S=B(i)-A(i,[1:i-1,i+1:n])*X([1:i-1,i+1:n]);
X1(i)=S/A(i,i);
end
g=abs(X1(i)-X(i))/X1(i)*100;
eror=norm(g);
X=X1;
abcd=[abcd;X'];
end
disp('abcd');
disp([abcd]);
disp('eror');
disp([eror]);

请先登录，再回答此问题。

回答（2 个）

Torsten
Torsten 2023-10-29
编辑：Torsten 2023-11-2
Ran in:
A(2,2) equals 0, but you divide by 0.
What is H ?
clc;
A=[1 3 1 -1; 2 0 1 1; 0 -1 4 1; 0 1 1 -5];
A=[A(2,:);A(1,:);A(3,:);A(4,:)];
B=[3;1;6;16];
B=[B(2);B(1);B(3);B(4)];
X=[0;0;0;0];
N=20;
[Xi,g,H]=jacobi(A,B,X,N)
abcd 0 0 0 0 0.5000 1.0000 1.5000 -3.2000 1.3500 -0.7333 2.5500 -2.7000 0.5750 -1.2000 1.9917 -2.8367 0.9225 -0.8011 1.9092 -3.0417 1.0663 -0.9578 2.0601 -2.9784 0.9591 -1.0349 2.0052 -2.9795 0.9872 -0.9813 1.9862 -3.0060 1.0099 -0.9931 2.0062 -2.9990 0.9964 -1.0050 2.0015 -2.9974 0.9980 -0.9984 1.9981 -3.0007 1.0013 -0.9989 2.0006 -3.0001 0.9997 -1.0006 2.0003 -2.9997 0.9997 -0.9999 1.9998 -3.0001 1.0002 -0.9998 2.0000 -3.0000 1.0000 -1.0001 2.0000 -3.0000 1.0000 -1.0000 2.0000 -3.0000 1.0000 -1.0000 2.0000 -3.0000 1.0000 -1.0000 2.0000 -3.0000 1.0000 -1.0000 2.0000 -3.0000 1.0000 -1.0000 2.0000 -3.0000 eror 2.3864e-05
Output argument "H" (and possibly others) not assigned a value in the execution with "solution>jacobi" function.
function[Xi,g,H]=jacobi(A,B,X,N)
Xi=inv(A)*B;
abcd=X';
n=4;
X1=X;
for k=1:N,
for i=1:n,
S=B(i)-A(i,[1:i-1,i+1:n])*X([1:i-1,i+1:n]);
X1(i)=S/A(i,i);
end
g=abs(X1(i)-X(i))/X1(i)*100;
eror=norm(g);
X=X1;
abcd=[abcd;X'];
end
disp('abcd');
disp([abcd]);
disp('eror');
disp([eror]);
end
John D'Errico
John D'Errico 2023-10-29
编辑：John D'Errico 2023-10-29
Ran in:
You write a function that allows you to pass in A and B, but then you define A anb B in the function? Sigh. I think you do not understand why functions exist. Anyway...
A=[1 3 1 -1; 2 0 1 1; 0 -1 4 1; 0 1 1 -5]
A = 4×4
1 3 1 -1 2 0 1 1 0 -1 4 1 0 1 1 -5
Next, you want to read about the Jacobi method.
https://en.wikipedia.org/wiki/Jacobi_method
Somewhere along the way, your teacher should have said the Jacobi method ONLY applies to strictly diagonally dominant matrices. We can find this statement from that Wiki link I have copied here.
"A sufficient (but not necessary) condition for the method to converge is that the matrix A is strictly or irreducibly diagonally dominant. Strict row diagonal dominance means that for each row, the absolute value of the diagonal term is greater than the sum of absolute values of other terms:"
What does that mean? A matrix with any zero element on the diagonal can NEVER be convergent using Jacobi iteration. That is because you will be always dividing by zero. Those Nans are evidence of what happens.
So what should you do? Get a different matrix. Or, get a different algorithm, since Jacobi cannot be used to solve this problem.
  5 个评论
显示 3更早的评论
John D'Errico
John D'Errico 2023-10-29
Ran in:
@Torsten
Even on the actual problem, it is not strictly diagonally dominant, but yes, it is sufficient, on that specific problem to use a row permutation.
The problem is when you suggest that row permutations are sufficient to solve a problem, but NOT explain enough about the issues that you will now cause further issues down the line. Others might actually believe you, and think you have made a valid general claim because you have a good reputation on the site.
A=[1 3 1 -1; 2 0 1 1; 0 -1 4 1; 0 1 1 -5]
A = 4×4
1 3 1 -1 2 0 1 1 0 -1 4 1 0 1 1 -5
B=[3;1;6;16];
P = [2 1 3 4]; % a row permutation. It must be applied to both A and B.
A(P,:)
ans = 4×4
2 0 1 1 1 3 1 -1 0 -1 4 1 0 1 1 -5
A you can see, now the diagonal elements are >= the sum of absolute values of the off diagonals.
[X,iter,del] = jaciter(A(P,:),B(P),1e-14,1000)
X = 4×1
1.0000 -1.0000 2.0000 -3.0000
iter = 50
del = 7.9813e-15
A\B
ans = 4×1
1.0000 -1.0000 2.0000 -3.0000
function [X,iter,del] = jaciter(A,b,tol,itmax)
% uses Jacobi iteration to solve the problem A*X==B.
%
% A: NxN array, presumed non-singular
% b: Nx1 vector
% tol: absolute tolerance for declaration of convergence
% itmax: maximum number of iterations allowed
%
% Note, the sizes of A and b are not checked for compatibility,
% nor is a check made for diagonal dominance.
n = size(A,1);
D = diag(A);
b = b(:);
LU = A - diag(D); % only the off-diagonal elements of A
Dinv = 1./D; % leave Dinv as a vector.
X0 = rand(n,1); % random start point for the iteration
del = inf;
iter = 0;
while (del>tol) && (iter <= itmax)
iter = iter + 1;
if iter >= itmax
warning('Maximum iterations exceeded')
end
% Note the use of .* below. This will cause an element-wise
% multiply to form X
X = Dinv.*(b - LU*X0);
del = norm(X - X0);
X0 = X;
end
if del > tol
warning('Convergence tolerance was not met at termination')
end
end
Beatrice
Beatrice 2023-11-2
编辑：Beatrice 2023-11-2
thank you for you're correction. Actually i get this matrix problem somewhre in a book called metode numerik by bambang triatmodjo page 85 (you can check it).
I do find it strange when i see this matrix, because one of the conditions for using the Jacobi method is that the diagonal side of the matrix must be larger than the surrounding columns/rows. But I think that there is another solution to this problem. Once again, thank you for your solution. I do learn more from you, and i hope you can help me again in the future =)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by