Python sympy.legendre函数代码示例

def decoding3dMaxRe(l) :
	if l is 0 : return numpy.array([1])
	import sympy
	x=sympy.Symbol("x")
	maxRoot = max((float(root) for root in sympy.solve(sympy.legendre(l+1,x),x)))
	return numpy.array([1.]+ [float(sympy.legendre(m, maxRoot))
		for m in xrange(1,l+1)])

def test_assoc_legendre():
    Plm = assoc_legendre
    Q = sqrt(1 - x**2)

    assert Plm(0, 0, x) == 1
    assert Plm(1, 0, x) == x
    assert Plm(1, 1, x) == -Q
    assert Plm(2, 0, x) == (3*x**2 - 1)/2
    assert Plm(2, 1, x) == -3*x*Q
    assert Plm(2, 2, x) == 3*Q**2
    assert Plm(3, 0, x) == (5*x**3 - 3*x)/2
    assert Plm(3, 1, x).expand() == (( 3*(1 - 5*x**2)/2 ).expand() * Q).expand()
    assert Plm(3, 2, x) == 15*x * Q**2
    assert Plm(3, 3, x) == -15 * Q**3

    # negative m
    assert Plm(1, -1, x) == -Plm(1, 1, x)/2
    assert Plm(2, -2, x) == Plm(2, 2, x)/24
    assert Plm(2, -1, x) == -Plm(2, 1, x)/6
    assert Plm(3, -3, x) == -Plm(3, 3, x)/720
    assert Plm(3, -2, x) == Plm(3, 2, x)/120
    assert Plm(3, -1, x) == -Plm(3, 1, x)/12

    n = Symbol("n")
    m = Symbol("m")

    X = Plm(n, m, x)
    assert isinstance(X, assoc_legendre)

    assert Plm(n, 0, x) == legendre(n, x)

    raises(ValueError, lambda: Plm(-1, 0, x))
    raises(ValueError, lambda: Plm(0, 1, x))

def gauss_lobatto_points(p):
    """
    Returns the list of Gauss-Lobatto points of the order 'p'.
    """
    x = Symbol("x")
    print "creating"
    e = legendre(p, x).diff(x)
    e = Poly(e, x)
    print "polydone"
    if e == 0:
        return []
    print "roots"
    if e == 1:
        r = []
    else:
        with workdps(40):
            r, err = polyroots(e.all_coeffs(), error=True)
        #if err > 1e-40:
        #    raise Exception("Internal Error: Root is not precise")
    print "done"
    p = []
    p.append("-1.0")
    for x in r:
        if abs(x) < 1e-40: x = 0
        p.append(str(x))
    p.append("1.0")
    return p

def test_jacobi():
    n = Symbol("n")
    a = Symbol("a")
    b = Symbol("b")

    assert jacobi(0, a, b, x) == 1
    assert jacobi(1, a, b, x) == a/2 - b/2 + x*(a/2 + b/2 + 1)

    assert jacobi(n, a, a, x) == RisingFactorial(a + 1, n)*gegenbauer(n, a + S(1)/2, x)/RisingFactorial(2*a + 1, n)
    assert jacobi(n, a, -a, x) == ((-1)**a*(-x + 1)**(-a/2)*(x + 1)**(a/2)*assoc_legendre(n, a, x)*
                                   factorial(-a + n)*gamma(a + n + 1)/(factorial(a + n)*gamma(n + 1)))
    assert jacobi(n, -b, b, x) == ((-x + 1)**(b/2)*(x + 1)**(-b/2)*assoc_legendre(n, b, x)*
                                   gamma(-b + n + 1)/gamma(n + 1))
    assert jacobi(n, 0, 0, x) == legendre(n, x)
    assert jacobi(n, S.Half, S.Half, x) == RisingFactorial(S(3)/2, n)*chebyshevu(n, x)/factorial(n + 1)
    assert jacobi(n, -S.Half, -S.Half, x) == RisingFactorial(S(1)/2, n)*chebyshevt(n, x)/factorial(n)

    X = jacobi(n, a, b, x)
    assert isinstance(X, jacobi)

    assert jacobi(n, a, b, -x) == (-1)**n*jacobi(n, b, a, x)
    assert jacobi(n, a, b, 0) == 2**(-n)*gamma(a + n + 1)*hyper((-b - n, -n), (a + 1,), -1)/(factorial(n)*gamma(a + 1))
    assert jacobi(n, a, b, 1) == RisingFactorial(a + 1, n)/factorial(n)

    m = Symbol("m", positive=True)
    assert jacobi(m, a, b, oo) == oo*RisingFactorial(a + b + m + 1, m)

    assert conjugate(jacobi(m, a, b, x)) == jacobi(m, conjugate(a), conjugate(b), conjugate(x))

    assert diff(jacobi(n,a,b,x), n) == Derivative(jacobi(n, a, b, x), n)
    assert diff(jacobi(n,a,b,x), x) == (a/2 + b/2 + n/2 + S(1)/2)*jacobi(n - 1, a + 1, b + 1, x)

