Python sympy.solve函数代码示例

def SequentialSolving(eqns):
    sorteq = SortEquations(eqns)
    seqsol = {}
    unkn = sorted(sorteq[0].atoms(sy.Symbol)) 
    num_unkn = len(unkn)
    while num_unkn == 1 and not sorteq == []:  
        val_unkn = sy.solve(sorteq[0], unkn[0], simplify = False) 
        seqsol[unkn[0]] = sy.Rational(str(round(val_unkn[0], 4)))
    unkn = sorted(sorteq[0].atoms(sy.Symbol))    
    num_unkn = len(unkn)
    while num_unkn == 1:    
        val_unkn = sy.solve(sorteq[0], unkn[0])   
        seqsol[unkn[0]] = val_unkn[0]
        i = 0
        for eq in sorteq:
            if unkn[0] in sorteq[i]:    
                val_unkn[0] = round(val_unkn[0], 5)    
                repl = eq.subs(unkn[0], sy.Rational(str(val_unkn[0])))
                if not isinstance(repl, sy.Float):
                    sorteq[i] = repl
                else:
                    sorteq[i] = repl
            i+=1 
        
        sorteq = RemoveZeros(sorteq)
        if sorteq == []:
            num_unkn = 0
        else:
            sorteq = SortEquations(sorteq)
            unkn = sorted(sorteq[0].atoms(sy.Symbol))
            num_unkn = len(unkn)
                
    return seqsol

def benefit_from_demand(x, p, demand):
    """Converts the demand curve to the benefit.  It assumes that the
    demand is a x=d(p) function, where the x= is implicit.

    >>> sp.var('x p')
    (x, p)
    >>> sp.simplify(benefit_from_demand(x, p, 10/p -1) -
    ...             10*sp.log(x+1))
    0
    >>> sp.simplify(benefit_from_demand(x, p, sp.Eq(x, 10/p -1)) -
    ...             10*sp.log(x+1))
    0
    >>> sp.simplify(benefit_from_demand(x, p, 100-p) -
    ...             (-x**2/2 + 100*x))
    0
    >>> benefit_from_demand(x, p, sp.Piecewise((0, p < 0),
    ...                                        (-p + 100, p <= 100),
    ...                                        (0, True)))
    -x**2/2 + 100*x
    """
    if isinstance(demand, sp.relational.Relational):
        return sp.integrate(sp.solve(demand, p)[0], (x, 0, x))
    substracting = sp.solve(demand-x, p)
    if substracting:
        toint = substracting[0]
    else:
        substracting = sp.solve(demand, p)
        if substracting:
            toint = substracting[0] - x
        else:
            return None

    return sp.integrate(toint, (x, 0, x))

def getfunc(p1, p2, p3, p4):
    """ Get a point-returning function for a cubic
        curve over four points, with domain [0 - 3].
    """
    # knowns
    points = p1, p2, p3, p4

    # unknowns
    a, b, c, d, e, f, g, h = [Symbol(n) for n in 'abcdefgh']

    # coefficients
    xco = solve([(a * i**3 + b * i**2 + c * i + d - p.x) for (i, p) in enumerate(points)], [a, b, c, d])
    yco = solve([(e * i**3 + f * i**2 + g * i + h - p.y) for (i, p) in enumerate(points)], [e, f, g, h])

    # shorter variable names
    a, b, c, d = [xco[n] for n in (a, b, c, d)]
    e, f, g, h = [yco[n] for n in (e, f, g, h)]

    def func(t, d1=False):
        """ Return a position for given t or velocity (1st derivative) if arg. is True.
        """
        if d1:
            # first derivative
            return Point(3 * a * t**2 + 2 * b * t + c,
                         3 * e * t**2 + 2 * f * t + g)

        else:
            # actual function
            return Point(a * t**3 + b * t**2 + c * t + d,
                         e * t**3 + f * t**2 + g * t + h)

    return func

