A question on symbolic Matrices - unexpected Decimals in algebraic entry

asked 2012-05-01 07:22:55 +0100

BrianLoudon gravatar image

Firstly, I'd just like to say that I am starting out with Sage adn Python at the moment, really impressed by the power of Sage and its potential! I am coming at this from an engineering and C++ background, and am currently trying to look at an engineering problem using some matrix linear algebra. I am setting up a from of stiffness matrix and am trying to extract the symbolic eigenvalues and eigenvectors from this. The stiffness matrix below is generated through a few different matrix multiplication operations, based on node and connectivity matrices. Despite trying to define these as being "SR" rings, I seem to end up with entries like (b + 1.0g + 1.0s) rather than the (b + g + s) I was expecting. I suspect this is due to me not understanding rings properly and how to apply them. I have tried the explanations in the Sage documentation, can anybody recommend some further reading or tutorials on this, because it does seem to be a vital topic in getting to grips with Sage?

I enclose my code below, in case anybody can spot a glaringly stupid thing that I've done...

many thanks,

Brian

N = matrix(QQ,[[1,0,-1,0],[0,1,0,-1]])

C = matrix(QQ,[[-1,0,1,0],[0,-1,0,1],[-1,1,0,0],[0,-1,1,0],[0,0,-1,1],[1,0,0,-1]])

D = C.transpose()

space = N.nrows()


M = N*D
M
II = identity_matrix(QQ,space)
b,s,g = var('b s g')
numMembers = M.ncols()

LVector  = II

#loop through members to create LVector
for i in range(numMembers):


    m1 = M.matrix_from_columns([i])
    mt = m1.transpose()
    mag = m1.norm()
    Coeff = m1*mt/(mag^2)
    size = m1.norm()

    if i<2:
        #for bars
        L = -g*(II - Coeff) + b*Coeff
    else:
        L = g*(II - Coeff) + s*Coeff

    L.factor()

    if i == 0:
        LVector = L
    else:
        LVector = LVector.block_sum(L)
    print L.str()

print LVector.str()

Pre = D.tensor_product(II)
Post = C.tensor_product(II)
K = Pre*LVector
K = K*Post

Ks = K.change_ring(SR)
print Ks.str()

Eigen = K.eigenvectors_left()
print Eigen.str()
edit retag flag offensive close merge delete

Comments

2

I don't have time to test this, but my guess is that mag.norm() returns something irrational, and every calculation involving it automatically converts to a real number, even if the irrational bit cancels later. In other words, you would expect 1.5/1.5 to be 1, but it is actually 1.0. And I wouldn't expect conversion to SR to convert reals to rationals; SR should support reals. But I am only guessing.

chaesloc2 gravatar imagechaesloc2 ( 2012-05-01 07:50:38 +0100 )edit

fantastic! thanks! I cast the result of mag^2 to an int and got the results I was expecting, great stuff! This has allowed me to verify my calc against an example and I can use this as a basis for what I'm trying to build up. Many thanks!

BrianLoudon gravatar imageBrianLoudon ( 2012-05-01 08:10:57 +0100 )edit

@chaesloc2; if you put your answer as an answer then Brian can accept it! :)

kcrisman gravatar imagekcrisman ( 2012-05-01 12:36:54 +0100 )edit