Added python scripts

scripts/gauss_chebyshev_t.py

@@ -0,0 +1,90 @@
+
+import numpy as np
+import sympy as sp
+from sympy.core import S, Dummy, pi
+from sympy.functions.elementary.trigonometric import sin, cos
+from sympy.functions.special.gamma_functions import gamma
+from sympy.polys.orthopolys import (legendre_poly, laguerre_poly,
+                                    hermite_poly, jacobi_poly)
+from sympy.polys.rootoftools import RootOf
+
+
+def gauss_chebyshev_t(n, n_digits):
+    r"""
+    Computes the Gauss-Chebyshev quadrature [1]_ points and weights of
+    the first kind.
+
+    The Gauss-Chebyshev quadrature of the first kind approximates the integral:
+
+    .. math::
+        \int_{-1}^{1} \frac{1}{\sqrt{1-x^2}} f(x)\,dx \approx
+            \sum_{i=1}^n w_i f(x_i)
+
+    The nodes `x_i` of an order `n` quadrature rule are the roots of `T_n`
+    and the weights `w_i` are given by:
+
+    .. math::
+        w_i = \frac{\pi}{n}
+
+    Parameters
+    ==========
+
+    n : the order of quadrature
+
+    n_digits : number of significant digits of the points and weights to return
+
+    Returns
+    =======
+
+    (x, w) : the ``x`` and ``w`` are lists of points and weights as Floats.
+             The points `x_i` and weights `w_i` are returned as ``(x, w)``
+             tuple of lists.
+
+    Examples
+    ========
+
+    >>> from sympy import S
+    >>> from sympy.integrals.quadrature import gauss_chebyshev_t
+    >>> x, w = gauss_chebyshev_t(3, 5)
+    >>> x
+    [0.86602, 0, -0.86602]
+    >>> w
+    [1.0472, 1.0472, 1.0472]
+
+    >>> x, w = gauss_chebyshev_t(6, 5)
+    >>> x
+    [0.96593, 0.70711, 0.25882, -0.25882, -0.70711, -0.96593]
+    >>> w
+    [0.5236, 0.5236, 0.5236, 0.5236, 0.5236, 0.5236]
+
+    See Also
+    ========
+
+    gauss_legendre, gauss_laguerre, gauss_hermite, gauss_gen_laguerre, gauss_chebyshev_u, gauss_jacobi, gauss_lobatto
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature
+    .. [2] http://people.sc.fsu.edu/~jburkardt/cpp_src/chebyshev1_rule/chebyshev1_rule.html
+    """
+    xi = []
+    w = []
+    for i in range(1, n+1):
+        xi.append((cos((2*i-S.One)/(2*n)*S.Pi)).n(n_digits))
+        w.append((S.Pi/n).n(n_digits))
+    return xi, w
+
+
+
+digits = 36
+npoints = list(range(1,17))
+
+for n in npoints:
+    xi, wi = gauss_chebyshev_t(n, digits)
+
+    print('\nN = %3d'%(n))
+    for i in range(0,n):
+        print('%3d %40.36f %40.36f'%(i+1,xi[i],wi[i]))
+
+

scripts/gauss_cheyshev_u.py