def print_assignment(eq, s, s0=0, pochoir=False):
	s1 = print_myccode(s, None, pochoir)
	if(s0==0):
		s2 = print_myccode(simplify(solve(eq,s)[0]), None, pochoir)
	else:
		s2 = print_myccode(simplify(solve(eq,s)[0] - s0) + s0, None, pochoir)
	return s1 + '=' + s2

 def __init__(self):
     h, g, h0, g0, a, d = sympy.symbols('h g h0 g0 a d')
     #g1 = sympy.Max(-d, g0 + h0 - 1/(4*a))
     g1 = g0 + h0 - 1/(4*a)
     h1 = h0 - 1/(2*a)
     parabola = d - a*h*h       # =g on boundary
     slope = g1 + h - h1        # =g on line with slope=1 through z1
     h2 = sympy.solve(parabola-slope,h)[0] # First solution is above z1
     g2 = h2 - h1 + g1          # Line has slope of 1
     g3 = g1
     h3 = sympy.solve((parabola-g).subs(g, g1),h)[1] # Second is on right
     r_a = sympy.Rational(1,24) # a=1/24 always
     self.h = tuple(x.subs(a,r_a) for x in (h0, h1, h2, h3))
     self.g = tuple(x.subs(a,r_a) for x in (g0, g1, g2, g3))
     def integrate(f):
         ia = sympy.integrate(
             sympy.integrate(f,(g, g1, g1+h-h1)),
             (h, h1, h2)) # Integral of f over right triangle
         ib = sympy.integrate(
             sympy.integrate(f,(g, g1, parabola)),
             (h, h2, h3)) # Integral of f over region against parabola
         return (ia+ib).subs(a, r_a)
     i0 = integrate(1)               # Area = integral of pie slice
     E = lambda f:(integrate(f)/i0)  # Expected value wrt Lebesgue measure
     sigma = lambda f,g:E(f*g) - E(f)*E(g)
     self.d = d
     self.Eh = collect(E(h),sympy.sqrt(d-g0-h0+6))
     self.Eg = E(g)
     self.Sigmahh = sigma(h,h)
     self.Sigmahg = sigma(h,g)
     self.Sigmagg = sigma(g,g)
     return

def solve_eq(rho_m, u_m, u_s):
    u_max, u_star, rho_max, rho_star, A, B = sympy.symbols("u_max u_star rho_max rho_star A B")

    eq1 = sympy.Eq(0, u_max * rho_max * (1 - A * rho_max - B * rho_max ** 2))
    eq2 = sympy.Eq(0, u_max * (1 - 2 * A * rho_star - 3 * B * rho_star ** 2))
    eq3 = sympy.Eq(u_star, u_max * (1 - A * rho_star - B * rho_star ** 2))

    eq4 = sympy.Eq(eq2.lhs - 3 * eq3.lhs, eq2.rhs - 3 * eq3.rhs)
    eq4.simplify()
    eq4.expand()

    rho_sol = sympy.solve(eq4, rho_star)[0]
    B_sol = sympy.solve(eq1, B)[0]

    quadA = eq2.subs([(rho_star, rho_sol), (B, B_sol)])
    quadA.simplify()

    A_sol = sympy.solve(quadA, A)[0]

    aval = A_sol.evalf(subs={u_star: u_s, u_max: u_m, rho_max: rho_m})
    bval = B_sol.evalf(subs={rho_max: rho_m, A: aval})

    rho_sol = sympy.solve(eq2, rho_star)[0]
    rho_val = rho_sol.evalf(subs={u_max: u_m, A: aval, B: bval})

    return aval, bval, rho_val

def PlotFor2Param(k1, k2, x, y, yVal, eq1, eq2, eqDet, eqTr, valK1m, valK3m, valK2 = 0.95, valK3 = 0.0032):
    eqk1_1Det = solve(eqDet.subs(x, eq2), k1)
    eqForK1Det = eq1 - eqk1_1Det[0]
    eqK2Det = solve(eqForK1Det.subs(km1, valK1m).subs(km3, valK3m).subs(k3, valK3), k2)
    funcK2Det = lambdify(y, eqK2Det[0])
    funcK1Det = lambdify(y, eqk1_1Det[0].subs(x, eq2).subs(k2, eqK2Det[0]).subs(km1, valK1m).subs(km3, valK3m).subs(k3, valK3))
    k1DetAll = funcK1Det(yVal)
    k2DetAll = funcK2Det(yVal)
    k1DetPos = [k1DetAll[i] for i in range(len(k1DetAll)) if k1DetAll[i] > 0 and k2DetAll[i] > 0]
    k2DetPos = [k2DetAll[i] for i in range(len(k2DetAll)) if k1DetAll[i] > 0 and k2DetAll[i] > 0]
    hopf, = plt.plot(k1DetPos, k2DetPos, linestyle = '--', color = 'r', label = 'hopf line')

    eqk1_1Tr = solve(eqTr.subs(x, eq2), k1)
    eqForK1Tr = eq1 - eqk1_1Tr[0]
    eqK2Tr = solve(eqForK1Tr.subs(km1, valK1m).subs(km3, valK3m).subs(k3, valK3), k2)
    funcK2Tr = lambdify(y, eqK2Tr[0])
    funcK1Tr = lambdify(y, eq1.subs(x, eq2).subs(k2, eqK2Tr[0]).subs(km1, valK1m).subs(km3, valK3m).subs(k3, valK3))
    k1AllTr = funcK1Tr(yVal)
    k2AllTr = funcK2Tr(yVal)
    k1PosTr = [k1AllTr[i] for i in range(len(k1AllTr)) if k1AllTr[i] > 0 and k2AllTr[i] > 0]
    print(len(k1PosTr))
    k2PosTr = [k2AllTr[i] for i in range(len(k2AllTr)) if k1AllTr[i] > 0 and k2AllTr[i] > 0]
    print(len(k2PosTr))
    sadle, = plt.plot(k1PosTr, k2PosTr, color = 'g', label = 'sadle-nodle line')

    plt.xlim([0, 2])
    plt.ylim([0,5])
    plt.legend(handles=[hopf, sadle])
    plt.xlabel('k1')
    plt.ylabel('k2')
    plt.show()

