Jacobi Iteration With Ruby
April 28 2010
When solving linear algebraic systems, it is often beneficial to use iterative methods. Iterative solvers are particularly useful when a matrix is diagonally dominant. One such iterative method that is particularly useful is Jacobi Iteration.
Given a matrix A, you first decompose it into 3 separate matrices: L, D, and U. These are the lower, diagonal, and upper sections of the matrix, respectively.
For example, let’s assume we’ve been given this matrix A:
We will then zero-out everything above and including the diagonal entries. This will give us our new L matrix.
Likewise, we’ll remove all entries below, and including the diagonal from matrix A to give us our matrix U.
Finally, we’ll zero-out all non-diagonal entries, to yield our new matrix D.
As stated previously, since A=L+D+U, we can deduce that Ax=Dx+(L+U)x. And since Ax=b, we can conclude that Dx=b-(L+U)x, which is the basis of Jacobi’s Method.
Given an estimate vector (x) (which is generally initialized to 0 if no educated guess can be provided), the next estimate is found by:
And since D is a diagonal matrix, the above equation can be rewritten in terms of the components of A and b.
So now that we understand how Jacobi Iteration works, we can write a program to do it for us! Ruby is not necessarily the best language to do this in, as these systems can often times be quite large. Compiled languages such as C++ are going to converge much faster.
To start out, we’re going to write a method that takes in our Matrix A, our right hand size b, the number of times we’d like to iterate, and our initial guess.
From here, we will need to fire up our estimate vector x. We’ll initialize its first element to be the initial “guess,” and the remaining entries to be 0.
Finally, we implement the equation above as nested loops to iteratively try to converge to a solution.
And there you have it! If you're doing this on large matrices, I would definitely suggest using a different language, such as C++ to handle this, since Ruby can not only be a bit slower, but it's also quite a memory hog when dealing with so many objects.