1
$\begingroup$

I implemented the Jacobi iteration using Matlab based on this paper, and the code is as follows:

function x = jacobi(A, b)
% Executes iterations of Jacobi's method to solve Ax = b.

assert(size(A, 1) == size(A, 2))
assert(size(A, 1) == size(b, 1))

% matrix with all zeros except the diagonal elements,
% which are those of A.
D = diag(diag(A));

% matrix A but with its diagonal elements equal zero
R = D - A;

tol = 1.e-8; 
max_num_of_iter = 1000;

% initialize relative error to large value
relative_error = inf; 
num_of_iter = 1;

% initial guess
x_prev = zeros(size(b, 1), 1);

% iterate until convergence or maximum number of iterations is reached
while relative_error > tol && num_of_iter < max_num_of_iter
    x = D \ (b + R * x_prev);

    relative_error = norm(x - x_prev, inf) / norm(x, inf);

    x_prev = x;

    num_of_iter = num_of_iter + 1;
end

Now what I don't understand is why the iteration x is found using D \ (b + R * x_prev), more specifically why do we have a + sign instead of a - sign in D \ (b + R * x_prev)?

According to [Wikipedia](https://en.wikipedia.org/wiki/Jacobi_method, the Jacobi method is as follows:

$$\mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -R\mathbf {x} ^{(k)})$$

where $D$ is an $n \times n$ matrix whose entries are all zeros except the diagonal entries which are the diagonal entries of $A$. $R = D - A$, i.e. a matrix whose diagonal elements are zeros and all other are the same as the corresponding elements in $A$.

Clearly, if we want to solve

$$Ax = b$$

and $$A = L + D + U =\\ D + L + U = \\ D + R \Rightarrow \\ A - D = R = \\ L + U $$
then

$$(D + L + U) x = b \\ Dx + Lx + Ux = b \Rightarrow\\ Dx = b - Lx - Ux \Rightarrow \\ Dx = b - (L + U) x$$

since we have defined $L + U = R$, then

$$Dx = b - (L + U) x \Rightarrow \\ Dx = b - Rx \Rightarrow \\ x = D^{-1} \times (b - Rx)$$

So the original formula of the iterative method $\mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -R\mathbf {x} ^{(k)})$ provided by Wikipedia seems to be consistent with my calculations. The problem is that if I use a - instead of + I don't find the right solutions. How do I know that? I'm using x = A \ b (i.e. guassian elimination) to check for the solution.

$\endgroup$
4
$\begingroup$

In your case (line 12 in your code) $$A = D-R.$$ Then $$ \begin{aligned} A\mathbf{x} {}&{}= D\mathbf{x}-R\mathbf{x}\\ {}&{}= \mathbf{b}\,, \end{aligned}$$ which rearranges to $$\mathbf{x} = D^{-1}\left( \mathbf{b} + R\mathbf{x} \right)\,,$$ which leads to the iteration which is the the first line in your while loop.

This only differ's from the Wikipedia version by the definition, i.e. Wikipedia's $$R' = A - D$$ to your $$R = D - A\,.$$

$\endgroup$
  • $\begingroup$ So why would the Wikipedia be wrong, since it doesn't give me the right results? $\endgroup$ – nbro Oct 27 '16 at 12:07
  • 1
    $\begingroup$ It's not wrong, it's just the definition is different. If you change your line 12 R = D-A to the wikipedia definition R = A-D, then the wikipedia iteration x = D \ (b - R * x_prev) should be correct. $\endgroup$ – Steve Oct 27 '16 at 12:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.