Python function.expand函数代码示例

def test(As,Bs,Cs,Ds,Us,sI_A,Js,Ys,B0,C0,C,Psi,Lamda):

    D1=BlockMatrix([[Us.inv()*Psi,zeros(Us.rows,1)],[-Lamda,eye(Lamda.rows)]]).as_mutable()
    D2=BlockMatrix([[sI_A,B0],[-C,Js]]).as_mutable()
    D3=BlockMatrix([[As,Bs],[-Cs,Ds]]).as_mutable()
    D4=BlockMatrix([[C0,-Ys],[zeros(Ys.cols,C0.cols),eye(Ys.cols)]]).as_mutable()
    return expand(simplify(D1*D2))==expand(D3*D4)

    def _eval_simplify(self, ratio=1.7, measure=None):
        from sympy.simplify.simplify import factor_sum, sum_combine
        from sympy.core.function import expand
        from sympy.core.mul import Mul

        # split the function into adds
        terms = Add.make_args(expand(self.function))
        s_t = [] # Sum Terms
        o_t = [] # Other Terms

        for term in terms:
            if term.has(Sum):
                # if there is an embedded sum here
                # it is of the form x * (Sum(whatever))
                # hence we make a Mul out of it, and simplify all interior sum terms
                subterms = Mul.make_args(expand(term))
                out_terms = []
                for subterm in subterms:
                    # go through each term
                    if isinstance(subterm, Sum):
                        # if it's a sum, simplify it
                        out_terms.append(subterm._eval_simplify())
                    else:
                        # otherwise, add it as is
                        out_terms.append(subterm)

                # turn it back into a Mul
                s_t.append(Mul(*out_terms))
            else:
                o_t.append(term)

        # next try to combine any interior sums for further simplification
        result = Add(sum_combine(s_t), *o_t)

        return factor_sum(result, limits=self.limits)

def test_expand_frac():
    assert expand((x + y)*y/x/(x + 1), frac=True) == \
        (x*y + y**2)/(x**2 + x)
    assert expand((x + y)*y/x/(x + 1), numer=True) == \
        (x*y + y**2)/(x*(x + 1))
    assert expand((x + y)*y/x/(x + 1), denom=True) == \
        y*(x + y)/(x**2 + x)
    eq = (x + 1)**2/y
    assert expand_numer(eq, multinomial=False) == eq

def _add_splines(c, b1, d, b2):
    """Construct c*b1 + d*b2."""
    if b1 == S.Zero or c == S.Zero:
        return expand(piecewise_fold(d * b2))
    if b2 == S.Zero or d == S.Zero:
        return expand(piecewise_fold(c * b1))
    new_args = []
    n_intervals = len(b1.args)
    assert n_intervals == len(b2.args)
    new_args.append((expand(c * b1.args[0].expr), b1.args[0].cond))
    for i in range(1, n_intervals - 1):
        new_args.append((expand(c * b1.args[i].expr + d * b2.args[i - 1].expr), b1.args[i].cond))
    new_args.append((expand(d * b2.args[-2].expr), b2.args[-2].cond))
    new_args.append(b2.args[-1])
    return Piecewise(*new_args)

def _qsimplify_fermion_product_site(e):
    """
    Assumes that all fermion operators have the same name.
    """
    work = [e]
    out = []
    while len(work):
        e = work.pop()
        if not isinstance(e, Mul):
            out.append(e)
            continue

        c, nc = e.args_cnc()
        ann = [f.is_annihilation for f in nc]

        # Find the first inverted pair of operators
        idx = next((i for i in range(len(ann)) if ann[i:i+2] == [True,False]), None)
        if idx is None:
            out.append(e)
            continue

        # Substitute the inverted pair for a cannonical subexpression
        e = Mul(*c) * (Mul(*nc[:idx])
                       * (Integer(1) - Mul(nc[idx+1],nc[idx]))
                       * Mul(*nc[idx+2:]))

        # Expand the new as a sum of fermion operator strings
        e = expand(e)

        # Recursive simplify
        if isinstance(e, Add):
            work.extend(e.args)
        else:
            work.append(e)
    return Add(*out)