def legendre_norm(i, x):
    """
    Returns the normalized integrated Legendre polynomial.
    """
    f = legendre(i, x)
    n = sqrt(integrate(f**2, (x, -1, 1)))
    return f/n

    def polyrs(self, uhat):
        "Compute the polynomial for the response surface."
        import sympy

        dim = len(self.params)
        var = []
        for i, p in enumerate(self.params):
            var.append(sympy.Symbol(str(p.name)))

        chaos = chaos_sequence(dim, self.level)
        chaos = np.int_(chaos)

        for d in range(0, dim):
            poly = np.array(map(lambda x: sympy.legendre(x, var[d]), chaos[:, d]))
            if d == 0:
                s = poly
            else:
                s *= poly

        eqn = np.sum(uhat * s)
        for i, p in enumerate(self.params):
            pmin, pmax = p.pdf.srange
            c = (pmax + pmin) / 2.0
            s = (pmax - pmin) / 2.0
            eqn = eqn.subs(var[i], (var[i] - c)/s)
        return sympy.sstr(eqn.expand().evalf())

def test_gegenbauer():
    n = Symbol("n")
    a = Symbol("a")

    assert gegenbauer(0, a, x) == 1
    assert gegenbauer(1, a, x) == 2*a*x
    assert gegenbauer(2, a, x) == -a + x**2*(2*a**2 + 2*a)
    assert gegenbauer(3, a, x) == \
        x**3*(4*a**3/3 + 4*a**2 + 8*a/3) + x*(-2*a**2 - 2*a)

    assert gegenbauer(-1, a, x) == 0
    assert gegenbauer(n, S(1)/2, x) == legendre(n, x)
    assert gegenbauer(n, 1, x) == chebyshevu(n, x)
    assert gegenbauer(n, -1, x) == 0

    X = gegenbauer(n, a, x)
    assert isinstance(X, gegenbauer)

    assert gegenbauer(n, a, -x) == (-1)**n*gegenbauer(n, a, x)
    assert gegenbauer(n, a, 0) == 2**n*sqrt(pi) * \
        gamma(a + n/2)/(gamma(a)*gamma(-n/2 + S(1)/2)*gamma(n + 1))
    assert gegenbauer(n, a, 1) == gamma(2*a + n)/(gamma(2*a)*gamma(n + 1))

    assert gegenbauer(n, Rational(3, 4), -1) == zoo

    m = Symbol("m", positive=True)
    assert gegenbauer(m, a, oo) == oo*RisingFactorial(a, m)

    assert conjugate(gegenbauer(n, a, x)) == gegenbauer(n, conjugate(a), conjugate(x))

    assert diff(gegenbauer(n, a, x), n) == Derivative(gegenbauer(n, a, x), n)
    assert diff(gegenbauer(n, a, x), x) == 2*a*gegenbauer(n - 1, a + 1, x)

def test_manualintegrate_orthogonal_poly():
    n = symbols('n')
    a, b = 7, S(5)/3
    polys = [jacobi(n, a, b, x), gegenbauer(n, a, x), chebyshevt(n, x),
        chebyshevu(n, x), legendre(n, x), hermite(n, x), laguerre(n, x),
        assoc_laguerre(n, a, x)]
    for p in polys:
        integral = manualintegrate(p, x)
        for deg in [-2, -1, 0, 1, 3, 5, 8]:
            # some accept negative "degree", some do not
            try:
                p_subbed = p.subs(n, deg)
            except ValueError:
                continue
            assert (integral.subs(n, deg).diff(x) - p_subbed).expand() == 0

        # can also integrate simple expressions with these polynomials
        q = x*p.subs(x, 2*x + 1)
        integral = manualintegrate(q, x)
        for deg in [2, 4, 7]:
            assert (integral.subs(n, deg).diff(x) - q.subs(n, deg)).expand() == 0

        # cannot integrate with respect to any other parameter
        t = symbols('t')
        for i in range(len(p.args) - 1):
            new_args = list(p.args)
            new_args[i] = t
            assert isinstance(manualintegrate(p.func(*new_args), t), Integral)

def legendre_int(i, x):
    """
    Returns the normalized integrated Legendre polynomial.
    """
    y = Symbol("y", dummy=True)
    f = legendre(i, y)
    n = sqrt(integrate(f**2, (y, -1, 1)))
    return integrate(f, (y, -1, x))/n

def acceptance( cosTheta1, phi, cosTheta2, i_max, j_max, k_max, numEvents, c):
	returnValue = 0.
        for i in range(i_max+1):
		for k in range(3):
			for j in range(0, 5, 2): # must have l >= k
 				if (j < k): continue
 				P_i  = N(legendre(i,    cosTheta2))
				Y_jk = abs(N(Zlm(j, k, acos(cosTheta1), phi)))
				#print returnValue, P_i, Y_jk, c[i][j][k]
 				returnValue += c[i][k][j]*(P_i * Y_jk);
	return returnValue

