Python Matrix.subs方法代码示例

# 需要导入模块: from sympy.matrices import Matrix [as 别名]
# 或者: from sympy.matrices.Matrix import subs [as 别名]
def ALGO4(As,Bs,Cs,Ds,do_test):
    
    
    
#-------------------STEP 1------------------------------
    if not is_row_proper(As):
        Us,Ast=row_proper(As)
        
    else:
        Us,Ast=eye(As.rows),As
    
    Bst=Us*Bs
    Bst=expand(Bst)
    r=Ast.cols
    #-------------------STEP 2------------------------------
    K=simplify(Ast.inv()*Bst)  #very important 
    Ys=zeros(K.shape)
    for i,j in  product(range(K.rows),range(K.cols)):
        Ys[i,j],q=div(numer(K[i,j]),denom(K[i,j]))       
    
    B_hat=Bst-Ast*Ys
    
    #-------------------END STEP 2------------------------------
    #-------------------STEP 3------------------------------
    Psi=diag(*[[s**( mc.row_degrees(Ast,s)[j]  -i -1) for i in range( mc.row_degrees(Ast,s)[j])] for j in range(r)]).T
    S=diag(*[s**(rho) for rho in mc.row_degrees(Ast,s)])
    Ahr=mc.highest_row_degree_matrix(Ast,s)
    Help=Ast-S*Ahr
    
    SOL={}
    numvar=Psi.rows*Psi.cols
    alr=symbols('a0:%d'%numvar)
    Alr=Matrix(Psi.cols,Psi.rows,alr)
    RHS=Psi*Alr
    for i,j in  product(range(Help.rows),range(Help.cols)):                 #diagonal explain later
        SOL.update(solve_undetermined_coeffs(Eq(Help[i,j],RHS[i,j]),alr,s))
    
    Alr=Alr.subs(SOL)    #substitute(SOL)
    
    
    Aoc=Matrix(BlockDiagMatrix(*[Matrix(rho, rho, lambda i,j: KroneckerDelta(i+1,j))for rho in mc.row_degrees(Ast,s)]))
    Boc=eye(sum(mc.row_degrees(Ast,s)))
    Coc=Matrix(BlockDiagMatrix(*[SparseMatrix(1,rho,{(x,0):1 if x==0 else 0 for x in range(rho)}) for rho in mc.row_degrees(Ast,s)]))

    A0=Aoc-Alr*Ahr.inv()*Matrix(Coc)
    C0=Ahr.inv()*Coc



    SOL={}
    numvar=Psi.cols*Bst.cols
    b0=symbols('b0:%d'%numvar)
    B0=Matrix(Psi.cols,Bst.cols,b0)
    RHS=Psi*B0

    for i,j in  product(range(B_hat.rows),range(B_hat.cols)):                 #diagonal explain later
        SOL.update(solve_undetermined_coeffs(Eq(B_hat[i,j],RHS[i,j]),b0,s))
    B0=B0.subs(SOL)    #substitute(SOL)

    LHS_matrix=simplify(Cs*C0)                                        #left hand side of the equation (1)

    sI_A=s*eye(A0.cols)- A0
    max_degree=mc.find_degree(LHS_matrix,s)                   #get the degree of the matrix at the LHS
                                                  #which is also the maximum degree for the coefficients of Λ(s)
    #---------------------------Creating Matrices Λ(s) and C -------------------------------------

    Lamda=[]
    numvar=((max_degree))*A0.cols
    a=symbols('a0:%d'%numvar)
    for i in range(A0.cols):                                        # paratirisi den douleuei to prin giat;i otra oxi diagonios
        p=sum(a[n +i*(max_degree)]*s**n for n in range(max_degree))  # we want variables one degree lower because we are multiplying by first order monomials
        Lamda.append(p)
        
    Lamda=Matrix(Cs.rows,A0.cols,Lamda)                                #convert the list to Matrix
    
    c=symbols('c0:%d'%(Lamda.rows*Lamda.cols))
    C=Matrix(Lamda.rows,Lamda.cols,c)
    #-----------------------------------------
    
    RHS_matrix=Lamda*sI_A +C                            #right hand side of the equation (1)
     
    '''
    -----------Converting equation (1) to a system of linear -----------
    -----------equations, comparing the coefficients of the  -----------
    -----------polynomials in both sides of the equation (1) -----------
    '''
     EQ=[Eq(LHS_matrix[i,j],expand(RHS_matrix[i,j])) for i,j in product(range(LHS_matrix.rows),range(LHS_matrix.cols)) ]

