Python sympy.roots函数代码示例

    def ZPK(self):
        """Convert to pole-zero-gain (PZK) form.

        See also canonical, general, mixedfrac, and partfrac"""

        N, D, delay = self._as_ratfun_delay()

        K = sym.cancel(N.LC() / D.LC())
        if delay != 0:
            K *= sym.exp(self.var * delay)

        zeros = sym.roots(N)
        poles = sym.roots(D)

        return self.__class__(_zp2tf(zeros, poles, K, self.var))

def get_world_segments(root_a, root_b, root_c, initial_t, final_t, intersection_radius):
    """
    The world consists of
    three axis lines,
    six circles marking the intersections,
    and a single parametric curve.
    @return: a collection of (p0, p1, style) triples
    """
    seg_length_min = 0.1
    segments = []
    # add the axis line segments
    d = 5
    f_x = pcurve.LineSegment(np.array([-d, 0, 0]), np.array([d, 0, 0]))
    f_y = pcurve.LineSegment(np.array([0, -d, 0]), np.array([0, d, 0]))
    f_z = pcurve.LineSegment(np.array([0, 0, -d]), np.array([0, 0, d]))
    x_axis_segs = pcurve.get_piecewise_curve(f_x, 0, 1, 10, seg_length_min)
    y_axis_segs = pcurve.get_piecewise_curve(f_y, 0, 1, 10, seg_length_min)
    z_axis_segs = pcurve.get_piecewise_curve(f_z, 0, 1, 10, seg_length_min)
    segments.extend((p0, p1, STYLE_X) for p0, p1 in x_axis_segs)
    segments.extend((p0, p1, STYLE_Y) for p0, p1 in y_axis_segs)
    segments.extend((p0, p1, STYLE_Z) for p0, p1 in z_axis_segs)
    # add the parametric curve
    roots = (root_a, root_b, root_c)
    polys = interlace.roots_to_differential_polys(roots)
    f_poly = interlace.Multiplex((sympyutils.WrappedUniPoly(p) for p in polys))
    poly_segs = pcurve.get_piecewise_curve(f_poly, initial_t, final_t, 10, seg_length_min)
    segments.extend((p0, p1, STYLE_CURVE) for p0, p1 in poly_segs)
    # add the intersection circles
    x_roots_symbolic = sympy.roots(polys[0])
    y_roots_symbolic = sympy.roots(polys[1])
    z_roots_symbolic = sympy.roots(polys[2])
    x_roots = [float(r) for r in x_roots_symbolic]
    y_roots = [float(r) for r in y_roots_symbolic]
    z_roots = [float(r) for r in z_roots_symbolic]
    for r in x_roots:
        f = pcurve.OrthoCircle(f_poly(r), intersection_radius, 0)
        segs = pcurve.get_piecewise_curve(f, 0, 1, 10, seg_length_min)
        segments.extend((p0, p1, STYLE_X) for p0, p1 in segs)
    for r in y_roots:
        f = pcurve.OrthoCircle(f_poly(r), intersection_radius, 1)
        segs = pcurve.get_piecewise_curve(f, 0, 1, 10, seg_length_min)
        segments.extend((p0, p1, STYLE_Y) for p0, p1 in segs)
    for r in z_roots:
        f = pcurve.OrthoCircle(f_poly(r), intersection_radius, 2)
        segs = pcurve.get_piecewise_curve(f, 0, 1, 10, seg_length_min)
        segments.extend((p0, p1, STYLE_Z) for p0, p1 in segs)
    # return the segments
    return segments

def test_legendre():
    raises(ValueError, lambda: legendre(-1, x))
    assert legendre(0, x) == 1
    assert legendre(1, x) == x
    assert legendre(2, x) == ((3*x**2 - 1)/2).expand()
    assert legendre(3, x) == ((5*x**3 - 3*x)/2).expand()
    assert legendre(4, x) == ((35*x**4 - 30*x**2 + 3)/8).expand()
    assert legendre(5, x) == ((63*x**5 - 70*x**3 + 15*x)/8).expand()
    assert legendre(6, x) == ((231*x**6 - 315*x**4 + 105*x**2 - 5)/16).expand()

    assert legendre(10, -1) == 1
    assert legendre(11, -1) == -1
    assert legendre(10, 1) == 1
    assert legendre(11, 1) == 1
    assert legendre(10, 0) != 0
    assert legendre(11, 0) == 0

    assert roots(legendre(4, x), x) == {
        sqrt(Rational(3, 7) - Rational(2, 35)*sqrt(30)): 1,
        -sqrt(Rational(3, 7) - Rational(2, 35)*sqrt(30)): 1,
        sqrt(Rational(3, 7) + Rational(2, 35)*sqrt(30)): 1,
        -sqrt(Rational(3, 7) + Rational(2, 35)*sqrt(30)): 1,
    }

    n = Symbol("n")

    X = legendre(n, x)
    assert isinstance(X, legendre)

    assert legendre(-n, x) == legendre(n - 1, x)
    assert legendre(n, -x) == (-1)**n*legendre(n, x)
    assert diff(legendre(n, x), x) == \
        n*(x*legendre(n, x) - legendre(n - 1, x))/(x**2 - 1)
    assert diff(legendre(n, x), n) == Derivative(legendre(n, x), n)

    def test_singularities(self):
        S1 = sympify(singularities(f1,x,y))
        S2 = sympify(singularities(f2,x,y))
        S3 = sympify(singularities(f3,x,y))
        S4 = sympify(singularities(f4,x,y))
        S5 = sympify(singularities(f5,x,y))
        S6 = sympify(singularities(f6,x,y))