def highest_row_degree_matrix(T,s):
    """
    Returns a matrix containing with the coefficients of the a polynomial 
    including those which the coefficient is the Zero matrix
    
    computes the highest row degree coefficient matrix of poynomial matrix A
    A : polynomial Matrix
    D :degree of matrix up to which we want to complete the list with zeroes
    Example:  TODO
    """
    s=symbols('s')
    RD=row_degrees(T,s)
    return  Matrix(T.rows,T.cols, lambda i,j: expand(T[i,j]).coeff(s,RD[i]))

def test_issues_5919_6830():
    # issue 5919
    n = -1 + 1/x
    z = n/x/(-n)**2 - 1/n/x
    assert expand(z) == 1/(x**2 - 2*x + 1) - 1/(x - 2 + 1/x) - 1/(-x + 1)

    # issue 6830
    p = (1 + x)**2
    assert expand_multinomial((1 + x*p)**2) == (
        x**2*(x**4 + 4*x**3 + 6*x**2 + 4*x + 1) + 2*x*(x**2 + 2*x + 1) + 1)
    assert expand_multinomial((1 + (y + x)*p)**2) == (
        2*((x + y)*(x**2 + 2*x + 1)) + (x**2 + 2*x*y + y**2)*
        (x**4 + 4*x**3 + 6*x**2 + 4*x + 1) + 1)
    A = Symbol('A', commutative=False)
    p = (1 + A)**2
    assert expand_multinomial((1 + x*p)**2) == (
        x**2*(1 + 4*A + 6*A**2 + 4*A**3 + A**4) + 2*x*(1 + 2*A + A**2) + 1)
    assert expand_multinomial((1 + (y + x)*p)**2) == (
        (x + y)*(1 + 2*A + A**2)*2 + (x**2 + 2*x*y + y**2)*
        (1 + 4*A + 6*A**2 + 4*A**3 + A**4) + 1)
    assert expand_multinomial((1 + (y + x)*p)**3) == (
        (x + y)*(1 + 2*A + A**2)*3 + (x**2 + 2*x*y + y**2)*(1 + 4*A +
        6*A**2 + 4*A**3 + A**4)*3 + (x**3 + 3*x**2*y + 3*x*y**2 + y**3)*(1 + 6*A
        + 15*A**2 + 20*A**3 + 15*A**4 + 6*A**5 + A**6) + 1)
    # unevaluate powers
    eq = (Pow((x + 1)*((A + 1)**2), 2, evaluate=False))
    # - in this case the base is not an Add so no further
    #   expansion is done
    assert expand_multinomial(eq) == \
        (x**2 + 2*x + 1)*(1 + 4*A + 6*A**2 + 4*A**3 + A**4)
    # - but here, the expanded base *is* an Add so it gets expanded
    eq = (Pow(((A + 1)**2), 2, evaluate=False))
    assert expand_multinomial(eq) == 1 + 4*A + 6*A**2 + 4*A**3 + A**4

    # coverage
    def ok(a, b, n):
        e = (a + I*b)**n
        return verify_numerically(e, expand_multinomial(e))

    for a in [2, S.Half]:
        for b in [3, S(1)/3]:
            for n in range(2, 6):
                assert ok(a, b, n)

    assert expand_multinomial((x + 1 + O(z))**2) == \
        1 + 2*x + x**2 + O(z)
    assert expand_multinomial((x + 1 + O(z))**3) == \
        1 + 3*x + 3*x**2 + x**3 + O(z)

    assert expand_multinomial(3**(x + y + 3)) == 27*3**(x + y)

def sum_simplify(s):
    """Main function for Sum simplification"""
    from sympy.concrete.summations import Sum
    from sympy.core.function import expand

    terms = Add.make_args(expand(s))
    s_t = [] # Sum Terms
    o_t = [] # Other Terms

    for term in terms:
        if isinstance(term, Mul):
            other = 1
            sum_terms = []

            if not term.has(Sum):
                o_t.append(term)
                continue

            mul_terms = Mul.make_args(term)
            for mul_term in mul_terms:
                if isinstance(mul_term, Sum):
                    r = mul_term._eval_simplify()
                    sum_terms.extend(Add.make_args(r))
                else:
                    other = other * mul_term
            if len(sum_terms):
                #some simplification may have happened
                #use if so
                s_t.append(Mul(*sum_terms) * other)
            else:
                o_t.append(other)
        elif isinstance(term, Sum):
            #as above, we need to turn this into an add list
            r = term._eval_simplify()
            s_t.extend(Add.make_args(r))
        else:
            o_t.append(term)


    result = Add(sum_combine(s_t), *o_t)

    return result

