Python sympy.expand函数代码示例

def main():
    print(symbols('a b c d'))
    a, b, c, d = symbols('a b c d')

    group4 = [a, b, c, d]
    print(group4)
    four = group_solve(group4)
    print(four)
    print(expand(four))
    print(four.collect(a))
    print(four.subs(a, 10))
    print()

    group4_ = add1(group4)
    print(group4_)
    four1 = group_solve(group4_)
    print(four1)
    print(expand(four1))
    print(four1.collect(a))
    print()

    print(13717421)

    for i in range(1, 10 ** 4):
        sum_power(i)

 def _solve_recur(expr, rhs=S(0)):
     if expr == var:
         return rhs
     expr = expand(expr)
     if isinstance(expr, Mul):
         lhs = S(1)
         # Try by rank
         l_coeff, var_expr, r_coeff = expr_coeff(expr, var)
         lhs = var_expr
         rhs = Inverse(l_coeff) * rhs * Inverse(r_coeff)
         return _solve_recur(lhs, rhs)
     if isinstance(expr, Add):
         lhs = 0
         for arg in expr.args:
             if var in arg:
                 lhs += arg
             else:
                 rhs -= arg
         if isinstance(lhs, Add):
             coeff = lhs.coeff(var)
             if expand(coeff * var) == lhs:
                 rhs /= coeff
                 lhs = var
         return _solve_recur(lhs, rhs)
     if isinstance(expr, Transpose):
         return _solve_recur(expr.args[0], Transpose(rhs))
     raise NotImplementedError("Can't handle expr of type %s" % type(expr))

def case0(f, N=3):
    B = 1 - x ** 3
    dBdx = sm.diff(B, x)

    # Compute basis functions and their derivatives
    phi = {0: [x ** (i + 1) * (1 - x) for i in range(N + 1)]}
    phi[1] = [sm.diff(phi_i, x) for phi_i in phi[0]]

    def integrand_lhs(phi, i, j):
        return phi[1][i] * phi[1][j]

    def integrand_rhs(phi, i):
        return f * phi[0][i] - dBdx * phi[1][i]

    Omega = [0, 1]

    u_bar = solve(integrand_lhs, integrand_rhs, phi, Omega, verbose=True, numint=False)
    u = B + u_bar
    print "solution u:", sm.simplify(sm.expand(u))

    # Calculate analytical solution

    # Solve -u''=f by integrating f twice
    f1 = sm.integrate(f, x)
    f2 = sm.integrate(f1, x)
    # Add integration constants
    C1, C2 = sm.symbols("C1 C2")
    u_e = -f2 + C1 * x + C2
    # Find C1 and C2 from the boundary conditions u(0)=0, u(1)=1
    s = sm.solve([u_e.subs(x, 0) - 1, u_e.subs(x, 1) - 0], [C1, C2])
    # Form the exact solution
    u_e = -f2 + s[C1] * x + s[C2]
    print "analytical solution:", u_e
    # print 'error:', u - u_e  # many terms - which cancel
    print "error:", sm.expand(u - u_e)

def alphaBetaToTransfer(adjMat, alphas, betas):
    s = sympy.Symbol('s')
    subs = {}

    # iterate over the columns of the adjacency matrix
    for i, col in enumerate(adjMat.getA()):
        expr = sympy.S(0)
        i += 1
        # iterate over the elements in each column of the adjacency matrix
        # this is to build the replacement expression for what were once "output to env"s
        # but are now the turnover rates of each node
        for j, elem in enumerate(col):
            j += 1
            if elem != 0:
                expr += (-sympy.var('a_%d%d' % (j,i)))

        subs[sympy.var('a_%d%d' % (i,i))] = expr

    #print "*** VAR->EXPR SUBSTITUTIONS: "
    #print subs

    for pos in range(0, len(betas)):
        orderedKeysBeta = getOrderedKeys(betas[pos])
        for key in orderedKeysBeta:
            betas[pos][s**key] = sympy.simplify(sympy.expand(betas[pos][s**key].xreplace(subs)))

        #finished all the betas so add the alphas, but only once
        orderedKeysAlpha = getOrderedKeys(alphas[pos])
        for key in orderedKeysAlpha:
            alphas[pos][s**key] = sympy.simplify(sympy.expand(alphas[pos][s**key].xreplace(subs)))

