Jacobi method: Error with index in the while loop.

N = 6;
I = eye(N);
A = toeplitz([2 -1 zeros(1, N-2)]);
b = [0 2 3 -1 2 1]';
x = [0 0 0 0 0 0]'; %initial guess
diagonal = diag(diag(A));
upper = triu(A);
lower = tril(A);
normValue = Inf;

B_jacobi = -inv(diagonal)*(lower+upper);
c_jacobi = inv(diagonal)*b;
eigen_value_jacobi = eig(B_jacobi);

B_gs = -inv(I+(inv(diagonal)*lower))*(inv(diagonal)*upper);
c_gs = inv(I+(inv(diagonal*lower)))*(inv(diagonal)*b);
eigen_value_gs = eig(B_gs);
%disp(A);

%Jacobi method
k = 1;
Tol = 0.001;

while normValue > Tol

      x(k) = B_jacobi*x(k-1) + c_jacboi;
      k= k+1;

      normValue= norm(x(k) - x(k-1));

end

[John BG]

Hello Mr Garcia

this is John BG <mailto:[e-mail] [e-mail]>

The Jacobi elliptic functions are already available with command

ellipj

function help here

https://uk.mathworks.com/help/matlab/ref/ellipj.html?s_tid=srchtitle

Answer 1

[Jan]

This cannot work:

k = 1;
while normValue > Tol
    x(k) = B_jacobi*x(k-1) + c_jacboi;

because in the 1st iteration k=1 you try to evaluate x(k-1), but k-1 is 0.

Simply start with k = 2.

Then replace "c_jacboi" by "c_jacobi".

Now you have the problem, that

B_jacobi*x(k-1) + c_jacobi

replies a matrix, but you try to assign it to the scalar x(k).

Maybe you want:

while normValue > Tol
    xNew      = B_jacobi*x + c_jacboi;
    k         = k+1;
    normValue = norm(x - xNew);
    x         = xNew;
end

Please read the documentation of inv:

doc inv

There you find the important hint, that the slash operator is much better.

B_jacobi = -diagonal \ (lower + upper);