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]]