Jacobi iterative method in matlab

I just started taking a course in numerical methods and I have an assignment to code the Jacobi iterative method in matlab. So this is my code (and it is working):

function x1 = jacobi2(a,b,x0,tol)
n = length(b);
for j = 1 : n
     x(j) = ((b(j) - a(j,[1:j-1,j+1:n]) * x0([1:j-1,j+1:n])) / a(j,j)); % the first iteration
end
x1 = x';
k = 1;
while norm(x1-x0,1) > tol
    for j = 1 : n
     x_ny(j) = ((b(j) - a(j,[1:j-1,j+1:n]) * x1([1:j-1,j+1:n])) / a(j,j));
    end
    x0 = x1;
    x1 = x_ny';
    k = k + 1;
end
k
x = x1';

I'm assuming there is alot I can do to make this code better since I'm new to matlab, and I would love som feedback on that. But my question is if I instead of what I have done should use the matrix method where we have xk+1 = inv(D) * (b - (L+U) * xk)). Is this a more effective method? And how should I think when deciding what method to use, how do I know what method is more effective?

If someone could help me it would be great!

[lubna ineirat]

what we can do in the fuzzy linear system?

Answer 1

[Bruno Pop-Stefanov]

1. Some feedback about your code

It's good practice to pre-allocate memory before a for loop. This is actually what Code Analyzer suggests for variables x and x_ny. Since you know that x will eventually contain n elements, I would add:

    x = zeros(n,1);

before the first for loop. That way you can also control that x will be a column vector (or a row vector if you use zeros(1,n)) and you do not need to transpose x after the loop.

Same remark for x_ny. Add

    x_ny = zeros(n,1);

before the second for loop.

You might actually be able to vectorize these for loops, if you find a way to rewrite lines 4 and 10.

Here are some general advice for performance:

http://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html

...and about vectorization in particular:

http://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

2. About the matrix method

I am not familiar with the Jacobi method, but I would avoid using inv. Calculating the inverse of a matrix numerically is a risky operation when the matrix is badly conditioned. It's also slower and less precise than other linear solvers. Instead, use mldivide to solve a system of linear equations. Based on how the system looks like, mldivide will choose an appropriate method.

    x(k+1) = D \ (b - (L+U)*x(k));

[Rafid Jabbar]

Dears, Please could one answer me, how I can solve below equation numerically by Jacobi method to get temperature distribution along z-axis, 1D problem, steady state: (

Answer 2

[Antonio Carlos R. Troyman]

It would be intersting to program the Jacobi Method for the generalized form of the eigenvalue problem (the one with separated stiffness and mass matrices). A good reference is the FORTRAN subroutine presented in the book "Numerical Methods in Finite Element Analysis" by Bathe & Wilson, 1976, Prentice-Hall, NJ, pages 458 - 460. If there is someone interested I have this routine in Visual Basic 6, so, please contact me in [e-mail].

Answer 3

[Prajakta pimpalkar]

function x1 = jacobi2(a,b,x0,tol)

n = length(b);

for j = 1 : n

x(j) = ((b(j) - a(j,[1:j-1,j+1:n]) * x0([1:j-1,j+1:n])) / a(j,j)); % the first iteration

end

x1 = x';

k = 1;

while norm(x1-x0,1) > tol

for j = 1 : n

x_ny(j) = ((b(j) - a(j,[1:j-1,j+1:n]) * x1([1:j-1,j+1:n])) / a(j,j));

end

x0 = x1;

x1 = x_ny';

k = k + 1;

end

k

x = x1';

[Okiki Akinsooto]

When used, its displaying error using - in line 15

While norm(x1-x0,1)>tol

Answer 4

[My]

function x = myjacobi(A,b,x0,N,TOL)

x = x0;

for k = 1:N
    x_new = x;
    
    for i = 1:length(x)
        x_new(i) = (b(i) - A(i,[1:i-1 i+1:end])*x([1:i-1 i+1:end])) / A(i,i);
    end
    
    if norm(x_new - x) < TOL
        x = x_new;
        return
    end
    
    x = x_new;
end

end

%Main

clc; clear;

A = [4 2 1;
    -1 2 0;
     2 1 4];

b  = [11; 3; 16];
x0 = [1; 1; 1];

x = myjacobi(A,b,x0,100,1e-3)

x = 3×1

1.0001 1.9999 3.0001

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Answer 5

[David Goodmanson]

the expression

x_new = (b - A*x + diag(A).*x)./diag(A)

works. I don't know if it might run into numerical accuracy issues due to adding and subtracting the same value twice, but I ran the example with very small tolerance and there was no difference between this and the for loop version.