@@ -0,0 +1,88 @@
+
+import numpy as np
+import sympy as sp
+from sympy.core import S, Dummy, pi
+from sympy.functions.elementary.trigonometric import sin, cos
+from sympy.functions.special.gamma_functions import gamma
+from sympy.polys.orthopolys import (legendre_poly, laguerre_poly,
+                                    hermite_poly, jacobi_poly)
+from sympy.polys.rootoftools import RootOf
+
+
+def gauss_chebyshev_u(n, n_digits):
+    r"""
+    Computes the Gauss-Chebyshev quadrature [1]_ points and weights of
+    the second kind.
+
+    The Gauss-Chebyshev quadrature of the second kind approximates the
+    integral:
+
+    .. math::
+        \int_{-1}^{1} \sqrt{1-x^2} f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)
+
+    The nodes `x_i` of an order `n` quadrature rule are the roots of `U_n`
+    and the weights `w_i` are given by:
+
+    .. math::
+        w_i = \frac{\pi}{n+1} \sin^2 \left(\frac{i}{n+1}\pi\right)
+
+    Parameters
+    ==========
+
+    n : the order of quadrature
+
+    n_digits : number of significant digits of the points and weights to return
+
+    Returns
+    =======
+
+    (x, w) : the ``x`` and ``w`` are lists of points and weights as Floats.
+             The points `x_i` and weights `w_i` are returned as ``(x, w)``
+             tuple of lists.
+
+    Examples
+    ========
+
+    >>> from sympy import S
+    >>> from sympy.integrals.quadrature import gauss_chebyshev_u
+    >>> x, w = gauss_chebyshev_u(3, 5)
+    >>> x
+    [0.70711, 0, -0.70711]
+    >>> w
+    [0.3927, 0.7854, 0.3927]
+
+    >>> x, w = gauss_chebyshev_u(6, 5)
+    >>> x
+    [0.90097, 0.62349, 0.22252, -0.22252, -0.62349, -0.90097]
+    >>> w
+    [0.084489, 0.27433, 0.42658, 0.42658, 0.27433, 0.084489]
+
+    See Also
+    ========
+
+    gauss_legendre, gauss_laguerre, gauss_hermite, gauss_gen_laguerre, gauss_chebyshev_t, gauss_jacobi, gauss_lobatto
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature
+    .. [2] http://people.sc.fsu.edu/~jburkardt/cpp_src/chebyshev2_rule/chebyshev2_rule.html
+    """
+    xi = []
+    w = []
+    for i in range(1, n+1):
+        xi.append((cos(i/(n+S.One)*S.Pi)).n(n_digits))
+        w.append((S.Pi/(n+S.One)*sin(i*S.Pi/(n+S.One))**2).n(n_digits))
+    return xi, w
+
+
+
+digits = 36
+npoints = list(range(1,17))
+
+for n in npoints:
+    xi, wi = gauss_chebyshev_u(n, digits)
+
+    print('\nN = %3d'%(n))
+    for i in range(0,n):
+        print('%3d %40.36f %40.36f'%(i+1,xi[i],wi[i]))

scripts/gauss_gen_laguerre.py

@@ -0,0 +1,100 @@
+
+import numpy as np
+import sympy as sp
+from sympy.core import S, Dummy
+from sympy.functions.special.gamma_functions import gamma
+from sympy.polys.orthopolys import (legendre_poly, laguerre_poly,
+                                    hermite_poly, jacobi_poly)
+from sympy.polys.rootoftools import RootOf
+
+
+
+def gauss_gen_laguerre(n, alpha, n_digits):
+    r"""
+    Computes the generalized Gauss-Laguerre quadrature [1]_ points and weights.
+
+    The generalized Gauss-Laguerre quadrature approximates the integral:
+
+    .. math::
+        \int_{0}^\infty x^{\alpha} e^{-x} f(x)\,dx \approx
+            \sum_{i=1}^n w_i f(x_i)
+
+    The nodes `x_i` of an order `n` quadrature rule are the roots of
+    `L^{\alpha}_n` and the weights `w_i` are given by:
+
+    .. math::
+        w_i = \frac{\Gamma(\alpha+n)}
+                {n \Gamma(n) L^{\alpha}_{n-1}(x_i) L^{\alpha+1}_{n-1}(x_i)}
+
+    Parameters
+    ==========
+
+    n : the order of quadrature
+
+    alpha : the exponent of the singularity, `\alpha > -1`
+
+    n_digits : number of significant digits of the points and weights to return
+
+    Returns
+    =======
+
+    (x, w) : the ``x`` and ``w`` are lists of points and weights as Floats.
+             The points `x_i` and weights `w_i` are returned as ``(x, w)``
+             tuple of lists.
+
+    Examples
+    ========
+
+    >>> from sympy import S
+    >>> from sympy.integrals.quadrature import gauss_gen_laguerre
+    >>> x, w = gauss_gen_laguerre(3, -S.Half, 5)
+    >>> x
+    [0.19016, 1.7845, 5.5253]
+    >>> w
+    [1.4493, 0.31413, 0.00906]
+
+    >>> x, w = gauss_gen_laguerre(4, 3*S.Half, 5)
+    >>> x
+    [0.97851, 2.9904, 6.3193, 11.712]
+    >>> w
+    [0.53087, 0.67721, 0.11895, 0.0023152]
+
+    See Also
+    ========
+
+    gauss_legendre, gauss_laguerre, gauss_hermite, gauss_chebyshev_t, gauss_chebyshev_u, gauss_jacobi, gauss_lobatto
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature
+    .. [2] http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.html
+    """
+    x = Dummy("x")
+    p = laguerre_poly(n, x, alpha=alpha, polys=True)
+    p1 = laguerre_poly(n-1, x, alpha=alpha, polys=True)
+    p2 = laguerre_poly(n-1, x, alpha=alpha+1, polys=True)
+    xi = []
+    w = []
+    for r in p.real_roots():
+        if isinstance(r, RootOf):
+            r = r.eval_rational(S(1)/10**(n_digits+2))
+        xi.append(r.n(n_digits))
+        w.append((gamma(alpha+n) /
+                 (n*gamma(n)*p1.subs(x, r)*p2.subs(x, r))).n(n_digits))
+    return xi, w
+
+
+
+alpha = 0
+digits = 36
+npoints = list(range(1,17))
+
+for n in npoints:
+    xi, wi = gauss_gen_laguerre(n, alpha, digits)
+
+    print('\nN = %3d'%(n))
+    for i in range(0,n):
+        print('%3d %40.36f %40.36f'%(i+1,xi[i],wi[i]))
+
+

