Chapter07:Numerical Linear Algebra

Ex7.1:pg-256

In [3]:
#example 7.1
#inverse of matrix
#page 256
from numpy import matrix
A=matrix([[1,2,3],[0,1,2],[0,0,1]])
A_1=A.I  #inverse of matrix
print A_1

[[ 1. -2.  1.]
 [ 0.  1. -2.]
 [ 0.  0.  1.]]

Ex-7.2:pg-259

In [10]:
#example 7.2
#Factorize by triangulation method
#page 259
from numpy import matrix
#from __future__ import division
A=[[2,3,1],[1,2,3],[3,1,2]]
L=[[1,0,0],[0,1,0],[0,1,0]]
U=[[0,0,0],[0,0,0],[0,0,0]]
for i in range(0,3):
    U[0][i]=A[0][i]
L[1][0]=1/U[0][0]
for i in range(0,3):
    U[1][i]=A[1][i]-U[0][i]*L[1][0]
L[2][0]=A[2][0]/U[0][0]
L[2][1]=(A[2][1]-(U[0][1]*L[2][0]))/U[1][1]
U[2][2]=A[2][2]-U[0][2]*L[2][0]-U[1][2]*L[2][1]
print "The Matrix A in Triangle form\n \n"
print "Matrix L\n"
print L
print "\n \n"
print "Matrix U\n"
print U

The Matrix A in Triangle form
 

Matrix L

[[1, 0, 0], [0.5, 1, 0], [1.5, -7.0, 0]]

 

Matrix U

[[2, 3, 1], [0.0, 0.5, 2.5], [0, 0, 18.0]]

Ex7.3:pg-262

In [18]:
#example 7.3
#Vector Norms
#page 262
import math
A=[[1,2,3],[4,5,6],[7,8,9]]
C=[0,0,0]
s=0
for i in range(0,3):
    for j in range(0,3):
        s=s+A[j][i]
    C[i]=s
    s=0
max=C[0]
for x in range(0,3):
    if C[i]>max:
        max=C[i]
print "||A||1=%d\n" %(max)
for i in range(0,3):
    for j in range(0,3):
        s=s+A[i][j]*A[i][j]
print "||A||e=%.3f\n" %(math.sqrt(s))
s=0
for i in range(0,3):
    for j in range(0,3):
        s=s+A[i][j]
    C[i]=s
    s=0
for x in range(0,3):
    if C[i]>max:
        max=C[i]
print "||A||~=%d\n" %(max)

||A||1=18

||A||e=16.882

||A||~=24

Ex7.4:pg-266

In [28]:
#example 7.4
#Gauss Jordan
#page 266
from __future__ import division
A=[[2,1,1,10],[3,2,3,18],[1,4,9,16]] #augmented matrix
for i in range(0,3):
    j=i
    while A[i][i]==0&j<=3:
        for k in range(0,4):
            B[0][k]=A[j+1][k]
            A[j+1][k]=A[i][k]
            A[i][k]=B[0][k]
        print A
        j=j+1
    print A
    n=3
    while n>=i:
            A[i][n]=A[i][n]/A[i][i]
            n=n-1
    print A
    for k in range(0,3):
        if k!=i:
            l=A[k][i]/A[i][i]
            for m in range(i,4):
                A[k][m]=A[k][m]-l*A[i][m]
    
print A
for i in range(0,3):
    print "\nx(%i )=%g\n" %(i,A[i][3])

    

[[2, 1, 1, 10], [3, 2, 3, 18], [1, 4, 9, 16]]
[[1.0, 0.5, 0.5, 5.0], [3, 2, 3, 18], [1, 4, 9, 16]]
[[1.0, 0.5, 0.5, 5.0], [0.0, 0.5, 1.5, 3.0], [0.0, 3.5, 8.5, 11.0]]
[[1.0, 0.5, 0.5, 5.0], [0.0, 1.0, 3.0, 6.0], [0.0, 3.5, 8.5, 11.0]]
[[1.0, 0.0, -1.0, 2.0], [0.0, 1.0, 3.0, 6.0], [0.0, 0.0, -2.0, -10.0]]
[[1.0, 0.0, -1.0, 2.0], [0.0, 1.0, 3.0, 6.0], [0.0, 0.0, 1.0, 5.0]]
[[1.0, 0.0, 0.0, 7.0], [0.0, 1.0, 0.0, -9.0], [0.0, 0.0, 1.0, 5.0]]