#         S7 = sympify(singularities(f7,x,y))
#         S8 = sympify(singularities(f8,x,y))
        S9 = sympify(singularities(f9,x,y))
        S10= sympify(singularities(f10,x,y))

        rt5 = sympy.roots(x**2 + 1, x).keys()

        S1act = sympify([
            ((0,0,1),(2,2,1)),
            ((0,1,0),(2,1,2))
            ])
        S2act = sympify([
            ((0,0,1),(3,4,2)),
            ((0,1,0),(4,9,1))
            ])
        S3act = sympify([
            ((0,0,1),(2,1,2)),
            ((1,-1,1),(2,1,2)),
            ((1,1,1),(2,1,2))
            ])
        S4act = sympify([
            ((0,0,1),(2,1,2))
            ])
        S5act = sympify([
            ((0,0,1),(3,3,3)),
            ((1,rt5[1],0),(3,3,3)),
            ((1,rt5[0],0),(3,3,3))
            ])
        S6act = sympify([
            ((0,0,1),(2,2,2)),
            ((1,0,0),(2,2,2))
            ])
#         S7act = sympify([((0,1,0),(3,6,3))])
#         S8act = sympify([
#             ((0,1,0),(6,21,3)),
#             ((1,0,0),(3,7,2))
#             ])
        S9act = sympify([((0,1,0),(5,12,1))])
        S10act= sympify([
            ((0,1,0),(3,6,1)),
            ((1,0,0),(4,6,4))
            ])

        self.assertItemsEqual(S1,S1act)
        self.assertItemsEqual(S2,S2act)
        self.assertItemsEqual(S3,S3act)
        self.assertItemsEqual(S4,S4act)
        self.assertItemsEqual(S5,S5act)
        self.assertItemsEqual(S6,S6act)
#         self.assertItemsEqual(S7,S7act)
#         self.assertItemsEqual(S8,S8act)
        self.assertItemsEqual(S9,S9act)
        self.assertItemsEqual(S10,S10act)

    def test_singular_points(self):
        S1 = singularities(f1,x,y)
        S2 = singularities(f2,x,y)
        S3 = singularities(f3,x,y)
        S4 = singularities(f4,x,y)
        S5 = singularities(f5,x,y)
        S6 = singularities(f6,x,y)
        S7 = singularities(f7,x,y)
        S8 = singularities(f8,x,y)
        S9 = singularities(f9,x,y)
        S10= singularities(f10,x,y)

        rt5 = [rt for rt,_ in sympy.roots(x**2 + 1, x).iteritems()]

        S1act = [(0,0,1),(0,1,0)]
        S2act = [(0,0,1),(0,1,0)]
        S3act = [(0,0,1),(1,-1,1),(1,1,1)]
        S4act = [(0,0,1)]
        S5act = [(0,0,1),(rt5[0],1,0),(rt5[1],1,0)]
        S6act = [(0,0,1),(1,0,0)]
        S7act = [(0,1,0)]
        S8act = [(0,1,0),(1,0,0)]
        S9act = [(0,1,0)]
        S10act= [(0,1,0),(1,0,0)]

        self.assertItemsEqual(S1,S1act)
        self.assertItemsEqual(S2,S2act)
        self.assertItemsEqual(S3,S3act)
        self.assertItemsEqual(S4,S4act)
        self.assertItemsEqual(S5,S5act)
        self.assertItemsEqual(S6,S6act)
        self.assertItemsEqual(S7,S7act)
        self.assertItemsEqual(S8,S8act)
        self.assertItemsEqual(S9,S9act)
        self.assertItemsEqual(S10,S10act)

    def ZPK(self):
        """Convert to pole-zero-gain (PZK) form.

        See also canonical, general, mixedfrac, and partfrac"""

        N, D, delay, undef = self.as_ratfun_delay_undef()

        var = self.var        
        Npoly = sym.Poly(N, var)
        Dpoly = sym.Poly(D, var)
        
        K = sym.cancel(Npoly.LC() / Dpoly.LC())
        if delay != 0:
            K *= sym.exp(self.var * delay)

        zeros = sym.roots(Npoly)
        poles = sym.roots(Dpoly)

        return _zp2tf(zeros, poles, K, self.var) * undef

 def test_roots_to_poly(self):
     roots = (1.0, 4.0, 5.0)
     p = sympyutils.roots_to_poly(roots)
     self.assertTrue(p.is_monic)
     root_to_count = sympy.roots(p)
     self.assertEqual(set(root_to_count.values()), set([1]))
     observed = set(root_to_count.keys())
     expected = set(roots)
     for r_observed, r_expected in zip(sorted(observed), sorted(expected)):
         self.assertAlmostEqual(r_observed, r_expected)