def test_issue_1572_1364_1368():
    assert solve((sqrt(x**2 - 1) - 2)) in ([sqrt(5), -sqrt(5)],
                                           [-sqrt(5), sqrt(5)])
    assert set(solve((2**exp(y**2/x) + 2)/(x**2 + 15), y)) == set([
        -sqrt(x)*sqrt(-log(log(2)) + log(log(2) + I*pi)),
        sqrt(x)*sqrt(-log(log(2)) + log(log(2) + I*pi))])

    C1, C2 = symbols('C1 C2')
    f = Function('f')
    assert solve(C1 + C2/x**2 - exp(-f(x)), f(x)) == [log(x**2/(C1*x**2 + C2))]
    a = Symbol('a')
    E = S.Exp1
    assert solve(1 - log(a + 4*x**2), x) in (
        [-sqrt(-a + E)/2, sqrt(-a + E)/2],
        [sqrt(-a + E)/2, -sqrt(-a + E)/2]
    )
    assert solve(log(a**(-3) - x**2)/a, x) in (
        [-sqrt(-1 + a**(-3)), sqrt(-1 + a**(-3))],
        [sqrt(-1 + a**(-3)), -sqrt(-1 + a**(-3))],)
    assert solve(1 - log(a + 4*x**2), x) in (
        [-sqrt(-a + E)/2, sqrt(-a + E)/2],
        [sqrt(-a + E)/2, -sqrt(-a + E)/2],)
    assert set(solve((
        a**2 + 1) * (sin(a*x) + cos(a*x)), x)) == set([-pi/(4*a), 3*pi/(4*a)])
    assert solve(3 - (sinh(a*x) + cosh(a*x)), x) == [2*atanh(S.Half)/a]
    assert set(solve(3 - (sinh(a*x) + cosh(a*x)**2), x)) == \
        set([
            2*atanh(-1 + sqrt(2))/a,
            2*atanh(S(1)/2 + sqrt(5)/2)/a,
            2*atanh(-sqrt(2) - 1)/a,
            2*atanh(-sqrt(5)/2 + S(1)/2)/a
             ])
    assert solve(atan(x) - 1) == [tan(1)]

    def generate_x(self, N=1000):
        '''
        this sampling method works wit all margins,
        but only frank and independent copulas can be used
        and it is limited to 2 dimensions
        '''
        # compute marginal prob of u1
        P_U1 = sy.simplify(sy.diff(self.C,self.U[0]))

        # invert marginal prob of u1
        y = sy.symbols('y')
        tmp = sy.solve(sy.Eq(P_U1,y),self.U[1])
        inv_P_U1 = sy.lambdify((self.U[0],y,self.D),tmp,'numpy')

        # invert margins
        inv_F = {}
        for m in [0,1]:
            u = sy.symbols('u')
            tmp = sy.solve(sy.Eq(self.F[m],u),self.X[m])[0]
            inv_F[m] = sy.lambdify((u,self.P[m]),tmp,'numpy')


        X = np.zeros((N,2))
        for i in range(N):
            u0, y = np.random.uniform(size=2)
            u1 = inv_P_U1(u0,y,self.C_para)
            X[i,0] = inv_F[0](u0,self.F_para[0])
            X[i,1] = inv_F[1](u1,self.F_para[1])

        return X

def test_minsolve_linear_system():
    def count(dic):
        return len([x for x in dic.itervalues() if x == 0])
    assert count(solve([x + y + z, y + z + a + t], minimal=True, quick=True)) == 3
    assert count(solve([x + y + z, y + z + a + t], minimal=True, quick=False)) == 3
    assert count(solve([x + y + z, y + z + a], minimal=True, quick=True)) == 1
    assert count(solve([x + y + z, y + z + a], minimal=True, quick=False)) == 2

