Jacobi iterative method problem
显示 更早的评论
I've implemented the Jacobi method in matlab but when i try it , the function give me wrongs results.. I don't know what I'm wrong. I post here my code and results.
function [x,niter,resrel] = jacobiSolver(a, b, varargin)
narginchk(2, 5); % check if the function receives the right number of input parameters
nargoutchk(0,3); % check if the function receives the right number of output parameters
checkAB(a,b) % priority control
% OPTIONAL PARAMETER CONTROL
switch nargin
case 4; [TOL, NMAX] = CheckNMAX(varargin{1}, varargin{2});
case 3; [TOL, NMAX] = CheckTOL(varargin{1}, 500);
case 2; TOL = 10^-6; NMAX = 500;
otherwise
end
% codifichiamo l'algoritmo => vedere se dobbiamo sostituire con spdiags
c = b./diag(a);
B = (-1./diag(a)')' .*(a-diag(diag(a)));
x0=zeros(length(b),1);
x=c;
nit = 0;
%TOLX = max(TOL*norm(x,1),realmin);
TOLX = TOL;
while norm(x-x0,1)>=TOLX && nit<NMAX
x0 = x ;
x = B*x0 + c;
TOLX = max(realmin, TOL*norm(x0,1));
nit = nit+1
end
% IF THE NMAX VALUE IS REACHED THEN THE USER WILL BE NOTIFIED WITH A WARNING
if(nit==NMAX)
warning("Limit reached, probably not convergence... calculation stopped and resrel = "+0);
end
if nargout > 1
niter = nit;
if nargout == 3
resrel = 0;
end
end
end
results obteined with simple matrix a and b:
Warning: Limit reached, probably not convergence... calculation stopped and resrel = 0
> In jacobiSolver (line 32)
ans =
1.0e+295 *
0.3636
1.2137
0.7255
0.9811
1.0341
correct results:
ans =
20.2162
-43.2162
-56.8378
102.1081
-21.7838
also are there better ways to implement this algorithm?
thanks a lot for helping!
7 个评论
Walter Roberson
2019-4-24
Please post sample a and b and optional arguments for us to test with.
Agostino Dorano
2019-4-24
Walter Roberson
2019-4-24
checkAB() ?
Agostino Dorano
2019-4-24
Walter Roberson
2019-4-24
编辑：Walter Roberson
2019-4-25
You calculate a mostly negative B and a positive c, and you have a recursion like
f(n) = B*f(n-1) + c
and you stop when norm(f(n)-f(n-1),1) < tolerance
but f(n) - f(n-1) is going to be B*f(n-1)+c - f(n-1) = (B-1)*f(n-1) + c and the 1 norm is the sum of the absolute value of the components, so sum( abs((B(K)-1)*x0(K) + c(K)) ) < tolerance . An x0 that would be true for would be x0 = -c./(B-1) or close by.
Anyhow, you start with an x0 and head for smaller and smaller x0 looking for one that will have the B vs c balance. However, your code with its 1-norm cannot tell which side of the balance you are on, so when you overshoot the balance, you just keep going in the same direction, getting further and further from the balance.
valerio auricchio
2019-4-25
Sorry Walter Roberson can you explain bettere this "An x0 that would be true for would be x0 = -c./B or close by.
Anyhow, you start with an x0 and head for smaller and smaller x0 looking for one that will have the B vs c balance. However, your code with its 1-norm cannot tell which side of the balance you are on, so when you overshoot the balance, you just keep going in the same direction, getting further and further from the balance." I don't think I caught well.
Walter Roberson
2019-4-25
If x0 exactly equaled -c./(B-1) then each abs((B(K)-1)*x0(K)+c(K)) component would be 0, and the sum would be 0. But you do not need an exact 0, so the x0 values can potentially be close to -c./(B-1) and still have the 1-norm less than the tolerance.
(Note: I made a correction for the -1)
采纳的回答
John D'Errico
2019-4-25
编辑：John D'Errico
2019-4-25
Note that the Jacobi method for solution of a linear system is only convergent for SOME matrices. It is NOT always convergent. You might want to do some reading, here:
The point is, if your system is not diagonally dominant, then the Jacobi iteration will not converge. It will diverge. You are trying this on a RANDOM matrix. But a divergent result tends to imply divergence has occurred. Digging down into your comments, I see this matrix. Sigh. I wonder if it is diagonally dominant? Any takers on that bet?
a = [5 1 1 1 4
5 2 5 3 1
1 3 5 5 5
5 5 3 4 5
4 5 5 5 4];
b =[76
75
40
66
18];
So, what does diagonally dominant mean? The requirement is that
abs(A(i,i)) >= sum(abs(A(i,j)))
where that sum is taken over all j~=i.
Now, look at the 5th row of a. a(5,5) is 4. The other elements of a on that row are [4 5 5 5]. They sum to 19. Is 19>4? (Yes.) Therefore your matrix is NOT diagonally dominant. Therefore, you would expect the Jacobi iteration on this problem to at least potentially be divergent. In fact, none of the rows of a satisfy the requirement.
Did it diverge? Yes. Case closed. When you see a problem with a method diverging, don't just give up. Ask if it diverged for a valid reason. Do some reading about the method in question.
Suppose we tried a different matrix instead? I'll just use the guts of your code, because I don't feel like finding and saving all of your associated functions.
a = a + eye(5)*30;
c = b./diag(a);
B = (-1./diag(a)')' .*(a-diag(diag(a)));
x0=zeros(length(b),1);
x=c;
nit = 0;
NMAX = 1000;
TOL = 1e-12;
TOLX = 1e-12;
while norm(x-x0,1)>=TOLX && nit<NMAX
x0 = x ;
x = B*x0 + c;
TOLX = max(realmin, TOL*norm(x0,1));
nit = nit+1;
end
Did it work?
x
x =
2.0931
1.7777
0.77996
1.3459
-0.2909
>> a\b
ans =
2.0931
1.7777
0.77996
1.3459
-0.2909
>> a
a =
35 1 1 1 4
5 32 5 3 1
1 3 35 5 5
5 5 3 34 5
4 5 5 5 34
Of course. It worked because the Jacobi method is convergent for that matrix. Note that if you are taking linear algebra, and they are teaching the Jacobi iteration method for linear systems, this should be one of the things they would have dicussed at some point. I suppose it MIGHT be the thing they discuss AFTER you try using it to solve a non-convergent problem.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Model Order Reduction 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
