美文网首页
sympy计算smith标准型

sympy计算smith标准型

作者: 一路向后 | 来源:发表于2020-12-24 22:33 被阅读0次

    1.源码实现

    import itertools
    from sympy import *
    
    def nchoosek(a, b, d=1, n=1):
        c = []
        for i in itertools.combinations(range(a,b+1,d),n):
            c.append(list(i))
        return c
    
    def smith_form(C):
        #1---求lambda矩阵的行列式因子(开始)
        n = C.shape[0]
        h = [];
    
        for k in range(1, n+1):
            y = [];
    
            #从n个元素(第一个元素到第n个元素的所有元素)中取出k项,其中所有可能的组合放在p中(组合数)
            p = Matrix(nchoosek(1, n, n=k));
            q = p;
            m = p.shape[0]
    
            #print("p=")
            #print(p)
    
            #1-b将所有的i阶行列式求出
            for i in range(0, m):
                a = p.row(i)
    
                for j in range(0, m):
                    b = q.row(j)
                    B = []
    
                    for s in range(0, k):
                        g = [];
    
                        for t in range(0, k):
                            u = a[s]
                            v = b[t]
                            g.append(C[u-1, v-1])
    
                        B.append(g)
    
                    c = Matrix(B).det()
    
                    if c != 0:
                        y.append(c);
            #1-b结束
            o = len(y)
    
            #1-a--求所有j阶行列式的最大公因子(开始)
            for j in range(0, o-1):
                y[j+1] = gcd(y[j],y[j+1]);
            #1-a--结束
    
            h.append(y[o-1]);
        #1---求lambda矩阵的行列式因子(结束)
    
        d = [];
        d.append(h[0])
    
        for i in range(1, n):
            d.append(expand(h[i]/h[i-1]));
    
        #2---将上面所求行列式因子进行因式分解,得到因式分解后的行列式因子
        k = len(d)
    
        for i in range(0, k):
            t = d[i]
            t = factor(t)
            d[i] = t
    
        return diag(d)
    
    
    #声明变量x
    x = symbols('x')
    
    #输入lambda矩阵
    C = Matrix([[-x+1, 2*x-1, x], [x, x**2, -x], [x**2+1, x**2+x-1, -x**2]])
    
    #计算smith标准型
    D = smith_form(C)
    
    print(C)
    print(D)
    

    2.运行及结果

    $ python3 example.py
    Matrix([[-x + 1, 2*x - 1, x], [x, x**2, -x], [x**2 + 1, x**2 + x - 1, -x**2]])
    Matrix([[1], [x], [-x*(x**2 + 1)]])
    

    相关文章

      网友评论

          本文标题:sympy计算smith标准型

          本文链接:https://www.haomeiwen.com/subject/njpvnktx.html