Jacobi iterative method in matlab

113 次查看（过去 30 天）
ribagorza
ribagorza 2014-10-7
编辑： David Goodmanson 2025-12-22，11:31
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!

请先登录，再回答此问题。

回答（5 个）

Bruno Pop-Stefanov
Bruno Pop-Stefanov 2014-10-8
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));
  1 个评论
Rafid Jabbar
Rafid Jabbar 2017-5-15
编辑：Rafid Jabbar 2017-5-15

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: (

请先登录，再进行评论。

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 troy@oceanica.ufrj.br.
Prajakta pimpalkar
Prajakta pimpalkar 2021-9-28
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';
  1 个评论
Okiki Akinsooto
Okiki Akinsooto 2023-6-10
When used, its displaying error using - in line 15
While norm(x1-x0,1)>tol

请先登录，再进行评论。

My
My 2025-12-21，20:03
编辑：Torsten 2025-12-21，21:54
Ran in:
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>
David Goodmanson
David Goodmanson 2025-12-22，3:56
编辑：David Goodmanson 2025-12-22，11:31
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.

类别

在 Help Center 和 File Exchange 中查找有关 Startup and Shutdown 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by