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
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
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
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 个评论
Antonio Carlos R. Troyman
2022-4-18
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.
0 个评论
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!