Jacobi iterative method in matlab

765 次查看(过去 30 天)
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!

回答(3 个)

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:
...and about vectorization in particular:
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: (

请先登录,再进行评论。


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

请先登录,再进行评论。


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.

类别

Help CenterFile Exchange 中查找有关 MATLAB 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by