Solving Ax=b with matrices (bigger than 2 by 2 matrices)

I have been working on creating a program that solves linear systems of equations for the Jacobi and Gauss-Seidel iterative methods. (Link to the methods: https://www.cis.upenn.edu/~cis515/cis515-12-sl5.pdf) The specific problem I have with my code is that it works for 2 by 2 matrices but does not work for anything larger than that. My sample 2x2 matrix is A=[[3,1],[2,-1]] and b=[5,0]. This works just fine but the 3x3 matrix I am trying to solve is A=[[1,1,1],[0,2,5],[2,5,-1]] and b=[6,-4,27] and it does not compute the correct answer.

import numpy as np
from math import *
from matplotlib import pyplot as plt
x_list = []
x_list_1 = []

it_1 = int(input("What is the maximum amount of iterations: "))
it= it_1+1
n = int(input("Enter the matrix size: "))

A= np.zeros((n,n))

for i in range(n):
A[i] = input().split(" ")

MatrixA = np.array(A)
print("A= ")
print (A)

for i in range(A.shape[0]):
row = ["{}*x{}".format(A[i, j], j + 1) for j in range(A.shape[1])]

b = np.zeros((n,1))

for i in range(n):
b[i] = input().split(" ")

MatrixB = np.array(b)
print("b= ")
print (b)

for i in range(b.shape[0]):
row = ["{}*x{}".format(A[i, j], j + 1) for j in range(b.shape[1])]

x = np.zeros_like(b)
for it_count in range(it):
print("Iteration", it_count, ":")
print(x)
x_new = np.zeros_like(x)

for i in range(A.shape[0]):
norm = np.linalg.norm(x)
x_list.append(norm)
s1 = np.dot(A[i, :i], x[:i])
s2 = np.dot(A[i, i + 1:], x[i + 1:])
x_new[i] = (b[i] - s1 - s2) / A[i, i]

if np.allclose(x, x_new, atol=1e-25, rtol=0):
break

x = x_new

print("Solution (Jacobi):")
print(x)

x = np.zeros_like(b)
def gauss(A, b, x, it):

counting = range(0,it)
lower_triangle = np.tril(A)
upper_triangle = A - lower_triangle
for it_count in range(it):
print("Iteration", it_count, ":")
print(x)
norm = np.linalg.norm(x)
x_list_1.append(norm)
x = np.dot(np.linalg.inv(lower_triangle), b - np.dot(upper_triangle, x))
return x

print("")
z = gauss(A, b, x, it)
print("Solution(Gauss-Seidel):")
print(z)

Jacobi_Magnitude_List = list(dict.fromkeys(x_list))
plt.plot(Jacobi_Magnitude_List, label="Jacobi Magnitude")
plt.plot(x_list_1, label = "Gauss-Seidel Magnitude")
plt.title("Jacobi and Gauss-Seidel Method Magnitude at each Iteration")
plt.xlabel("Iteration Number")
plt.ylabel("Magnitude")
plt.legend()
plt.show()


I am not sure what in my code is limiting it to only solving 2 by 2 matrices.

• Not to detract from the post, but aren’t there built-in commands that can do this for you? – Mnifldz Nov 19 '19 at 2:22
• Tangential comment, but when you're in the debugging phase you shouldn't read in the matrix from user input like this. Just define $A$ in the code, with no user input. It will make debugging much faster. – littleO Nov 19 '19 at 2:27