diff options
| -rw-r--r-- | numpy/polynomial/chebyshev.py | 72 | ||||
| -rw-r--r-- | numpy/polynomial/hermite.py | 93 | ||||
| -rw-r--r-- | numpy/polynomial/hermite_e.py | 92 | ||||
| -rw-r--r-- | numpy/polynomial/laguerre.py | 89 | ||||
| -rw-r--r-- | numpy/polynomial/legendre.py | 92 |
5 files changed, 433 insertions, 5 deletions
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py index c80123bd3..dbadaa67b 100644 --- a/numpy/polynomial/chebyshev.py +++ b/numpy/polynomial/chebyshev.py @@ -42,6 +42,9 @@ Misc Functions - `chebvander` -- Vandermonde-like matrix for Chebyshev polynomials. - `chebvander2d` -- Vandermonde-like matrix for 2D power series. - `chebvander3d` -- Vandermonde-like matrix for 3D power series. +- `chebgauss` -- Gauss-Chebyshev quadrature, points and weights. +- `chebweight` -- Chebyshev weight function. +- `chebcompanion` -- symmetrized companion matrix in Chebyshev form. - `chebfit` -- least-squares fit returning a Chebyshev series. - `chebpts1` -- Chebyshev points of the first kind. - `chebpts2` -- Chebyshev points of the second kind. @@ -95,7 +98,8 @@ __all__ = ['chebzero', 'chebone', 'chebx', 'chebdomain', 'chebline', 'chebval', 'chebder', 'chebint', 'cheb2poly', 'poly2cheb', 'chebfromroots', 'chebvander', 'chebfit', 'chebtrim', 'chebroots', 'chebpts1', 'chebpts2', 'Chebyshev', 'chebval2d', 'chebval3d', - 'chebgrid2d', 'chebgrid3d', 'chebvander2d','chebvander3d'] + 'chebgrid2d', 'chebgrid3d', 'chebvander2d','chebvander3d', + 'chebcompanion', 'chebgauss', 'chebweight'] chebtrim = pu.trimcoef @@ -1724,6 +1728,72 @@ def chebroots(cs): return r +def chebgauss(deg): + """Gauss Chebyshev quadrature. + + Computes the sample points and weights for Gauss-Chebyshev quadrature. + These sample points and weights will correctly integrate polynomials of + degree ``2*deg - 1`` or less over the interval ``[-1, 1]`` with the + weight function ``f(x) = 1/sqrt(1 - x**2)``. + + Parameters + ---------- + deg : int + Number of sample points and weights. It must be >= 1. + + Returns + ------- + x : ndarray + 1-D ndarray containing the sample points. + y : ndarray + 1-D ndarray containing the weights. + + Notes + ----- + The results have only been tested up to degree 100. Higher degrees may + be problematic. There are closed form solutions for the sample points + and weights. If ``n = deg``, then + + ``x_i = cos(pi*(2*i - 1)/(2*n))`` + ``w_i = pi/n`` + + where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th + root of ``L_n``, and then scaling the results to get the right value + when integrating 1. + + """ + ideg = int(deg) + if ideg != deg or ideg < 1: + raise ValueError("deg must be a non-negative integer") + + x = np.cos(np.pi * np.arange(1, 2*ideg, 2) / (2.0*ideg)) + w = np.ones(ideg)*(np.pi/ideg) + + return x, w + + +def chebweight(x): + """Weight function of the Chebyshev polynomials. + + The weight function for which the Chebyshev polynomials are orthogonal. + In this case the weight function is ``1/(1 - x**2)``. Note that the + Chebyshev polynomials are not normalized. + + Parameters + ---------- + x : array_like + Values at which the weight function will be computed. + + Returns + ------- + w : ndarray + The weight function at `x`. + + """ + w = 1./(np.sqrt(1. + x) * np.sqrt(1. - x)) + return w + + def chebpts1(npts): """Chebyshev points of the first kind. diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py index 5babfb7b1..58cc655d7 100644 --- a/numpy/polynomial/hermite.py +++ b/numpy/polynomial/hermite.py @@ -39,6 +39,9 @@ Misc Functions - `hermvander` -- Vandermonde-like matrix for Hermite polynomials. - `hermvander2d` -- Vandermonde-like matrix for 2D power series. - `hermvander3d` -- Vandermonde-like matrix for 3D power series. +- `hermgauss` -- Gauss-Hermite quadrature, points and weights. +- `hermweight` -- Hermite weight function. +- `hermcompanion` -- symmetrized companion matrix in Hermite form. - `hermfit` -- least-squares fit returning a Hermite series. - `hermtrim` -- trim leading coefficients from a Hermite series. - `hermline` -- Hermite series of given straight line. @@ -67,7 +70,8 @@ __all__ = ['hermzero', 'hermone', 'hermx', 'hermdomain', 'hermline', 'hermval', 'hermder', 'hermint', 'herm2poly', 'poly2herm', 'hermfromroots', 'hermvander', 'hermfit', 'hermtrim', 'hermroots', 'Hermite', 'hermval2d', 'hermval3d', 'hermgrid2d', 'hermgrid3d', - 'hermvander2d', 'hermvander3d'] + 'hermvander2d', 'hermvander3d', 'hermcompanion', 'hermgauss', + 'hermweight'] hermtrim = pu.trimcoef @@ -1507,6 +1511,93 @@ def hermroots(cs): return r +def hermgauss(deg): + """Gauss Hermite quadrature. + + Computes the sample points and weights for Gauss-Hermite quadrature. + These sample points and weights will correctly integrate polynomials of + degree ``2*deg - 1`` or less over the interval ``[-inf, inf]`` with the + weight function ``f(x) = exp(-x**2)``. + + Parameters + ---------- + deg : int + Number of sample points and weights. It must be >= 1. + + Returns + ------- + x : ndarray + 1-D ndarray containing the sample points. + y : ndarray + 1-D ndarray containing the weights. + + Notes + ----- + The results have only been tested up to degree 100. Higher degrees may + be problematic. The weights are determined by using the fact that + + w = c / (H'_n(x_k) * H_{n-1}(x_k)) + + where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th + root of ``H_n``, and then scaling the results to get the right value + when integrating 1. + + """ + ideg = int(deg) + if ideg != deg or ideg < 1: + raise ValueError("deg must be a non-negative integer") + + # first approximation of roots. We use the fact that the companion + # matrix is symmetric in this case in order to obtain better zeros. + c = np.array([0]*deg + [1]) + m = hermcompanion(c) + x = la.eigvals(m) + x.sort() + + # improve roots by one application of Newton + dy = hermval(x, c) + df = hermval(x, hermder(c)) + x -= dy/df + + # compute the weights. We scale the factor to avoid possible numerical + # overflow. + fm = hermval(x, c[1:]) + fm /= np.abs(fm).max() + df /= np.abs(df).max() + w = 1/(fm * df) + + # for Hermite we can also symmetrize + w = (w + w[::-1])/2 + x = (x - x[::-1])/2 + + # scale w to get the right value + w *= np.sqrt(np.pi) / w.sum() + + return x, w + + +def hermweight(x): + """Weight function of the Hermite polynomials. + + The weight function for which the Hermite polynomials are orthogonal. + In this case the weight function is ``exp(-x**2)``. Note that the + Hermite polynomials are not normalized. + + Parameters + ---------- + x : array_like + Values at which the weight function will be computed. + + Returns + ------- + w : ndarray + The weight function at `x`. + + """ + w = np.exp(-x**2) + return w + + # # Hermite series class # diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py index 97ef366bc..c8d641734 100644 --- a/numpy/polynomial/hermite_e.py +++ b/numpy/polynomial/hermite_e.py @@ -39,6 +39,9 @@ Misc Functions - `hermevander` -- Vandermonde-like matrix for Hermite_e polynomials. - `hermevander2d` -- Vandermonde-like matrix for 2D power series. - `hermevander3d` -- Vandermonde-like matrix for 3D power series. +- `hermegauss` -- Gauss-Hermite_e quadrature, points and weights. +- `hermeweight` -- Hermite_e weight function. +- `hermecompanion` -- symmetrized companion matrix in Hermite_e form. - `hermefit` -- least-squares fit returning a Hermite_e series. - `hermetrim` -- trim leading coefficients from a Hermite_e series. - `hermeline` -- Hermite_e series of given straight line. @@ -68,7 +71,7 @@ __all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline', 'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots', 'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'HermiteE', 'hermeval2d', 'hermeval3d', 'hermegrid2d', 'hermegrid3d', 'hermevander2d', - 'hermevander3d'] + 'hermevander3d', 'hermecompanion', 'hermegauss', 'hermeweight'] hermetrim = pu.trimcoef @@ -1503,6 +1506,93 @@ def hermeroots(cs): return r +def hermegauss(deg): + """Gauss Hermite_e quadrature. + + Computes the sample points and weights for Gauss-Hermite_e quadrature. + These sample points and weights will correctly integrate polynomials of + degree ``2*deg - 1`` or less over the interval ``[-inf, inf]`` with the + weight function ``f(x) = exp(-.5*x**2)``. + + Parameters + ---------- + deg : int + Number of sample points and weights. It must be >= 1. + + Returns + ------- + x : ndarray + 1-D ndarray containing the sample points. + y : ndarray + 1-D ndarray containing the weights. + + Notes + ----- + The results have only been tested up to degree 100. Higher degrees may + be problematic. The weights are determined by using the fact that + + w = c / (He'_n(x_k) * He_{n-1}(x_k)) + + where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th + root of ``He_n``, and then scaling the results to get the right value + when integrating 1. + + """ + ideg = int(deg) + if ideg != deg or ideg < 1: + raise ValueError("deg must be a non-negative integer") + + # first approximation of roots. We use the fact that the companion + # matrix is symmetric in this case in order to obtain better zeros. + c = np.array([0]*deg + [1]) + m = hermecompanion(c) + x = la.eigvals(m) + x.sort() + + # improve roots by one application of Newton + dy = hermeval(x, c) + df = hermeval(x, hermeder(c)) + x -= dy/df + + # compute the weights. We scale the factor to avoid possible numerical + # overflow. + fm = hermeval(x, c[1:]) + fm /= np.abs(fm).max() + df /= np.abs(df).max() + w = 1/(fm * df) + + # for Hermite_e we can also symmetrize + w = (w + w[::-1])/2 + x = (x - x[::-1])/2 + + # scale w to get the right value + w *= np.sqrt(2*np.pi) / w.sum() + + return x, w + + +def hermeweight(x): + """Weight function of the Hermite_e polynomials. + + The weight function for which the Hermite_e polynomials are orthogonal. + In this case the weight function is ``exp(-.5*x**2)``. Note that the + Hermite_e polynomials are not normalized. + + Parameters + ---------- + x : array_like + Values at which the weight function will be computed. + + Returns + ------- + w : ndarray + The weight function at `x`. + + """ + w = np.exp(-.5*x**2) + return w + + # # HermiteE series class # diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py index 5ab2d23ae..872829fbc 100644 --- a/numpy/polynomial/laguerre.py +++ b/numpy/polynomial/laguerre.py @@ -39,6 +39,9 @@ Misc Functions - `lagvander` -- Vandermonde-like matrix for Laguerre polynomials. - `lagvander2d` -- Vandermonde-like matrix for 2D power series. - `lagvander3d` -- Vandermonde-like matrix for 3D power series. +- `laggauss` -- Gauss-Laguerre quadrature, points and weights. +- `lagweight` -- Laguerre weight function. +- `lagcompanion` -- symmetrized companion matrix in Laguerre form. - `lagfit` -- least-squares fit returning a Laguerre series. - `lagtrim` -- trim leading coefficients from a Laguerre series. - `lagline` -- Laguerre series of given straight line. @@ -66,7 +69,8 @@ __all__ = ['lagzero', 'lagone', 'lagx', 'lagdomain', 'lagline', 'lagadd', 'lagsub', 'lagmulx', 'lagmul', 'lagdiv', 'lagpow', 'lagval', 'lagder', 'lagint', 'lag2poly', 'poly2lag', 'lagfromroots', 'lagvander', 'lagfit', 'lagtrim', 'lagroots', 'Laguerre', 'lagval2d', - 'lagval3d', 'laggrid2d', 'laggrid3d', 'lagvander2d', 'lagvander3d'] + 'lagval3d', 'laggrid2d', 'laggrid3d', 'lagvander2d', 'lagvander3d', + 'lagcompanion', 'laggauss', 'lagweight'] lagtrim = pu.trimcoef @@ -1504,6 +1508,89 @@ def lagroots(cs): return r +def laggauss(deg): + """Gauss Laguerre quadrature. + + Computes the sample points and weights for Gauss-Laguerre quadrature. + These sample points and weights will correctly integrate polynomials of + degree ``2*deg - 1`` or less over the interval ``[0, inf]`` with the + weight function ``f(x) = exp(-x)``. + + Parameters + ---------- + deg : int + Number of sample points and weights. It must be >= 1. + + Returns + ------- + x : ndarray + 1-D ndarray containing the sample points. + y : ndarray + 1-D ndarray containing the weights. + + Notes + ----- + The results have only been tested up to degree 100. Higher degrees may + be problematic. The weights are determined by using the fact that + + w = c / (L'_n(x_k) * L_{n-1}(x_k)) + + where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th + root of ``L_n``, and then scaling the results to get the right value + when integrating 1. + + """ + ideg = int(deg) + if ideg != deg or ideg < 1: + raise ValueError("deg must be a non-negative integer") + + # first approximation of roots. We use the fact that the companion + # matrix is symmetric in this case in order to obtain better zeros. + c = np.array([0]*deg + [1]) + m = lagcompanion(c) + x = la.eigvals(m) + x.sort() + + # improve roots by one application of Newton + dy = lagval(x, c) + df = lagval(x, lagder(c)) + x -= dy/df + + # compute the weights. We scale the factor to avoid possible numerical + # overflow. + fm = lagval(x, c[1:]) + fm /= np.abs(fm).max() + df /= np.abs(df).max() + w = 1/(fm * df) + + # scale w to get the right value, 1 in this case + w /= w.sum() + + return x, w + + +def lagweight(x): + """Weight function of the Laguerre polynomials. + + The weight function for which the Laguerre polynomials are orthogonal. + In this case the weight function is ``exp(-x)``. Note that the Laguerre + polynomials are not normalized, indeed, may be much greater than the + normalized versions. + + Parameters + ---------- + x : array_like + Values at which the weight function will be computed. + + Returns + ------- + w : ndarray + The weight function at `x`. + + """ + w = np.exp(-x) + return w + # # Laguerre series class # diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py index b9a74501c..93ba1d348 100644 --- a/numpy/polynomial/legendre.py +++ b/numpy/polynomial/legendre.py @@ -40,6 +40,9 @@ Misc Functions - `legvander` -- Vandermonde-like matrix for Legendre polynomials. - `legvander2d` -- Vandermonde-like matrix for 2D power series. - `legvander3d` -- Vandermonde-like matrix for 3D power series. +- `leggauss` -- Gauss-Legendre quadrature, points and weights. +- `legweight` -- Legendre weight function. +- `legcompanion` -- symmetrized companion matrix in Legendre form. - `legfit` -- least-squares fit returning a Legendre series. - `legtrim` -- trim leading coefficients from a Legendre series. - `legline` -- Legendre series representing given straight line. @@ -67,7 +70,8 @@ __all__ = ['legzero', 'legone', 'legx', 'legdomain', 'legline', 'legadd', 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow', 'legval', 'legder', 'legint', 'leg2poly', 'poly2leg', 'legfromroots', 'legvander', 'legfit', 'legtrim', 'legroots', 'Legendre','legval2d', - 'legval3d', 'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d'] + 'legval3d', 'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d', + 'legcompanion', 'leggauss', 'legweight'] legtrim = pu.trimcoef @@ -1510,6 +1514,92 @@ def legroots(cs): return r +def leggauss(deg): + """Gauss Legendre quadrature. + + Computes the sample points and weights for Gauss-Legendre quadrature. + These sample points and weights will correctly integrate polynomials of + degree ``2*deg - 1`` or less over the interval ``[-1, 1]`` with the + weight function ``f(x) = 1``. + + Parameters + ---------- + deg : int + Number of sample points and weights. It must be >= 1. + + Returns + ------- + x : ndarray + 1-D ndarray containing the sample points. + y : ndarray + 1-D ndarray containing the weights. + + Notes + ----- + The results have only been tested up to degree 100. Higher degrees may + be problematic. The weights are determined by using the fact that + + w = c / (L'_n(x_k) * L_{n-1}(x_k)) + + where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th + root of ``L_n``, and then scaling the results to get the right value + when integrating 1. + + """ + ideg = int(deg) + if ideg != deg or ideg < 1: + raise ValueError("deg must be a non-negative integer") + + # first approximation of roots. We use the fact that the companion + # matrix is symmetric in this case in order to obtain better zeros. + c = np.array([0]*deg + [1]) + m = legcompanion(c) + x = la.eigvals(m) + x.sort() + + # improve roots by one application of Newton + dy = legval(x, c) + df = legval(x, legder(c)) + x -= dy/df + + # compute the weights. We scale the factor to avoid possible numerical + # overflow. + fm = legval(x, c[1:]) + fm /= np.abs(fm).max() + df /= np.abs(df).max() + w = 1/(fm * df) + + # for Legendre we can also symmetrize + w = (w + w[::-1])/2 + x = (x - x[::-1])/2 + + # scale w to get the right value + w *= 2. / w.sum() + + return x, w + + +def legweight(x): + """Weight function of the Legendre polynomials. + + The weight function for which the Legendre polynomials are orthogonal. + In this case the weight function is simply one. Note that the Legendre + polynomials are not normalized. + + Parameters + ---------- + x : array_like + Values at which the weight function will be computed. + + Returns + ------- + w : ndarray + The weight function at `x`. + + """ + w = x*0.0 + 1.0 + return w + # # Legendre series class # |