def test_expand_non_commutative():
    A = Symbol('A', commutative=False)
    B = Symbol('B', commutative=False)
    C = Symbol('C', commutative=False)
    a = Symbol('a')
    b = Symbol('b')
    i = Symbol('i', integer=True)
    n = Symbol('n', negative=True)
    m = Symbol('m', negative=True)
    p = Symbol('p', polar=True)
    np = Symbol('p', polar=False)

    assert (C*(A + B)).expand() == C*A + C*B
    assert (C*(A + B)).expand() != A*C + B*C
    assert ((A + B)**2).expand() == A**2 + A*B + B*A + B**2
    assert ((A + B)**3).expand() == (A**2*B + B**2*A + A*B**2 + B*A**2 +
                                     A**3 + B**3 + A*B*A + B*A*B)
    # issue 6219
    assert ((a*A*B*A**-1)**2).expand() == a**2*A*B**2/A
    # Note that (a*A*B*A**-1)**2 is automatically converted to a**2*(A*B*A**-1)**2
    assert ((a*A*B*A**-1)**2).expand(deep=False) == a**2*(A*B*A**-1)**2
    assert ((a*A*B*A**-1)**2).expand() == a**2*(A*B**2*A**-1)
    assert ((a*A*B*A**-1)**2).expand(force=True) == a**2*A*B**2*A**(-1)
    assert ((a*A*B)**2).expand() == a**2*A*B*A*B
    assert ((a*A)**2).expand() == a**2*A**2
    assert ((a*A*B)**i).expand() == a**i*(A*B)**i
    assert ((a*A*(B*(A*B/A)**2))**i).expand() == a**i*(A*B*A*B**2/A)**i
    # issue 6558
    assert (A*B*(A*B)**-1).expand() == A*B*(A*B)**-1
    assert ((a*A)**i).expand() == a**i*A**i
    assert ((a*A*B*A**-1)**3).expand() == a**3*A*B**3/A
    assert ((a*A*B*A*B/A)**3).expand() == \
        a**3*A*B*(A*B**2)*(A*B**2)*A*B*A**(-1)
    assert ((a*A*B*A*B/A)**-3).expand() == \
        a**-3*(A*B*(A*B**2)*(A*B**2)*A*B*A**(-1))**-1
    assert ((a*b*A*B*A**-1)**i).expand() == a**i*b**i*(A*B/A)**i
    assert ((a*(a*b)**i)**i).expand() == a**i*a**(i**2)*b**(i**2)
    e = Pow(Mul(a, 1/a, A, B, evaluate=False), S(2), evaluate=False)
    assert e.expand() == A*B*A*B
    assert sqrt(a*(A*b)**i).expand() == sqrt(a*b**i*A**i)
    assert (sqrt(-a)**a).expand() == sqrt(-a)**a
    assert expand((-2*n)**(i/3)) == 2**(i/3)*(-n)**(i/3)
    assert expand((-2*n*m)**(i/a)) == (-2)**(i/a)*(-n)**(i/a)*(-m)**(i/a)
    assert expand((-2*a*p)**b) == 2**b*p**b*(-a)**b
    assert expand((-2*a*np)**b) == 2**b*(-a*np)**b
    assert expand(sqrt(A*B)) == sqrt(A*B)
    assert expand(sqrt(-2*a*b)) == sqrt(2)*sqrt(-a*b)