def test_one_dof():
    # This is for a 1 dof spring-mass-damper case.
    # It is described in more detail in the Kane docstring.
    q, u = dynamicsymbols('q u')
    qd, ud = dynamicsymbols('q u', 1)
    m, c, k = symbols('m c k')
    N = ReferenceFrame('N')
    P = Point('P')
    P.set_vel(N, u * N.x)

    kd = [qd - u]
    FL = [(P, (-k * q - c * u) * N.x)]
    pa = Particle()
    pa.mass = m
    pa.point = P
    BL = [pa]

    KM = Kane(N)
    KM.coords([q])
    KM.speeds([u])
    KM.kindiffeq(kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    assert expand(rhs[0]) == expand(-(q * k + u * c) / m)
    assert KM.linearize() == (Matrix([[0, 1], [k, c]]), Matrix([]))

def _interpolate_Function(expr):
    path = 'None'
    factor = 'None'
    change = False
    summand_0 = 'None'
    summand_1 = 'None'
    res = expr
    for i in np.arange(len(expr.args)):
        argument = sympy.expand(expr.args[i])
        if argument.func == sympy.Add:
            for j in np.arange(len(argument.args)):
                summand = argument.args[j]
                if summand.func == sympy.Mul:
                    for k in np.arange(len(summand.args)):
                        temp = 0
                        if summand.args[k] == sympy.Symbol('a'):
                            temp = sympy.Mul(sympy.Mul(*summand.args[:k]),
                                             sympy.Mul(*summand.args[k+1:]))
                            #print(temp)
                        if not temp == int(temp):
                            #print('found one')
                            factor = (temp)
                            path = np.array([i, j, k])
    if not factor == 'None':
        change = True

        sign = np.sign(factor)
        offsets = np.array([int(factor), sign * (int(sign * factor) + 1)])
        weights = 1/np.abs(offsets - factor)
        weights = weights/np.sum(weights)

        res = (  weights[0] * _interchange(expr, offsets[0] * sympy.Symbol('a'), path[:-1])
               + weights[1] * _interchange(expr, offsets[1] * sympy.Symbol('a'), path[:-1]))

    return sympy.expand(res), change

def _interpolate_expression(expr):
    change = False
    expr = sympy.expand(expr)
    res = expr

    if isinstance(expr, sympy.Function):# and not change:
        path = np.array([])
        temp, change = _interpolate_Function(_follow_path(expr, path))
        res = _interchange(expr, temp, path)

    for i in np.arange(len(expr.args)):
        path = np.array([i])
        if isinstance(_follow_path(expr, path), sympy.Function) and not change:
            temp, change = _interpolate_Function(_follow_path(expr, path))
            res = _interchange(expr, temp, path)

        for j in np.arange(len(expr.args[i].args)):
            path = np.array([i, j])
            if isinstance(_follow_path(expr, path), sympy.Function) and not change:
                temp, change = _interpolate_Function(_follow_path(expr, path))
                res = _interchange(expr, temp, path)

    if change:
        res = _interpolate_expression(res)

    return sympy.expand(res)

def test_SVCVS_laplace_d3_n1():
    """Test VCCS with a laplace defined transfer function with second order
    numerator and third order denominator
    """

    pycircuit.circuit.circuit.default_toolkit = symbolic
    cir = SubCircuit()

    n1,n2 = cir.add_nodes('1','2')

    b0,a0,a1,a2,a3,Gdc = [sympy.Symbol(symname, real=True) for
                                symname in 'b0,a0,a1,a2,a3,Gdc'
                                .split(',')]

    s = sympy.Symbol('s', complex=True)

    cir['VS']   = VS( n1, gnd, vac=1)
    cir['VCVS'] = SVCVS( n1, gnd, n2, gnd,
                        denominator = [a0, a1, a2, a3],
                        numerator   = [b0, 0, 0])

    res = AC(cir, toolkit=symbolic).solve(s, complexfreq=True)

    assert_equal(sympy.cancel(sympy.expand(res.v(n2,gnd))),
                 sympy.expand((b0*s*s)/(a0*s*s*s+a1*s*s+a2*s+a3)))

def tests():
    x, y, z = symbols('x,y,z')
    #print(x + x + 1)
    expr = x**2 - y**2
    factors = factor(expr)
    print(factors, " | ", expand(factors))
    pprint(expand(factors))

def test_one_dof():
    # This is for a 1 dof spring-mass-damper case.
    # It is described in more detail in the KanesMethod docstring.
    q, u = dynamicsymbols('q u')
    qd, ud = dynamicsymbols('q u', 1)
    m, c, k = symbols('m c k')
    N = ReferenceFrame('N')
    P = Point('P')
    P.set_vel(N, u * N.x)

    kd = [qd - u]
    FL = [(P, (-k * q - c * u) * N.x)]
    pa = Particle('pa', P, m)
    BL = [pa]

    KM = KanesMethod(N, [q], [u], kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    assert expand(rhs[0]) == expand(-(q * k + u * c) / m)
    assert (KM.linearize(A_and_B=True, new_method=True)[0] ==
            Matrix([[0, 1], [-k/m, -c/m]]))

    # Ensure that the old linearizer still works and that the new linearizer
    # gives the same results. The old linearizer is deprecated and should be
    # removed in >= 0.7.7.
    M_old = KM.mass_matrix_full
    # The old linearizer raises a deprecation warning, so catch it here so
    # it doesn't cause py.test to fail.
    with warnings.catch_warnings():
        warnings.filterwarnings("ignore", category=SymPyDeprecationWarning)
        F_A_old, F_B_old, r_old = KM.linearize()
    M_new, F_A_new, F_B_new, r_new = KM.linearize(new_method=True)
    assert simplify(M_new.inv() * F_A_new - M_old.inv() * F_A_old) == zeros(2)

def test_pend():
    q, u = dynamicsymbols('q u')
    qd, ud = dynamicsymbols('q u', 1)
    m, l, g = symbols('m l g')
    N = ReferenceFrame('N')
    P = Point('P')
    P.set_vel(N, -l * u * sin(q) * N.x + l * u * cos(q) * N.y)
    kd = [qd - u]

    FL = [(P, m * g * N.x)]
    pa = Particle()
    pa.mass = m
    pa.point = P
    BL = [pa]

    KM = Kane(N)
    KM.coords([q])
    KM.speeds([u])
    KM.kindiffeq(kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    rhs.simplify()
    assert expand(rhs[0]) == expand(-g / l * sin(q))

def T_exact_symbolic(verbose=False):
    """Compute the exact solution formula via sympy."""
    # sol1: solution for t < t_star,
    # sol2: solution for t > t_star
    import sympy as sym
    T0 = sym.symbols('T0')
    k = sym.symbols('k', positive=True)
    # Piecewise linear T_sunction
    t, t_star, C0, C1 = sym.symbols('t t_star C0 C1')
    T_s = C0
    I = sym.integrate(sym.exp(k*t)*T_s, (t, 0, t))
    sol1 = T0*sym.exp(-k*t) + k*sym.exp(-k*t)*I
    sol1 = sym.simplify(sym.expand(sol1))
    if verbose:
        # Some debugging print
        print 'solution t < t_star:', sol1
        #print sym.latex(sol1)
    T_s = C1
    I = sym.integrate(sym.exp(k*t)*C0, (t, 0, t_star)) + \
        sym.integrate(sym.exp(k*t)*C1, (t, t_star, t))
    sol2 = T0*sym.exp(-k*t) + k*sym.exp(-k*t)*I
    sol2 = sym.simplify(sym.expand(sol2))
    if verbose:
        print 'solution t > t_star:', sol2
        #print sym.latex(sol2)

    # Convert to numerical functions
    exact0 = sym.lambdify([t, C0, k, T0],
                          sol1, modules='numpy')
    exact1 = sym.lambdify([t, C0, C1, t_star, k, T0],
                          sol2, modules='numpy')
    return exact0, exact1

 def gen_errs(self, a, b, n, x, y):
     """
     What is the best thing to do here?
     """
     errs = set()
     
     # Generate som obvious candidates
     errs.add(sym.expand(self.gen_prob(-a, -b, n, x, y)))
     if a < 0:
         errs.add(sym.expand(self.gen_prob(-a, b, n, x, y)))
     if b < 0:
         errs.add(sym.expand(self.gen_prob(a, -b, n, x, y)))
    
     
     expr = sym.expand(self.gen_prob(a, b, n, x, y))
     coeffs = sym.Poly(expr).coeffs()
     
     for i in range(4):
         
         expr1 = self.gen_poly(map(lambda l: random.choice([-1, 1]) * l, coeffs), x, y)
         expr2 = self.gen_poly(map(lambda l: random.choice([-2, -1, 0, 1, 2]) * l, coeffs), x, y)
         errs.update([expr1, expr2, expr - (expr2 - sym.LM(expr))])
     
     errs = list(errs)
     errs_ = [err for err in errs if err != expr]
     random.shuffle(errs_)
     return errs

def test_two_dof():
    # This is for a 2 d.o.f., 2 particle spring-mass-damper.
    # The first coordinate is the displacement of the first particle, and the
    # second is the relative displacement between the first and second
    # particles. Speeds are defined as the time derivatives of the particles.
    q1, q2, u1, u2 = dynamicsymbols('q1 q2 u1 u2')
    q1d, q2d, u1d, u2d = dynamicsymbols('q1 q2 u1 u2', 1)
    m, c1, c2, k1, k2 = symbols('m c1 c2 k1 k2')
    N = ReferenceFrame('N')
    P1 = Point('P1')
    P2 = Point('P2')
    P1.set_vel(N, u1 * N.x)
    P2.set_vel(N, (u1 + u2) * N.x)
    kd = [q1d - u1, q2d - u2]

    # Now we create the list of forces, then assign properties to each
    # particle, then create a list of all particles.
    FL = [(P1, (-k1 * q1 - c1 * u1 + k2 * q2 + c2 * u2) * N.x), (P2, (-k2 *
        q2 - c2 * u2) * N.x)]
    pa1 = Particle('pa1', P1, m)
    pa2 = Particle('pa2', P2, m)
    BL = [pa1, pa2]

    # Finally we create the KanesMethod object, specify the inertial frame,
    # pass relevant information, and form Fr & Fr*. Then we calculate the mass
    # matrix and forcing terms, and finally solve for the udots.
    KM = KanesMethod(N, q_ind=[q1, q2], u_ind=[u1, u2], kd_eqs=kd)
    KM.kanes_equations(FL, BL)
    MM = KM.mass_matrix
    forcing = KM.forcing
    rhs = MM.inv() * forcing
    assert expand(rhs[0]) == expand((-k1 * q1 - c1 * u1 + k2 * q2 + c2 * u2)/m)
    assert expand(rhs[1]) == expand((k1 * q1 + c1 * u1 - 2 * k2 * q2 - 2 *
                                    c2 * u2) / m)

    def test_apply_substitutions(self):

        def apply_correct_substitutions(monomial, substitutions):
            if isinstance(monomial, int) or isinstance(monomial, float):
                return monomial
            original_monomial = monomial
            changed = True
            while changed:
                for lhs, rhs in substitutions.items():
                    monomial = monomial.subs(lhs, rhs)
                if original_monomial == monomial:
                    changed = False
                original_monomial = monomial
            return monomial

        length, h, U, t = 2, 3.8, -6, 1
        fu = generate_operators('fu', length)
        fd = generate_operators('fd', length)
        _b = flatten([fu, fd])
        hamiltonian = 0
        for j in range(length):
            hamiltonian += U * (Dagger(fu[j])*Dagger(fd[j]) * fd[j]*fu[j])
            hamiltonian += -h/2*(Dagger(fu[j])*fu[j] - Dagger(fd[j])*fd[j])
            for k in get_neighbors(j, len(fu), width=1):
                hamiltonian += -t*Dagger(fu[j])*fu[k]-t*Dagger(fu[k])*fu[j]
                hamiltonian += -t*Dagger(fd[j])*fd[k]-t*Dagger(fd[k])*fd[j]
        substitutions = fermionic_constraints(_b)
        monomials = expand(hamiltonian).as_coeff_mul()[1][0].as_coeff_add()[1]
        substituted_hamiltonian = sum([apply_substitutions(monomial,
                                                           substitutions)
                                       for monomial in monomials])
        correct_hamiltonian = sum([apply_correct_substitutions(monomial,
                                                               substitutions)
                                   for monomial in monomials])
        self.assertTrue(substituted_hamiltonian == expand(correct_hamiltonian))