Author: David Merkley
Language: Python
Links to Code:
Task 1
Description/Purpose: This is the code that is used for all of the vector operations and I added it to my software manual.
Implementation/Code:
import random
def printMatrix(A, ind=0):
n = len(A[0])
for row in A:
printVector(row, ind)
def printVector(row, ind=0):
print(" "*ind, end="")
for x in row:
print("{:10.5f}".format(x), end="")
print()
def rv(n):
return [random.random() for i in range(n)]
def vectorAdd(x, y):
return [x[i] + y[i] for i in range(len(x))]
def vectorSub(x, y):
return [x[i] - y[i] for i in range(len(x))]
def vectorScalarMultiply(x, a):
return [a * xi for xi in x]
def vectorDotProduct(x, y):
return sum([x[i] * y[i] for i in range(len(x))])
def vectorOuterProduct(x, y):
return [[xi * yj for yj in y] for xi in x]
x = [1, 2, 3]
y = [1, 1, 1]
printVector(x)
printVector(y)
printVector(vectorAdd(x, y))
printVector(vectorSub(x, y))
printVector(vectorScalarMultiply(x, 3))
printVector([vectorDotProduct(x, y)])
printMatrix(vectorOuterProduct(x, y))
Output:
1.00000 2.00000 3.00000
1.00000 1.00000 1.00000
2.00000 3.00000 4.00000
0.00000 1.00000 2.00000
3.00000 6.00000 9.00000
6.00000
1.00000 1.00000 1.00000
2.00000 2.00000 2.00000
3.00000 3.00000 3.00000
Task 2
Description/Purpose: This is the code used for all of the norm operations. I added this to my software manual as well.
Implementation/Code:
from math import sqrt
from math import isinf
def vectorSub(x, y):
return [x[i] + y[i] for i in range(len(x))]
def l1(x):
return sum([abs(xi) for xi in x])
def l22(x):
return sum([abs(xi)*abs(xi) for xi in x])
def l2(x):
return sqrt(sum([abs(xi)*abs(xi) for xi in x]))
def linfty(x):
return max([abs(xi) for xi in x])
def lp(x, p):
if p == 1:
return l1(x)
elif p == 2:
return l2(x)
elif isinf(p):
return linfty(x)
else:
return sum([abs(xi)**p for xi in x])**(1/p)
def lperr(x, y, p):
return lp(vectorSub(x, y), p)
def l1err(x, y):
return lperr(x, y, 1)
def l22err(x, y):
return l22(vectorSub(x, y))
def l2err(x, y):
return lperr(x, y, 2)
def linftyerr(x, y):
return lperr(x, y, float('inf'))
def printVector(row, ind=0):
print(" "*ind, end="")
for x in row:
print("{:10.5f}".format(x), end="")
print()
x = [1, 2, 3]
y = [1, 1, 1]
printVector(x)
printVector(y)
print(l1(x))
print(l22(x))
print(l2(x))
print(linfty(x))
print(lp(x, 5))
print(lperr(x, y, 5))
print(l1err(x, y))
print(l22err(x, y))
print(l2err(x, y))
print(linftyerr(x, y))
Output:
1.00000 2.00000 3.00000
1.00000 1.00000 1.00000
6
14
3.7416573867739413
3
3.077384885394063
4.194902104167999
9
29
5.385164807134504
4
Task 3
Description/Purpose: Verifies quadratic convergence using Newton’s method.
Implementation/Code:
def vectorAdd(x, y):
return [x[i] + y[i] for i in range(len(x))]
def vectorSub(x, y):
return [x[i] + y[i] for i in range(len(x))]
def vectorScalarMultiply(x, a):
return [a * xi for xi in x]
def printMatrix(A, ind=0):
n = len(A[0])
for row in A:
printVector(row, ind)
def printVector(row, ind=0):
print(" "*ind, end="")
for x in row:
print("{:10.5f}".format(x), end="")
print()
def matrixAddition(A, B):
return [vectorAdd(A[i], B[i]) for i in range(len(A))]
def matrixSubtraction(A, B):
return [vectorSub(A[i], B[i]) for i in range(len(A))]
def matrixScalarMultiply(A, a):
return [vectorScalarMultiply(A[i], a) for i in range(len(A))]
def transpose(A, doMutate=False):
if doMutate:
A = [[A[j][i] for j in range(len(A))] for i in range(len(A))]
return A
else:
return [[A[j][i] for j in range(len(A))] for i in range(len(A))]
def matrixVectorMultiply(A, b):
m = len(A)
n = len(A[0])
ans = [0 ] *m
for i in range(m):
for j in range(n):
ans[i] += A[i][j] * b[j]
return ans
def matrixMatrixMultiply(A, B):
m = len(A)
n = len(A[0])
r = len(B[0])
C = [[0 ] *r for i in range(m)]
for i in range(m):
for j in range(n):
for k in range(r):
C[i][k] += A[i][j] * B[j][k]
return C
A = [[1, 2], [3, 4]]
B = [[1, 1], [1, 1]]
C = [[1, 2], [3, 4], [5, 6]]
c = [1, 2]
print("Matrix Scalar Multiply")
printMatrix(matrixScalarMultiply(A, 2))
print("Matrix Add")
printMatrix(matrixAddition(A, B))
print("Matrix Sub")
printMatrix(matrixSubtraction(A, B))
print("Matrix Transpose")
printMatrix(transpose(A))
print("Matrix Vector Multiply")
printVector(matrixVectorMultiply(C, c))
print("Matrix Matrix Multiply")
printMatrix(matrixMatrixMultiply(C, A))\
Output:
Matrix Scalar Multiply
2.00000 4.00000
6.00000 8.00000
Matrix Add
2.00000 3.00000
4.00000 5.00000
Matrix Sub
2.00000 3.00000
4.00000 5.00000
Matrix Transpose
1.00000 3.00000
2.00000 4.00000
Matrix Vector Multiply
5.00000 11.00000 17.00000
Matrix Matrix Multiply
7.00000 10.00000
15.00000 22.00000
23.00000 34.00000
Task 4
Description/Purpose: This is the code used for the Jacobi Iteration. Added to my software manual.
Implementation/Code:
def l22err(x, y):
return l22(vectorSub(x, y))
def vectorSub(x, y):
return [x[i] - y[i] for i in range(len(x))]
def l22(x):
return sum([abs(xi)*abs(xi) for xi in x])
def residual(A, b, x):
n = len(b)
r = [1.0 * ri for ri in b]
for i in range(n):
for j in range(n):
r[i] -= A[i][j] * x[j]
return r
def JacobiIteration(A, b, x0, tol=1e-10, maxiter=1e3):
n = len(b)
tol2 = tol * tol
err2 = 2*tol
it = 0
while err2 > tol2 and it < maxiter:
r = residual(A, b, x0)
x1 = [x0[i] + r[i] / A[i][i] for i in range(n)]
err2 = l22err(x1, x0)
x0 = x1
it += 1
if it == maxiter:
print("Max Iter")
return x1
Task 5
Description/Purpose: This is code to generate a random diagonally dominant matrix. Then errors for Jacobi and Gaussian elimination are below.
Implementation/Code:
def generateSquareRandomSymmetricDiagonallyDominant(n):
A = [[0]*n for i in range(n)]
for i in range(n):
A[i][i] = n
for j in range(i + 1, n):
A[i][j] = random.random()
for j in range(i):
A[i][j] = A[j][i]
return A
A = generateSquareRandomSymmetricDiagonallyDominant(100)
x = [1]*100
b = matrixVectorMultiply(A, x)
x0 = rv(100)
xj = JacobiIteration(A, b, x0, 1e-13, 1e3)
xe = rowEchelonFormAndBackSub(A, b, doMutate=True, leaveUpperTriangularImplicit=True, scalePartialPivoting=False)
print("Jacobi Error:", lperr(xj, x, 4))
print("Euclid Error:", lperr(xe, x, 4))
Output:
Jacobi Error: 8.430620745549493e-15
Euclid Error: 3.215571698321486e-15
Task 6
In the wikipedia for the Jacobi method it says, “the Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations.” So Jacobi method is for a matrix that we know is diagonally dominant. While the Gauss-Seidel method is more general and can solve the system of equations without it needing to be diagonally dominant. This is super interesting to think about in terms of length of time to compute each of these.