def classify_pde(eq, func=None, dict=False, **kwargs):
    """
    Returns a tuple of possible pdsolve() classifications for a PDE.

    The tuple is ordered so that first item is the classification that
    pdsolve() uses to solve the PDE by default.  In general,
    classifications at the near the beginning of the list will produce
    better solutions faster than those near the end, thought there are
    always exceptions.  To make pdsolve use a different classification,
    use pdsolve(PDE, func, hint=<classification>).  See also the pdsolve()
    docstring for different meta-hints you can use.

    If ``dict`` is true, classify_pde() will return a dictionary of
    hint:match expression terms. This is intended for internal use by
    pdsolve().  Note that because dictionaries are ordered arbitrarily,
    this will most likely not be in the same order as the tuple.

    You can get help on different hints by doing help(pde.pde_hintname),
    where hintname is the name of the hint without "_Integral".

    See sympy.pde.allhints or the sympy.pde docstring for a list of all
    supported hints that can be returned from classify_pde.


    Examples
    ========
    >>> from sympy.solvers.pde import classify_pde
    >>> from sympy import Function, diff, Eq
    >>> from sympy.abc import x, y
    >>> f = Function('f')
    >>> u = f(x, y)
    >>> ux = u.diff(x)
    >>> uy = u.diff(y)
    >>> eq = Eq(1 + (2*(ux/u)) + (3*(uy/u)))
    >>> classify_pde(eq)
    ('1st_linear_constant_coeff_homogeneous',)
    """

    prep = kwargs.pop('prep', True)

    if func and len(func.args) != 2:
        raise NotImplementedError("Right now only partial "
            "differential equations of two variables are supported")

    if prep or func is None:
        prep, func_ = _preprocess(eq, func)
        if func is None:
            func = func_

    if isinstance(eq, Equality):
        if eq.rhs != 0:
            return classify_pde(eq.lhs - eq.rhs, func)
        eq = eq.lhs

    f = func.func
    x = func.args[0]
    y = func.args[1]
    fx = f(x,y).diff(x)
    fy = f(x,y).diff(y)

    # TODO : For now pde.py uses support offered by the ode_order function
    # to find the order with respect to a multi-variable function. An
    # improvement could be to classify the order of the PDE on the basis of
    # individual variables.
    order = ode_order(eq, f(x,y))

    # hint:matchdict or hint:(tuple of matchdicts)
    # Also will contain "default":<default hint> and "order":order items.
    matching_hints = {'order': order}

    if not order:
        if dict:
            matching_hints["default"] = None
            return matching_hints
        else:
            return ()

    eq = expand(eq)

    a = Wild('a', exclude = [f(x,y)])
    b = Wild('b', exclude = [f(x,y), fx, fy, x, y])
    c = Wild('c', exclude = [f(x,y), fx, fy, x, y])
    d = Wild('d', exclude = [f(x,y), fx, fy, x, y])
    e = Wild('e', exclude = [f(x,y), fx, fy])
    n = Wild('n', exclude = [x, y])
    # Try removing the smallest power of f(x,y)
    # from the highest partial derivatives of f(x,y)
    reduced_eq = None
    if eq.is_Add:
        var = set(combinations_with_replacement((x,y), order))
        dummyvar = deepcopy(var)
        power = None
        for i in var:
            coeff = eq.coeff(f(x,y).diff(*i))
            if coeff != 1:
                match = coeff.match(a*f(x,y)**n)
                if match and match[a]:
                    power = match[n]
                    dummyvar.remove(i)
                    break
#.........这里部分代码省略.........

def solveset_complex(f, symbol):
    """ Solve a complex valued equation.

    Parameters
    ==========

    f : Expr
        The target equation
    symbol : Symbol
        The variable for which the equation is solved

    Returns
    =======

    Set
        A set of values for `symbol` for which `f` equal to
        zero. An `EmptySet` is returned if no solution is found.

    `solveset_complex` claims to be complete in the solution set that
    it returns.

    Raises
    ======

    NotImplementedError
        The algorithms for to find the solution of the given equation are
        not yet implemented.
    ValueError
        The input is not valid.
    RuntimeError
        It is a bug, please report to the github issue tracker.

    See Also
    ========

    solveset_real: solver for real domain

    Examples
    ========

    >>> from sympy import Symbol, exp
    >>> from sympy.solvers.solveset import solveset_complex
    >>> from sympy.abc import x, a, b, c
    >>> solveset_complex(a*x**2 + b*x +c, x)
    {-b/(2*a) - sqrt(-4*a*c + b**2)/(2*a), -b/(2*a) + sqrt(-4*a*c + b**2)/(2*a)}

    Due to the fact that complex extension of my real valued functions are
    multivariate even some simple equations can have infinitely many solution.

    >>> solveset_complex(exp(x) - 1, x)
    ImageSet(Lambda(_n, 2*_n*I*pi), Integers())

    """
    if not symbol.is_Symbol:
        raise ValueError(" %s is not a symbol" % (symbol))

    f = sympify(f)
    original_eq = f
    if not isinstance(f, (Expr, Number)):
        raise ValueError(" %s is not a valid sympy expression" % (f))

    f = together(f)
    # Without this equations like a + 4*x**2 - E keep oscillating
    # into form  a/4 + x**2 - E/4 and (a + 4*x**2 - E)/4
    if not fraction(f)[1].has(symbol):
        f = expand(f)

    if f.is_zero:
        raise NotImplementedError("S.Complex set is not yet implemented")
    elif not f.has(symbol):
        result = EmptySet()
    elif f.is_Mul and all([_is_finite_with_finite_vars(m) for m in f.args]):
        result = Union(*[solveset_complex(m, symbol) for m in f.args])
    else:
        lhs, rhs_s = invert_complex(f, 0, symbol)
        if lhs == symbol:
            result = rhs_s
        elif isinstance(rhs_s, FiniteSet):
            equations = [lhs - rhs for rhs in rhs_s]
            result = EmptySet()
            for equation in equations:
                if equation == f:
                    result += _solve_as_rational(equation, symbol,
                                                 solveset_solver=solveset_complex,
                                                 as_poly_solver=_solve_as_poly_complex)
                else:
                    result += solveset_complex(equation, symbol)
        else:
            raise NotImplementedError

    if isinstance(result, FiniteSet):
        result = [s for s in result
                  if isinstance(s, RootOf)
                  or domain_check(original_eq, symbol, s)]
        return FiniteSet(*result)
    else:
        return result

