i was calculate using gauss jacobi method with 4x4 matrice, why does NaN appear? what should i do?
10 次查看(过去 30 天)
显示 更早的评论
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]);
0 个评论
回答(2 个)
Torsten
2023-10-29
编辑:Torsten
2023-11-2
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)
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
2023-10-29
编辑:John D'Errico
2023-10-29
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]
Next, you want to read about the 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 个评论
John D'Errico
2023-10-29
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]
B=[3;1;6;16];
P = [2 1 3 4]; % a row permutation. It must be applied to both A and B.
A(P,:)
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)
A\B
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!