413 - Householder Symbolic Explicit example

1322 days ago by Professor413

A = matrix(QQ,[[2,0,4,7],[-1,3,5,2],[3,4,-1,-3],[-4,2,0,1]]) A 
       
[ 2  0  4  7]
[-1  3  5  2]
[ 3  4 -1 -3]
[-4  2  0  1]
[ 2  0  4  7]
[-1  3  5  2]
[ 3  4 -1 -3]
[-4  2  0  1]
x = A.transpose()[0].column() print(x) e1 = identity_matrix(4)[0].column() v = x + sqrt(x[0][0]^2+x[1][0]^2+x[2][0]^2+x[3][0]^2)*e1 print print(v) normsq = (v.transpose()*v)[0][0] print(normsq) H1 = identity_matrix(4)-(2/normsq)*v*v.transpose() show(H1) A1 = H1*A show(A1.simplify_full()) 
       
[ 2]
[-1]
[ 3]
[-4]

[sqrt(30) + 2]
[          -1]
[           3]
[          -4]
(sqrt(30) + 2)^2 + 26

                                
                            
[ 2]
[-1]
[ 3]
[-4]

[sqrt(30) + 2]
[          -1]
[           3]
[          -4]
(sqrt(30) + 2)^2 + 26

                                
x = A1.transpose()[1].column() # print(x) e2 = identity_matrix(4)[1].column() x[0] = 0 # get rid of the stuff above the diagonal on the resulting matrix v = (x + sqrt(x[0][0]^2+x[1][0]^2+x[2][0]^2+x[3][0]^2)*e2).simplify_full() # print(v) normsq = (v.transpose()*v)[0][0] print(normsq) H2 = (identity_matrix(4)-(2/normsq)*v*v.transpose()).simplify_full() # print(H2) A2 = (H2*A1).simplify_full() show(A2) 
       
1/4*(8*sqrt(30) + 117)^2/(sqrt(30) + 15)^2 + 4*(sqrt(30) +
16)^2/(sqrt(30) + 15)^2 + 1/900*(2*sqrt(869)*(sqrt(30) + 15) +
91*sqrt(30) + 180)^2/(sqrt(30) + 2)^2
1/4*(8*sqrt(30) + 117)^2/(sqrt(30) + 15)^2 + 4*(sqrt(30) + 16)^2/(sqrt(30) + 15)^2 + 1/900*(2*sqrt(869)*(sqrt(30) + 15) + 91*sqrt(30) + 180)^2/(sqrt(30) + 2)^2
x = A2.transpose()[2].column() # print(x) e3 = identity_matrix(4)[2].column() x[0] = 0 # get rid of the stuff above the diagonal on the resulting matrix x[1] = 0 print(x[1][0]) v = (x + sqrt(x[0][0]^2+x[1][0]^2+x[2][0]^2+x[3][0]^2)*e3).simplify_full() print(n(v)) normsq = (v.transpose()*v)[0][0] print(normsq) H3 = (identity_matrix(4)-(2/normsq)*v*v.transpose()).simplify_full() # print(H3) A3 = (H3*A2).simplify_full() show(n(A3)) 
       
0
[  0.000000000000000]
[  0.000000000000000]
[0.00503058044563949]
[  0.248698744891207]
4*(sqrt(869)*(sqrt(79)*(754617436114418005990290661*sqrt(30) +
4133209920515324789529634200) -
3*sqrt(83)*(237450130239941837016927953*sqrt(30) +
1300567926427990058397130155)) +
11*sqrt(79)*(1976220427111604587721903923*sqrt(30) +
10824205065269680533138795825) -
2607*sqrt(83)*(8418336037819938344170181*sqrt(30) +
46109125435874690854944360))^2/(sqrt(869)*sqrt(79)*(23745013023994183701\
6927953*sqrt(30) + 1300567926427990058397130155) +
869*sqrt(79)*(8418336037819938344170181*sqrt(30) +
46109125435874690854944360))^2 +
16*(sqrt(869)*(8071784101750996010827993*sqrt(30) +
44210981885743553923674510) + 652101505389614136285206074*sqrt(30) +
3571707055859153537822790030)^2/(sqrt(869)*(237450130239941837016927953*\
sqrt(30) + 1300567926427990058397130155) +
7315534016865526421083887289*sqrt(30) + 40068830003775106352946648840)^2
0
[  0.000000000000000]
[  0.000000000000000]
[0.00503058044563949]
[  0.248698744891207]
4*(sqrt(869)*(sqrt(79)*(754617436114418005990290661*sqrt(30) + 4133209920515324789529634200) - 3*sqrt(83)*(237450130239941837016927953*sqrt(30) + 1300567926427990058397130155)) + 11*sqrt(79)*(1976220427111604587721903923*sqrt(30) + 10824205065269680533138795825) - 2607*sqrt(83)*(8418336037819938344170181*sqrt(30) + 46109125435874690854944360))^2/(sqrt(869)*sqrt(79)*(237450130239941837016927953*sqrt(30) + 1300567926427990058397130155) + 869*sqrt(79)*(8418336037819938344170181*sqrt(30) + 46109125435874690854944360))^2 + 16*(sqrt(869)*(8071784101750996010827993*sqrt(30) + 44210981885743553923674510) + 652101505389614136285206074*sqrt(30) + 3571707055859153537822790030)^2/(sqrt(869)*(237450130239941837016927953*sqrt(30) + 1300567926427990058397130155) + 7315534016865526421083887289*sqrt(30) + 40068830003775106352946648840)^2