How can I fix this for loop to solve a system of linear equations through the Jacobi iteration method?

Question

I'm trying to write a function that goes through the Jacobi iteration method for solving a system of linear equations. I've got most of it down, I just need to figure out how to iterate the last for loop either 1000 times or until the break condition is met. How can I make it so the value of x updates each iteration?

import numpy as np

def Jacobi(A,b,err):
    n = A.shape
    k = np.zeros(n[0])
    D = np.zeros(n)
    U = np.zeros(n)
    L = np.zeros(n)
    for i in range(n[0]):
        for j in range(n[0]):
            if i == j:
                D[i,j] = A[i,j]
            elif i < j:
                U[i,j] = A[i,j]
            else:
                L[i,j] = A[i,j]
    w = []
    for i in range(1000):
        x = np.linalg.inv(D)*(U+L)*x +np.linalg.inv(D)*b
        w.append(x)

        if abs(w[-1] - w[-2]) < err:
            break
    return w[-1]

For reference, my error statement says a list index in the if clause is out of range. I assume this is because there's only one element in w since I don't know how to make the for loop. Thanks in advance for any help.

Lutz Lehmann · Answer 1 · 2020-03-31T06:40:48.313

I'm pretty sure you missed the intent of that exercise, if you can use inv, then you can also use linalg.inv(A) or better linalg.solve(A,b). Note that you have sign errors and that the multiplication * is not the matrix multiplication between numpy arrays. (Your declaration of the arrays is incompatible with their later use.)

Your specific problem can be solved by adding an additional test

    if i>1 and abs(w[-1] - w[-2]) < err:

when the first condition fails the second is not evaluated.

You should contemplate if it is a waste of memory to construct the w list when all you ever need is the last two entries.

x_last, x = x, jacobi_step(A,b,x)

would also work to have these available.

The preparation can be reduced to

D=np.diag(A); A_reduced = A-np.diag(D);

then the Jacobi step is simply, using that the arithmetic operations are applied element-wise by default

x_last, x = x, (b-A_reduced.dot(x))/D

It would appear that the construction of the list is not needed since he only returns the last element. — sokato, Mar 31 '20 at 06:22
@sokato : You need the last vector to compute the update. This can be done in different ways, but the extra vector can not be avoided. — Lutz Lehmann, Mar 31 '20 at 06:29

How can I fix this for loop to solve a system of linear equations through the Jacobi iteration method?

1 Answers1