def run_copper_smith(f, N, h, k):
    X = s.S(s.ceiling((s.Pow(2, - 1 / 2) * s.Pow(h * k, -1 / (h * k - 1))) *
                      s.Pow(N, (h - 1) / (h * k - 1))) - 1)
    gen = lll(generate(f, N, h, k, X), s.S(0.75))
    final_poly = [gen[0][i] / s.Pow(X, i) for i in range(len(gen[0]))][::-1]

    roots = s.roots(final_poly)
    for r in roots:
        print(r)

    return roots

def singular_term(F,X,Y,L,I,version):
    """
    Computes a single set of singular terms of the Puiseux expansions.

    For `I=1`, the function computes the first term of each finite 
    `\mathbb{K}` expansion. For `I=2`, it computes the other terms of the
    expansion.
    """
    T = []

    # if the curve is singular then compute the singular tuples
    # otherwise, use the standard newton polygon method
    if is_singular(F,X,Y):
        for (q,m,l,Phi) in desingularize(F,X,Y):
            for eta in Phi.all_roots(radicals=True):
                tau = (q,1,m,1,sympy.together(eta))
                T.append((tau,0,1))
    else:
        # each side of the newton polygon corresponds to a K-term.
        for (q,m,l,Phi) in polygon(F,X,Y,I):
            # the rational method
            if version == 'rational':
                u,v = _bezout(q,m)      

            # each newton polygon side has a characteristic
            # polynomial. For each square-free factor, each root
            # corresponds to a K-term
            for (Psi,r) in _square_free(Phi,_Z):
                Psi = sympy.Poly(Psi,_Z)
                
                # compute the roots of Psi. Use the RootOf construct if 
                # possible. In the case when Psi is over EX (i.e. when
                # RootOf doesn't work) then compute symbolic roots.
                try:
                    roots = Psi.all_roots(radicals=True)
                except NotImplementedError:
                    roots = sympy.roots(Psi,_Z).keys()

                for xi in roots:
                    # the classical version returns the "raw" roots
                    if version == 'classical':
                        P = sympy.Poly(_U**q-xi,_U)
                        beta = RootOf(P,0,radicals=True)
                        tau = (q,1,m,beta,1)
                        T.append((tau,l,r))
                    # the rational version rescales parameters so as
                    # to include ony rational terms in the Puiseux
                    # expansions.
                    if version == 'rational':
                        mu = xi**(-v)
                        beta = sympy.together(xi**u)
                        tau = (q,mu,m,beta,1)
                        T.append((tau,l,r))
    return T

def _singular_points_finite(f,x,y):
    """
    Returns the finite singular points of f.
    """
    S = []

    # compute the finite singularities: use the resultant
    p  = sympy.Poly(f,[x,y])
    n  = p.degree(y)

    res = sympy.Poly(sympy.resultant(p,p.diff(y),y),x)
    for xk,deg in sympy.roots(res,x).iteritems():
        if deg > 1:
            fxk = sympy.Poly(f.subs({x:xk,y:_Z}), _Z)
            for ykj,_ in sympy.roots(fxk, _Z).iteritems():
                fx = f.diff(x)
                fy = f.diff(y)
                subs = {x:xk,y:ykj}
                if (fx.subs(subs) == 0) and (fy.subs(subs) == 0):
                    S.append((xk,ykj,1))

    return S