def test_issue_2813():
    assert solve(x ** 2 - x - 0.1, rational=True) == [S(1) / 2 + sqrt(35) / 10, -sqrt(35) / 10 + S(1) / 2]
    # [-0.0916079783099616, 1.09160797830996]
    ans = solve(x ** 2 - x - 0.1, rational=False)
    assert len(ans) == 2 and all(a.is_Number for a in ans)
    ans = solve(x ** 2 - x - 0.1)
    assert len(ans) == 2 and all(a.is_Number for a in ans)

def test_CRootOf___eval_Eq__():
    f = Function('f')
    eq = x**3 + x + 3
    r = rootof(eq, 2)
    r1 = rootof(eq, 1)
    assert Eq(r, r1) is S.false
    assert Eq(r, r) is S.true
    assert Eq(r, x) is S.false
    assert Eq(r, 0) is S.false
    assert Eq(r, S.Infinity) is S.false
    assert Eq(r, I) is S.false
    assert Eq(r, f(0)) is S.false
    assert Eq(r, f(0)) is S.false
    sol = solve(eq)
    for s in sol:
        if s.is_real:
            assert Eq(r, s) is S.false
    r = rootof(eq, 0)
    for s in sol:
        if s.is_real:
            assert Eq(r, s) is S.true
    eq = x**3 + x + 1
    sol = solve(eq)
    assert [Eq(rootof(eq, i), j) for i in range(3) for j in sol] == [
        False, False, True, False, True, False, True, False, False]
    assert Eq(rootof(eq, 0), 1 + S.ImaginaryUnit) == False

def test_solve_inequalities():
    system = [Lt(x ** 2 - 2, 0), Gt(x ** 2 - 1, 0)]

    assert solve(system) == And(
        Or(And(Lt(-sqrt(2), re(x)), Lt(re(x), -1)), And(Lt(1, re(x)), Lt(re(x), sqrt(2)))), Eq(im(x), 0)
    )
    assert solve(system, assume=Q.real(x)) == Or(And(Lt(-sqrt(2), x), Lt(x, -1)), And(Lt(1, x), Lt(x, sqrt(2))))

def test_RootOf___eval_Eq__():
    f = Function('f')
    r = RootOf(x**3 + x + 3, 2)
    r1 = RootOf(x**3 + x + 3, 1)
    assert Eq(r, r1) is S.false
    assert Eq(r, r) is S.true
    assert Eq(r, x) is S.false
    assert Eq(r, 0) is S.false
    assert Eq(r, S.Infinity) is S.false
    assert Eq(r, I) is S.false
    assert Eq(r, f(0)) is S.false
    assert Eq(r, f(0)) is S.false
    sol = solve(r.expr)
    for s in sol:
        if s.is_real:
            assert Eq(r, s) is S.false
    r = RootOf(r.expr, 0)
    for s in sol:
        if s.is_real:
            assert Eq(r, s) is S.true
    eq = (x**3 + x + 1)
    assert [Eq(RootOf(eq, i), j) for i in range(3) for j in solve(eq)] == [
        False, False, True, False, True, False, True, False, False
    ]
    assert Eq(RootOf(eq, 0), 1 + S.ImaginaryUnit) == False

    def __init__(self, width, height, x_radius=None,
                 curve=CIRCULAR):
        self.width = width
        self.height = height
        self.x_radius = x_radius
        self.curve = curve

        self.x, self.y, dx, dy = sympy.symbols("x y dx dy", real=True)
        a, b = sympy.symbols("a b", positive=True, real=True)

        self.ellipse = ((self.x - dx) / a) ** 2 + ((self.y - dy) / b) ** 2 - 1

        self.ellipse = self.ellipse.subs([(dx, self.width),
                                          (dy, b)])

        if curve == Transition.CIRCULAR:
            self.ellipse = self.ellipse.subs([(a, b)])

        ellipse = self.ellipse.subs([(self.x, 0),
                                     (self.y, self.height)])

        if curve == Transition.CIRCULAR:
            self.x_radius = self.y_radius = max(sympy.solve(ellipse, b))
            self.angle = math.asin(self.width / self.x_radius)
            self.arc_length = self.x_radius * self.angle

        elif curve == Transition.ELLIPTICAL:
            self.y_radius = max(sympy.solve(ellipse, b))
            # TODO: Fix me
            self.angle = 0
            self.arc_length = 0


        self.ellipse = self.ellipse.subs([(a, self.x_radius),
                                          (b, self.y_radius)])