meyers hw 7.ipynb

Kernel: SageMath 8.2

def height(m):
  return len(m)

def gauss_seidel(m, x0=None, eps=1e-5, max_iteration=100):
  """
  Parameters
  ----------
  m  : list of list of floats : coefficient matrix
  x0 : list of floats : initial guess
  eps: float : error tolerance
  max_iteration: int
  
  Returns
  -------  
  list of floats
      solution to the system of linear equation
  
  Raises
  ------
  ValueError
      Solution does not converge
  """
  n  = height(m)
  x0 = [0] * n if x0 == None else x0
  x1 = x0[:]

  for __ in range(max_iteration):
    for i in range(n):
      s = sum(-m[i][j] * x1[j] for j in range(n) if i != j) 
      x1[i] = (m[i][n] + s) / m[i][i]
    if all(abs(x1[i]-x0[i]) < eps for i in range(n)):
      return x1 
    x0 = x1[:]    
  raise ValueError('Solution does not converge')

if __name__ == '__main__':
  m = [[3,-1,5], [-1,2,4]]
  print(gauss_seidel(m))

[14108771/5038848, 34264163/10077696]

def jacobi(A,b,N=25,x=None):
    """Solves the equation Ax=b via the Jacobi iterative method."""
    # Create an initial guess if needed                                                                                                                                                            
    if x is None:
        x = zeros(len(A[0]))

    # Create a vector of the diagonal elements of A                                                                                                                                                
    # and subtract them from A                                                                                                                                                                     
    D = np.diag(A)
    R = A - np.diagflat(D)

    # Iterate for N times                                                                                                                                                                          
    for i in range(N):
        x = (b - np.dot(R,x)) / D
    return x

A = np.array([[2.0,1.0],[5.0,7.0]])
b = np.array([11.0,13.0])
guess = np.array([1.0,1.0])

sol = jacobi(A,b,N=25,x=guess)

print "A:"
print(A)

print "b:"
print(b)

print "x:"
print(sol)

A:
[[ 2.  1.]
 [ 5.  7.]]
b:
[ 11.  13.]
x:
[ 7.11110202 -3.22220342]

def generate(condition_number):
    import numpy as np
    i = 0
    j = 0
    A = np.zeros((condition_number, condition_number))
    while i< condition_number:
        while (j <condition_number):
            if i == j:
                A[i][j]=3
            elif np.abs(i-j)==1:
                A[i][j]=-1
            else:
                j = j + 1
        print(A[i])
        i = i + 1
        j = 0
    b = np.ones((condition_number,1))
    
    
    #return (b)
generate(6)        

#Had trouble generating the matrix, so I gave up since I can clearly computer it above