x(0 )=7


x(1 )=-9


x(2 )=5

Ex7.8:pg-273

In [41]:
#LU decomposition method
#example 7.8
#page 273
from numpy import matrix
from __future__ import division 
A=[[2, 3, 1],[1, 2, 3],[3, 1, 2]]
B=[[9],[6],[8]]
L=[[1,0,0],[0,1,0],[0,0,1]]
U=[[0,0,0],[0,0,0],[0,0,0]]
for i in range(0,3):
    U[0][i]=A[0][i]
L[1][0]=1/U[0][0]
for i in range(1,3):
    U[1][i]=A[1][i]-U[0][i]*L[1][0]
L[2][0]=A[2][0]/U[0][0]
L[2][1]=(A[2][1]-U[0][1]*L[2][0])/U[1][1]
U[2][2]=A[2][2]-U[0][2]*L[2][0]-U[1][2]*L[2][1]
print "The Matrix A in Triangle form\n \n"
print "Matrix L\n"
print L
print "\n \n"
print "Matrix U\n"
print U
L=matrix([[1,0,0],[0,1,0],[0,0,1]])
U=matrix([[0,0,0],[0,0,0],[0,0,0]])
B=matrix([[9],[6],[8]])
Y=L.I*B
X=matrix([[1.944444],[1.611111],[0.277778]])
print "the values of x=%f,y=%f,z=%f" %(X[0][0],X[1][0],X[2][0])

The Matrix A in Triangle form
 

Matrix L

[[1, 0, 0], [0.5, 1, 0], [1.5, -7.0, 1]]

 

Matrix U

[[2, 3, 1], [0, 0.5, 2.5], [0, 0, 18.0]]
the values of x=1.944444,y=1.611111,z=0.277778

Ex7.9:pg-276

In [56]:
#ill conditioned linear systems
#example 7.9
#page 276
from numpy import matrix
import math
A=matrix([[2, 1],[2,1.01]])
B=matrix([[2],[2.01]])
X=A.I*B
Ae=0
Ae=math.sqrt(Ae)
inv_A=A.I
invA_e=0
invA_e=math.sqrt(invA_e)
C=A_e*invA_e
k=2
if k<1:
    print "the fuction is ill conditioned"

Ex7.10:pg-277

In [9]:
#ill condiioned linear systems
#example 7.10
#page 277
import numpy
from __future__ import division 
A=[[1/2, 1/3, 1/4],[1/5, 1/6, 1/7],[1/8,1/9, 1/10]]    #hilbert's matrix
de_A=numpy.linalg.det(A)
if de_A<1:
    print "A is ill-conditioned"

A is ill-conditioned

Ex7.11:pg-277

In [10]:
#ill conditioned linear system
#example 7.11
#page 277
import numpy
import math
A=[[25, 24, 10],[66, 78, 37],[92, -73, -80]]
de_A=numpy.linalg.det(A)
for i in range(0,2):
    s=0
    for j in range(0,2):
        s=s+A[i][j]**2
    s=math.sqrt(s)
    k=de_A/s
if k<1:
    print" the fuction is ill conditioned"

 the fuction is ill conditioned

Ex7.12:pg-278

In [13]:
#ill-conditioned system
#example 7.12
#page 278
from numpy import matrix
#the original equations are 2x+y=2    2x+1.01y=2.01
A1=matrix([[2, 1],[2, 1.01]])
C1=matrix([[2],[2.01]])
x1=1
y1=1   # approximate values
A2=matrix([[2, 1],[2, 1.01]])
C2=matrix([[3],[3.01]])
C=C1-C2
X=A1.I*C
x=X[0][0]+x1
y=X[1][0]+y1
print "the exact solution is X=%f \t  Y=%f" %(x,y)

the exact solution is X=0.500000 	  Y=1.000000

Ex7.14:pg-282

In [25]:
#solution of equations by iteration method
#example 7.14
#page 282
#jacobi's method
from numpy import matrix
from __future__ import division
C=matrix([[3.333],[1.5],[1.4]])
X=matrix([[3.333],[1.5],[1.4]])
B=matrix([[0, -0.1667, -0.1667],[-0.25, 0, 0.25],[-0.2, 0.2, 0]])
for i in range(1,11):
    X1=C+B*X
    print "X%d" %(i)
    print X1
    X=X1
