Write a program in Python to solve a linear system of the form Ax = b by Gaussia

Question

Write a program in Python to solve a linear system of the form Ax = b by Gaussian elimination with scaled partial pivoting. You should pass the matrix A and the right hand side vector b, and the value of n (the number of rows and columns in the matrix). You will need to have vectors for the scaling factors and the index of the rows internal to the program. You can assume A is a square matrix. Your calling statements should be: def gauss (A, b, n) Your output statement should be of the form: The solution vector is XXXXXX. Test your program on the following: The Larger Numerical Example on page 87 of the text. The system in exercise #13b on page 99. The system Ax^rightarrow = b^rightarrow, where A = [a_i, j] and a_i, j = 1/1 + j + 1 for I, j = 0 ... n - 1 and b[i] = Sigma^n - 1_ j = 0 1/i + j + 1 (which is the sum of the entries in the row of A (so the solution is all ones). The matrix A is the Hilbert matrix and the way the b vector is constructed makes the solution obvious. Run your program for n = 6 and n = 15, which are 7x7 and 16x16 matrices. Compute the 2-norms of the residual vector, r^OverBar, and the error vector e^OverBar = x - x^OverBar. Use Maple to compute the condition number of both matrices and report those numbers in your Python window as a comment. Are the solutions accurate or not? Could we predict this from the condition numbers? Discuss your findings in complete sentences.

Explanation / Answer

Answer: See the code below:

------------------------------------------------

import numpy as np
def gaussy(A,b,n):
    #forward elimination phase of procedure
    l=[] #index vector
    s=[] #scaling vector
    r=0.0 #residual
    smax=0.0
    i=0
    j=0
    k=0
    xmult=0.0
    for i in range(n):
        l[i]=i
        smax=0.0
        for j in range(n):
            det_A=np.linalg.det(A)
            smax=max(smax,det_A)
        s[i]=smax
    rmax=0
    for k in range(n-1):
        rmax=0
        for i in range(k,n):
            r=abs(A[l[i]][k]/s[l[i]])
            if r > rmax:
                rmax=r
                j=i
        l[j]=l[k]
        for i in range(k+1,n):
            xmult=A[l[i]][k]/A[l[k]][k]
            A[l[i]][k]=xmult
            for j in range(k+1,n):
                A[l[i]][j]=A[l[i]][j]-xmult*A[l[k]][j]

    #back substitution phase of procedure
    sum=0.0
    x=[] #root vector
    for k in range(n-1):
        for i in range(k+1,n):
            b[l[i]]=b[l[i]]-A[l[i]][k]*b[l[k]]
    x[n]=b[l[n]]/A[l[n]][n]
    for i in range(n-1,1):
        sum=b[l[i]]
        for j in range(i+1,n):
            sum=sum-A[l[i]][j]*x[j]
        x[i]=sum/A[l[i]][i]
    return x  
#Read n
n=input("Enter value of n:")
n=int(n)

#Create matrix
A = [[0 for i in range(0,n+1)] for j in range(0,n+1)]
b=[0 for i in range(0,n+1)]

print("Creating Hilbert matrix A...")
for i in range(n+1):
    for j in range(n+1):
        A[i][j]=1.0/float(i+j+1)
      

print("Hilbert matrix A:")
for i in range(n+1):
    for j in range(n+1):
        print(str(A[i][j])+" ")
    print(" ")

print("Creating matrix b...")
for i in range(n+1):  
    for j in range(n+1):
        b[i]=b[i]+A[i][j]  

print("Matrix b:")
for i in range(n+1):
    print(str(b[i]))

x=gaussy(A,b,n)
disp(x)

-----------------------------------------------

Related Questions

Navigate

• Browse (All)
• Browse W
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.