def test_expand_log():
    t = Symbol('t', positive=True)
    # after first expansion, -2*log(2) + log(4); then 0 after second
    assert expand(log(t**2) - log(t**2/4) - 2*log(2)) == 0

def solveset_real(f, symbol):
    """ Solves a real valued equation.

    Parameters
    ==========

    f : Expr
        The target equation
    symbol : Symbol
        The variable for which the equation is solved

    Returns
    =======

    Set
        A set of values for `symbol` for which `f` is equal to
        zero. An `EmptySet` is returned if no solution is found.
        A `ConditionSet` is returned as unsolved object if algorithms
        to evaluate complete solutions are not yet implemented.

    `solveset_real` claims to be complete in the set of the solution it
    returns.

    Raises
    ======

    NotImplementedError
        Algorithms to solve inequalities in complex domain are
        not yet implemented.
    ValueError
        The input is not valid.
    RuntimeError
        It is a bug, please report to the github issue tracker.


    See Also
    =======

    solveset_complex : solver for complex domain

    Examples
    ========

    >>> from sympy import Symbol, exp, sin, sqrt, I
    >>> from sympy.solvers.solveset import solveset_real
    >>> x = Symbol('x', real=True)
    >>> a = Symbol('a', real=True, finite=True, positive=True)
    >>> solveset_real(x**2 - 1, x)
    {-1, 1}
    >>> solveset_real(sqrt(5*x + 6) - 2 - x, x)
    {-1, 2}
    >>> solveset_real(x - I, x)
    EmptySet()
    >>> solveset_real(x - a, x)
    {a}
    >>> solveset_real(exp(x) - a, x)
    {log(a)}

    * In case the equation has infinitely many solutions an infinitely indexed
      `ImageSet` is returned.

    >>> solveset_real(sin(x) - 1, x)
    ImageSet(Lambda(_n, 2*_n*pi + pi/2), Integers())

    * If the equation is true for any arbitrary value of the symbol a `S.Reals`
      set is returned.

    >>> solveset_real(x - x, x)
    (-oo, oo)

    """
    if not symbol.is_Symbol:
        raise ValueError(" %s is not a symbol" % (symbol))

    f = sympify(f)
    if not isinstance(f, (Expr, Number)):
        raise ValueError(" %s is not a valid sympy expression" % (f))

    original_eq = f
    f = together(f)

    # In this, unlike in solveset_complex, expression should only
    # be expanded when fraction(f)[1] does not contain the symbol
    # for which we are solving
    if not symbol in fraction(f)[1].free_symbols and f.is_rational_function():
        f = expand(f)

    if f.has(Piecewise):
        f = piecewise_fold(f)
    result = EmptySet()

    if f.expand().is_zero:
        return S.Reals
    elif not f.has(symbol):
        return EmptySet()
    elif f.is_Mul and all([_is_finite_with_finite_vars(m) for m in f.args]):
        # if f(x) and g(x) are both finite we can say that the solution of
        # f(x)*g(x) == 0 is same as Union(f(x) == 0, g(x) == 0) is not true in
        # general. g(x) can grow to infinitely large for the values where
        # f(x) == 0. To be sure that we are not silently allowing any
#.........这里部分代码省略.........

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