# 需要导入模块: from sympy.matrices import Matrix [as 别名]
# 或者: from sympy.matrices.Matrix import subs [as 别名]
def ALGO21(As,Bs,Cs,Ds,do_test):
    
    #-------------------STEP 1------------------------------
    r=As.rows
    p=Cs.rows
    m=Ds.cols

    Ts=BlockMatrix([[As,Bs,zeros(As.rows,p)],
                    [-Cs,Ds,eye(p)],
                    [zeros(m,Cs.cols),-eye(m),zeros(m,p)]]).as_mutable()
                  
    rt=r+p+m
    T_over_I=BlockMatrix([[Ts],[eye(rt)]]).as_mutable()
    H=reduction.column_proper(T_over_I)
    Qs_over_Rs=BlockMatrix([[H[1][0:rt,0:rt]],[H[0]]]).as_mutable()    
    #-------------------END STEP 1------------------------------
    
    #-------------------STEP 2------------------------------
    qi=mc.column_degrees(Qs_over_Rs,s)
    Ss=BlockDiagMatrix(*[Matrix(q+1,1,[s**(q-k) for k in range(q+1)]) for q in qi]).as_mutable() ## POS??POSS?   
    
    #find Qc
    numvar=rt*Ss.rows
    a=symbols('a0:%d'%numvar)
    Qc=Matrix(rt,Ss.rows,a)
    RHS_matrix=(Qc*Ss).as_mutable()
    Qs=Qs_over_Rs[0:rt,0:rt]

    SOL=dict()
    for i,j in  product(range(Qs.rows),range(Qs.cols)):                 
        SOL.update(solve_undetermined_coeffs(Eq(Qs[i,j],RHS_matrix[i,j]),a,s))
    Qc=Qc.subs(SOL) 

    #find Rc
    numvar=rt*Ss.rows
    Rc=Matrix(rt,Ss.rows,a)
    RHS_matrix=(Rc*Ss).as_mutable()
    Rs=H[0]
    SOL=dict()
    for i,j in  product(range(Rs.rows),range(Rs.cols)):                 
        SOL.update(solve_undetermined_coeffs(Eq(Rs[i,j],RHS_matrix[i,j]),a,s))
    Rc=Rc.subs(SOL) 
    #-------------------END STEP 2------------------------------
    
    #-------------------STEP 3------------------------------

    def E0block(i,j):
        if i==j and i!=0:
            return 1
        else:
            return 0
    
    E0=BlockDiagMatrix(*[Matrix(q+1,q+1, lambda i,j: E0block(i,j)) for q in qi]).as_mutable()
    A0=BlockDiagMatrix(*[Matrix(q+1,q+1, lambda i,j: KroneckerDelta(i-1,j)) for q in qi]).as_mutable()
    B0=BlockDiagMatrix(*[Matrix(q+1,1, lambda i,j: KroneckerDelta(i,j)) for q in qi]).as_mutable()
    n=sum(qi)+rt
    C0=eye(n)
    V=BlockMatrix([[zeros(p,r),zeros(p,m),eye(p)]]).as_mutable()
    U=BlockMatrix([[zeros(r,m)],[zeros(p,m)],[eye(m)]]).as_mutable()
    
    E=SparseMatrix(E0)
    A=SparseMatrix(A0-B0*Qc)
    B=SparseMatrix(B0*U)
    C=SparseMatrix(V*Rc)
    #E=E0.as_mutable()
    #pprint(B0*Qc)
    #A=(simplify(A0-B0*Qc)).as_mutable()
    #B=(B0*U).as_mutable()
    #C=(V*Rc).as_mutable()
     #-------------------END STEP 3------------------------------

    if do_test==True:
        test_result=test(B0,Ts,E,A,B,C,U,V,Rc,Rs,Ss,p,m,rt)
    else:
        test_result= 'not done' 
    
    D=zeros(p,m) 
    return E,A,B,C,D,test_result