def fund_ries(drc):
    a2,a1,a0 = drc
    r = roots(a2*t**2 + a1*t +a0,t)
    kl = list(r.keys())
    if len(kl) == 2:  # jednoduche korene
        r1,r2 = kl        
        re_k,im_k = r1.as_real_imag()
        if im_k == 0: # jednoduche realne korene
            return exp(r1*t), exp(r2*t)
        else:         # komplexne zdruz.korene
            return exp(re_k*t)*sin(im_k*t), exp(re_k*t)*cos(im_k*t)
    else: # dvojnasobny koren
        return exp(kl[0]*t), t*exp(kl[0]*t)

def _singular_points_infinite(f,x,y):
    """
    Returns the singular points of f at infinity.
    """
    S = []

    # compute homogenous polynomial
    F, d = homogenize(f,x,y,_z)

    # find singular points at infinity
    F0 = sympy.Poly(F.subs([(_z,0)]),[x,y])
    domain = sympy.QQ[sympy.I]
    solsX1 = sympy.roots(F0.subs({x:1,y:_Z}),_Z).keys()
    solsY1 = sympy.roots(F0.subs({y:1,x:_Z}),_Z).keys()

    all_sols = [ (1,yi,0) for yi in solsX1 ]
    all_sols.extend( [ (xi,1,0) for xi in solsY1 ] )

    # these points are in projective space, so filter out equal points
    # such as (1,I,0) == (-I,1,0). We normalize these projective
    # points such that 1 appears in either the x- or y- coordinate
    sols = []
    for xi,yi,zi in all_sols:
        normalized = (1,yi/xi,0) if xi != sympy.S(0) else (xi/yi,1,0)
        if not normalized in sols:
            sols.append(normalized)

    # Filter out any points that are not singular by checking if the
    # gradient vanishes.
    grad = [F.diff(var) for var in [x,y,_z]]
    for xi,yi,zi in sols:
        fsub = lambda e,x=x,y=y,_z=_z,xi=xi,yi=yi,zi=zi:  \
               e.subs({x:xi,y:yi,_z:zi}).simplify() != sympy.S(0)
        if not any(map(fsub,grad)):
            S.append((xi,yi,zi))

    return S

def tustin(tf, sRate):
	s =	getMapping(tf)['s']
	T = 1/sRate
	poles = sympy.roots(sympy.denom(tf), s, multiple=True)
	centroid = np.mean(np.abs(poles))
	invz = sympy.Symbol('invz', real=True)
	# normalized center frequency derived from poles of filter SISO transfer function
	w0 = 2*np.pi*centroid/(sRate*2*np.pi)
	# modified bilinear transform w/ frequency warping
	bt = w0/sympy.tan(w0*T/2) * ((1-invz)/(1+invz))
	dt = sympy.simplify(tf.subs({s: bt}))
	b = sympy.Poly(sympy.numer(dt)).all_coeffs()[::-1]
	a = sympy.Poly(sympy.denom(dt)).all_coeffs()[::-1]
	normalize = lambda x: float(x/a[0])
	return (map(normalize, b), map(normalize, a))

 def test_ndiff_roots_symbolic(self):
     roots = (1.25, 4.5, 5.75)
     a, b, c = roots
     polys = roots_to_differential_polys(roots)
     # compute linear root manually
     r = float(a + b + c) / 3
     # check linear root
     observed = sorted(float(r) for r in sympy.roots(polys[0]))
     expected = sorted([r])
     self.assertTrue(np.allclose(observed, expected))
     # compute quadratic roots manually
     A = a*a + b*b + c*c
     B = a*b + a*c + b*c
     S = a + b + c
     r0 = float(S + math.sqrt(A - B)) / 3
     r1 = float(S - math.sqrt(A - B)) / 3
     # check quadratic roots
     observed = sorted(float(r) for r in sympy.roots(polys[1]))
     expected = sorted([r0, r1])
     self.assertTrue(np.allclose(observed, expected))
     # check cubic roots
     observed = sorted(float(r) for r in sympy.roots(polys[2]))
     expected = sorted(roots)
     self.assertTrue(np.allclose(observed, expected))

def get_segmentation(p, t0, t1):
    """
    A segmentation is a sequence of triples (left, right, sign).
    @param p: a sympy Poly
    @param t0: initial time
    @param t1: final time
    @return: a segmentation
    """
    roots = sorted(float(r) for r in sympy.roots(p))
    points = [t0] + roots + [t1]
    segmentation = []
    for left, right in iterutils.pairwise(points):
        mid = (left + right) / 2
        sign = -1 if p.eval(mid) <= 0 else 1
        seg = (left, right, sign)
        segmentation.append(seg)
    return segmentation