scripts/gauss_hermite.py

@@ -0,0 +1,97 @@
+
+import numpy as np
+import sympy as sp
+from sympy.core import S, Dummy, pi
+from sympy.functions.combinatorial.factorials import factorial
+from sympy.functions.special.gamma_functions import gamma
+from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.polys.orthopolys import (legendre_poly, laguerre_poly,
+                                    hermite_poly, jacobi_poly)
+from sympy.polys.rootoftools import RootOf
+
+
+
+def gauss_hermite(n, n_digits):
+    r"""
+    Computes the Gauss-Hermite quadrature [1]_ points and weights.
+
+    The Gauss-Hermite quadrature approximates the integral:
+
+    .. math::
+        \int_{-\infty}^{\infty} e^{-x^2} f(x)\,dx \approx
+            \sum_{i=1}^n w_i f(x_i)
+
+    The nodes `x_i` of an order `n` quadrature rule are the roots of `H_n`
+    and the weights `w_i` are given by:
+
+    .. math::
+        w_i = \frac{2^{n-1} n! \sqrt{\pi}}{n^2 \left(H_{n-1}(x_i)\right)^2}
+
+    Parameters
+    ==========
+
+    n : the order of quadrature
+
+    n_digits : number of significant digits of the points and weights to return
+
+    Returns
+    =======
+
+    (x, w) : the ``x`` and ``w`` are lists of points and weights as Floats.
+             The points `x_i` and weights `w_i` are returned as ``(x, w)``
+             tuple of lists.
+
+    Examples
+    ========
+
+    >>> from sympy.integrals.quadrature import gauss_hermite
+    >>> x, w = gauss_hermite(3, 5)
+    >>> x
+    [-1.2247, 0, 1.2247]
+    >>> w
+    [0.29541, 1.1816, 0.29541]
+
+    >>> x, w = gauss_hermite(6, 5)
+    >>> x
+    [-2.3506, -1.3358, -0.43608, 0.43608, 1.3358, 2.3506]
+    >>> w
+    [0.00453, 0.15707, 0.72463, 0.72463, 0.15707, 0.00453]
+
+    See Also
+    ========
+
+    gauss_legendre, gauss_laguerre, gauss_gen_laguerre, gauss_chebyshev_t, gauss_chebyshev_u, gauss_jacobi, gauss_lobatto
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Gauss-Hermite_Quadrature
+    .. [2] http://people.sc.fsu.edu/~jburkardt/cpp_src/hermite_rule/hermite_rule.html
+    .. [3] http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_hermite_rule/gen_hermite_rule.html
+    """
+    x = Dummy("x")
+    p = hermite_poly(n, x, polys=True)
+    p1 = hermite_poly(n-1, x, polys=True)
+    xi = []
+    w = []
+    for r in p.real_roots():
+        if isinstance(r, RootOf):
+            r = r.eval_rational(S(1)/10**(n_digits+2))
+        xi.append(r.n(n_digits))
+        w.append(((2**(n-1) * factorial(n) * sqrt(pi)) /
+                 (n**2 * p1.subs(x, r)**2)).n(n_digits))
+    return xi, w
+
+
+
+digits = 36
+npoints = list(range(1,17))
+
+for n in npoints:
+    xi, wi = gauss_hermite(n, digits)
+
+    print('\nN = %3d'%(n))
+    for i in range(0,n):
+        print('%3d %40.36f %40.36f'%(i+1,xi[i],wi[i]))
+
+