def gauss_lobatto_jacobi_quadruature_points_weights(n=5):
    """
    compute Gauss-Lobatto-Jacobi (GLJ) quadrature weights and points [-1, 1]
    using the sympy symbolic package and analytical definitions of these points

    somehow seems not to work for even numbers n >= 10, I assume this is a bug
    in sympy

    :param n: number of integration points (order + 1)
    :type n: integer

    :returns: tuple of two numpy arrays of floats containing the points and
        weights
    """
    x = sp.symbols('x')

    # GLL points are defined as the roots of the derivative of a sum of
    # Legendre Polynomials and include the boundaries
    Pn_bar = (sp.legendre(n-1, x) + sp.legendre(n, x)) / (1 + x)
    dPn_bar = sp.cancel((1 - x ** 2) * sp.diff(Pn_bar, x))

    # evaluate and sort
    pointsl = []
    for i in np.arange(n):
        p = sp.RootOf(dPn_bar, i)
        pointsl.append(p.evalf())

    pointsl.sort()

    # weights
    weightsl = []
    for p in pointsl:
        if p == -1.:
            wi = 8. / (n ** 2 * (n + 1) * (n - 1))
        else:
            wi = 4. / ((n + 1) * (n - 1) * Pn_bar.subs(x, p) ** 2)
        weightsl.append(wi)

    return np.array(pointsl, dtype='float'), np.array(weightsl, dtype='float')

    def test_five_points(self):
        import numpy as np
        import sympy

        t_symb = sympy.Symbol('t')
        num_points = 5
        p4 = sympy.legendre(num_points - 1, t_symb)
        p4_prime = p4.diff(t_symb)
        start = -1.0
        stop = 1.0
        all_nodes = self._call_func_under_test(start, stop, num_points)
        self.assertEqual(all_nodes[0], start)
        self.assertEqual(all_nodes[-1], stop)
        inner_nodes = all_nodes[1:-1]
        # Make sure the computed roots are actually roots.
        self.assertTrue(np.allclose(
            [float(p4_prime.subs({t_symb: root}))
             for root in inner_nodes], 0))

def gauss_lobatto_points(p):
    """
    Returns the list of Gauss-Lobatto points of the order 'p'.
    """
    x = Symbol("x")
    e = (1-x**2)*legendre(p, x).diff(x)
    e = Poly(e, x)
    if e == 0:
        return []
    with workdps(40):
        r, err = polyroots(e.all_coeffs(), error=True)
    if err > 1e-40:
        raise Exception("Internal Error: Root is not precise")
    p = []
    for x in r:
        if abs(x) < 1e-40: x = 0
        p.append(str(x))
    return p

def gauss_quadruature_points_weights(n=5):
    """
    compute Gaussian quadrature weights and points [-1, 1] using the sympy
    symbolic package and analytical definitions of these points

    somehow seems not to work for even numbers n >= 10, I assume this is a bug
    in sympy

    :param n: number of integration points (order + 1)
    :type n: integer

    :returns: tuple of two numpy arrays of floats containing the points and
        weights
    """
    x = sp.symbols('x')

    # Gauss points are defined as the roots of the Legendre Polynomials
    Pn = sp.legendre(n, x)

    # evaluate and sort
    pointsl = []
    for i in np.arange(n):
        p = sp.RootOf(Pn, i)
        pointsl.append(p.evalf())

    pointsl.sort()

    # weights are found through the derivative of the Legendre Polynomials
    dPn = sp.diff(Pn)

    # evaluate
    weightsl = []
    for p in pointsl:
        wi = 2. / ((1 - p ** 2) * dPn.subs(x, p) ** 2)
        weightsl.append(wi)

    return np.array(pointsl, dtype='float'), np.array(weightsl, dtype='float')

    def __init__(self, n=None, verbose=False):
        self.n = n
        prec = dps = Float.getdps()
        self.prec = prec

        # Reuse old nodes
        if n in _gausscache:
            cacheddps, xs, ws = _gausscache[n]
            if cacheddps >= dps:
                self.x = [x for x in xs]
                self.w = [w for w in ws]
                return

        if verbose:
            print ("calculating nodes for degree-%i Gauss-Legendre "
                "quadrature..." % n)

        Float.store()
        Float.setdps(2*prec + 5)

        self.x = [None] * n
        self.w = [None] * n

        pf = polyfunc(legendre(n, 'x'), True)

        for k in xrange(n//2 + 1):
            if verbose and k % 4 == 0:
                print "  node", k, "of", n//2
            x, w = _gaussnode(pf, k, n)
            self.x[k] = x
            self.x[n-k-1] = -x
            self.w[k] = self.w[n-k-1] = w

        _gausscache[n] = (dps, self.x, self.w)

        Float.revert()