626.sagews

def jacobi_find_root(m, b, x0, num_iterations=10, full_output=True):
    """
    Use the Jacobi method to approximate a solution to the system of
    equations represented by m*x = b with initial guess x=x0
    """  
    
    #create diagonal matrix Q from input matrix m
    #using only the values of m on the main diagonal
    Q = diagonal_matrix([m[i,i] for i in range(m.ncols())])
    #calculate these outside of loop for efficiency purposes
    Qinv = Q.inverse()
    QinvTimesB = Qinv*b
    IMinusQinvTimesM = identity_matrix(RR,2) - Qinv*m1
    for i in range(num_iterations):
        x0 = (IMinusQinvTimesM)*x0 + QinvTimesB
        if full_output == True: 
            print x0

    return x0

m1 = matrix([[3,1],[1,2]])
b = vector([5,5])
x0 = vector([0.0, 0.0])
jacobi_find_root(m1, b,x0)

(1.66666666666667, 2.50000000000000) (0.833333333333333, 1.66666666666667) (1.11111111111111, 2.08333333333333) (0.972222222222222, 1.94444444444444) (1.01851851851852, 2.01388888888889) (0.995370370370370, 1.99074074074074) (1.00308641975309, 2.00231481481481) (0.999228395061728, 1.99845679012346) (1.00051440329218, 2.00038580246914) (0.999871399176955, 1.99974279835391) (0.999871399176955, 1.99974279835391)