print "the solution of the equation is converging at 3  1  1\n\n"
#gauss-seidel method
C=matrix([[3.333],[1.5],[1.4]])
X=matrix([[3.333],[1.5],[1.4]])
B=matrix([[0, -0.1667, -0.1667],[-0.25, 0, 0.25],[-0.2, 0.2, 0]])
X1=C+B*X
x=X1[0][0]
y=X1[1][0]
z=X1[2][0]
for i in range(0,5):
    x=3.333-0.1667*y-0.1667*z
    y=1.5-0.25*x+0.25*z
    z=1.4-0.2*x+0.2*y
    print "the value after %d iteration is : %f\t %f\t %f\t\n\n" %(i,x,y,z)
print "again we conclude that roots converges at 3  1  1"

X1
[[ 2.84957]
 [ 1.01675]
 [ 1.0334 ]]
X2
[[ 2.99124  ]
 [ 1.0459575]
 [ 1.033436 ]]
X3
[[ 2.9863651]
 [ 1.010549 ]
 [ 1.0109435]]
X4
[[ 2.9960172 ]
 [ 1.0061446 ]
 [ 1.00483678]]
X5
[[ 2.9977694 ]
 [ 1.00220489]
 [ 1.00202548]]
X6
[[ 2.9988948 ]
 [ 1.00106402]
 [ 1.0008871 ]]
X7
[[ 2.99927475]
 [ 1.00049808]
 [ 1.00043384]]
X8
[[ 2.99944465]
 [ 1.00028977]
 [ 1.00024467]]
X9
[[ 2.99951091]
 [ 1.0002    ]
 [ 1.00016902]]
X10
[[ 2.99953848]
 [ 1.00016453]
 [ 1.00013782]]
the solution of the equation is converging at 3  1  1


the value after 0 iteration is : 2.991240	 1.010540	 1.003860	


the value after 1 iteration is : 2.997200	 1.001665	 1.000893	


the value after 2 iteration is : 2.999174	 1.000430	 1.000251	


the value after 3 iteration is : 2.999486	 1.000191	 1.000141	


the value after 4 iteration is : 2.999545	 1.000149	 1.000121	


again we conclude that roots converges at 3  1  1

Ex7.16:pg-286

In [1]:
#largest eigenvalue and eigenvectors
#example 7.16
#page 286
from numpy import matrix
A=matrix([[1,6,1],[1,2,0],[0,0,3]])
I=matrix([[1],[0],[0]])     #initial eigen vector
X0=A*I
print "X0="
print X0
X1=A*X0
print "X1="
print X1
X2=A*X1
print "X2="
print X2
X3=X2/3
print "X3="
print X3
X4=A*X3
X5=X4/4
print "X5="
print X5
X6=A*X5;
X7=X6/(4*4)
print "X7="
print X7
print "as it can be seen that highest eigen value is 4 \n\n the eigen vector is %d %d %d" %(X7[0],X7[1],X7[2])

X0=
[[1]
 [1]
 [0]]
X1=
[[7]
 [3]
 [0]]
X2=
[[25]
 [13]
 [ 0]]
X3=
[[8]
 [4]
 [0]]
X5=
[[8]
 [4]
 [0]]
X7=
[[2]
 [1]
 [0]]
as it can be seen that highest eigen value is 4 

 the eigen vector is 2 1 0

Ex7.17:pg-290

In [35]:
#housrholder's method
#example 7.17
#page 290
from numpy import matrix
from __future__ import division
import math
A=[[1, 3, 4],[3, 2, -1],[4, -1, 1]]
print A[1][1]
S=math.sqrt(A[0][1]**2+A[0][2]**2)
v2=math.sqrt((1+A[0][1]/S)/2)
v3=A[0][2]/(2*S)
v3=v3/v2
V=matrix([[0],[v2],[v3]])
P1=matrix([[1, 0, 0],[0, 1-2*v2**2, -2*v2*v3],[0, -2*v2*v3, 1-2*v3**2]])
A1=P1*A*P1
print "the reduced matrix is \n\n"
print A1

2
the reduced matrix is 


[[  1.00000000e+00  -5.00000000e+00  -8.88178420e-16]
 [ -5.00000000e+00   4.00000000e-01   2.00000000e-01]
 [ -8.88178420e-16   2.00000000e-01   2.60000000e+00]]