diff options
Diffstat (limited to 'numpy/polynomial')
| -rw-r--r-- | numpy/polynomial/chebyshev.py | 427 | ||||
| -rw-r--r-- | numpy/polynomial/hermite.py | 443 | ||||
| -rw-r--r-- | numpy/polynomial/hermite_e.py | 444 | ||||
| -rw-r--r-- | numpy/polynomial/laguerre.py | 443 | ||||
| -rw-r--r-- | numpy/polynomial/legendre.py | 441 | ||||
| -rw-r--r-- | numpy/polynomial/polynomial.py | 359 |
6 files changed, 1408 insertions, 1149 deletions
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py index dbadaa67b..eb0087395 100644 --- a/numpy/polynomial/chebyshev.py +++ b/numpy/polynomial/chebyshev.py @@ -108,7 +108,7 @@ chebtrim = pu.trimcoef # functions and do minimal error checking. # -def _cseries_to_zseries(cs) : +def _cseries_to_zseries(c) : """Covert Chebyshev series to z-series. Covert a Chebyshev series to the equivalent z-series. The result is @@ -118,7 +118,7 @@ def _cseries_to_zseries(cs) : Parameters ---------- - cs : 1-d ndarray + c : 1-d ndarray Chebyshev coefficients, ordered from low to high Returns @@ -127,9 +127,9 @@ def _cseries_to_zseries(cs) : Odd length symmetric z-series, ordered from low to high. """ - n = cs.size - zs = np.zeros(2*n-1, dtype=cs.dtype) - zs[n-1:] = cs/2 + n = c.size + zs = np.zeros(2*n-1, dtype=c.dtype) + zs[n-1:] = c/2 return zs + zs[::-1] @@ -148,14 +148,14 @@ def _zseries_to_cseries(zs) : Returns ------- - cs : 1-d ndarray + c : 1-d ndarray Chebyshev coefficients, ordered from low to high. """ n = (zs.size + 1)//2 - cs = zs[n-1:].copy() - cs[1:n] *= 2 - return cs + c = zs[n-1:].copy() + c[1:n] *= 2 + return c def _zseries_mul(z1, z2) : @@ -338,7 +338,7 @@ def poly2cheb(pol) : Returns ------- - cs : ndarray + c : ndarray 1-d array containing the coefficients of the equivalent Chebyshev series. @@ -372,7 +372,7 @@ def poly2cheb(pol) : return res -def cheb2poly(cs) : +def cheb2poly(c) : """ Convert a Chebyshev series to a polynomial. @@ -383,7 +383,7 @@ def cheb2poly(cs) : Parameters ---------- - cs : array_like + c : array_like 1-d array containing the Chebyshev series coefficients, ordered from lowest order term to highest. @@ -418,17 +418,17 @@ def cheb2poly(cs) : """ from polynomial import polyadd, polysub, polymulx - [cs] = pu.as_series([cs]) - n = len(cs) + [c] = pu.as_series([c]) + n = len(c) if n < 3: - return cs + return c else: - c0 = cs[-2] - c1 = cs[-1] + c0 = c[-2] + c1 = c[-1] # i is the current degree of c1 for i in range(n - 1, 1, -1) : tmp = c0 - c0 = polysub(cs[i - 2], c1) + c0 = polysub(c[i - 2], c1) c1 = polyadd(tmp, polymulx(c1)*2) return polyadd(c0, polymulx(c1)) @@ -489,14 +489,24 @@ def chebline(off, scl) : def chebfromroots(roots) : """ - Generate a Chebyshev series with the given roots. + Generate a Chebyshev series with given roots. + + The function returns the coefficients of the polynomial - Return the array of coefficients for the C-series whose roots (a.k.a. - "zeros") are given by *roots*. The returned array of coefficients is - ordered from lowest order "term" to highest, and zeros of multiplicity - greater than one must be included in *roots* a number of times equal - to their multiplicity (e.g., if `2` is a root of multiplicity three, - then [2,2,2] must be in *roots*). + .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), + + in Chebyshev form, where the `r_n` are the roots specified in `roots`. + If a zero has multiplicity n, then it must appear in `roots` n times. + For instance, if 2 is a root of multiplicity three and 3 is a root of + multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The + roots can appear in any order. + + If the returned coefficients are `c`, then + + .. math:: p(x) = c_0 + c_1 * T_1(x) + ... + c_n * T_n(x) + + The coefficient of the last term is not generally 1 for monic + polynomials in Chebyshev form. Parameters ---------- @@ -506,28 +516,15 @@ def chebfromroots(roots) : Returns ------- out : ndarray - 1-d array of the C-series' coefficients, ordered from low to - high. If all roots are real, ``out.dtype`` is a float type; - otherwise, ``out.dtype`` is a complex type, even if all the - coefficients in the result are real (see Examples below). + 1-D array of coefficients. If all roots are real then `out` is a + real array, if some of the roots are complex, then `out` is complex + even if all the coefficients in the result are real (see Examples + below). See Also -------- - polyfromroots - - Notes - ----- - What is returned are the :math:`c_i` such that: - - .. math:: - - \\sum_{i=0}^{n} c_i*T_i(x) = \\prod_{i=0}^{n} (x - roots[i]) - - where ``n == len(roots)`` and :math:`T_i(x)` is the `i`-th Chebyshev - (basis) polynomial over the domain `[-1,1]`. Note that, unlike - `polyfromroots`, due to the nature of the C-series basis set, the - above identity *does not* imply :math:`c_n = 1` identically (see - Examples). + polyfromroots, legfromroots, lagfromroots, hermfromroots, + hermefromroots. Examples -------- @@ -652,16 +649,16 @@ def chebsub(c1, c2): return pu.trimseq(ret) -def chebmulx(cs): +def chebmulx(c): """Multiply a Chebyshev series by x. - Multiply the polynomial `cs` by x, where x is the independent + Multiply the polynomial `c` by x, where x is the independent variable. Parameters ---------- - cs : array_like + c : array_like 1-d array of Chebyshev series coefficients ordered from low to high. @@ -675,17 +672,17 @@ def chebmulx(cs): .. versionadded:: 1.5.0 """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) # The zero series needs special treatment - if len(cs) == 1 and cs[0] == 0: - return cs - - prd = np.empty(len(cs) + 1, dtype=cs.dtype) - prd[0] = cs[0]*0 - prd[1] = cs[0] - if len(cs) > 1: - tmp = cs[1:]/2 + if len(c) == 1 and c[0] == 0: + return c + + prd = np.empty(len(c) + 1, dtype=c.dtype) + prd[0] = c[0]*0 + prd[1] = c[0] + if len(c) > 1: + tmp = c[1:]/2 prd[2:] = tmp prd[0:-2] += tmp return prd @@ -806,16 +803,16 @@ def chebdiv(c1, c2): return quo, rem -def chebpow(cs, pow, maxpower=16) : +def chebpow(c, pow, maxpower=16) : """Raise a Chebyshev series to a power. - Returns the Chebyshev series `cs` raised to the power `pow`. The - arguement `cs` is a sequence of coefficients ordered from low to high. + Returns the Chebyshev series `c` raised to the power `pow`. The + arguement `c` is a sequence of coefficients ordered from low to high. i.e., [1,2,3] is the series ``T_0 + 2*T_1 + 3*T_2.`` Parameters ---------- - cs : array_like + c : array_like 1d array of chebyshev series coefficients ordered from low to high. pow : integer @@ -837,21 +834,21 @@ def chebpow(cs, pow, maxpower=16) : -------- """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) power = int(pow) if power != pow or power < 0 : raise ValueError("Power must be a non-negative integer.") elif maxpower is not None and power > maxpower : raise ValueError("Power is too large") elif power == 0 : - return np.array([1], dtype=cs.dtype) + return np.array([1], dtype=c.dtype) elif power == 1 : - return cs + return c else : # This can be made more efficient by using powers of two # in the usual way. - zs = _cseries_to_zseries(cs) + zs = _cseries_to_zseries(c) prd = zs for i in range(2, power + 1) : prd = np.convolve(prd, zs) @@ -1013,8 +1010,8 @@ def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0): Note that the result of each integration is *multiplied* by `scl`. Why is this important to note? Say one is making a linear change of variable :math:`u = ax + b` in an integral relative to `x`. Then - :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a` - - perhaps not what one would have first thought. + .. math::`dx = du/a`, so one will need to set `scl` equal to + :math:`1/a`- perhaps not what one would have first thought. Also note that, in general, the result of integrating a C-series needs to be "re-projected" onto the C-series basis set. Thus, typically, @@ -1086,48 +1083,48 @@ def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0): def chebval(x, c, tensor=True): """ - Evaluate a Chebyshev series. + Evaluate a Chebyshev series at points x. - If `c` is of length ``n + 1``, this function returns the value: + If `c` is of length `n + 1`, this function returns the value: - ``p(x) = c[0]*T_0(x) + c[1]*T_1(x) + ... + c[n]*T_n(x)`` + .. math:: p(x) = c_0 * T_0(x) + c_1 * T_1(x) + ... + c_n * T_n(x) - If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)`` - will have the same shape as `x`. If `x` is a algebra_like object that - supports multiplication and addition with itself and the values in `c`, - then an object of the same type is returned. + The parameter `x` is converted to an array only if it is a tuple or a + list, otherwise it is treated as a scalar. In either case, either `x` + or its elements must support multiplication and addition both with + themselves and with the elements of `c`. - In the case where c is multidimensional, the shape of the result - depends on the value of `tensor`. If tensor is true the shape of the - return will be ``c.shape[1:] + x.shape``, where the shape of a scalar - is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if - `x` is broadcast compatible with that. + If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If + `c` is multidimensional, then the shape of the result depends on the + value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + + x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that + scalars have shape (,). - If there are trailing zeros in the coefficients they still take part in - the evaluation, so they should be avoided if efficiency is a concern. + Trailing zeros in the coefficients will be used in the evaluation, so + they should be avoided if efficiency is a concern. Parameters ---------- - x : array_like, algebra_like - If x is a list or tuple, it is converted to an ndarray. Otherwise - it is left unchanged and if it isn't an ndarray it is treated as a - scalar. In either case, `x` or any element of an ndarray must - support addition and multiplication with itself and the elements of - `c`. + x : array_like, compatible object + If `x` is a list or tuple, it is converted to an ndarray, otherwise + it is left unchanged and treated as a scalar. In either case, `x` + or its elements must support addition and multiplication with + with themselves and with the elements of `c`. c : array_like Array of coefficients ordered so that the coefficients for terms of - degree n are contained in ``c[n]``. If `c` is multidimesional the - remaining indices enumerate multiple Chebyshev series. In the two + degree n are contained in c[n]. If `c` is multidimesional the + remaining indices enumerate multiple polynomials. In the two dimensional case the coefficients may be thought of as stored in the columns of `c`. tensor : boolean, optional - If true, the coefficient array shape is extended with ones on the - right, one for each dimension of `x`. Scalars are treated as having - dimension 0 for this action. The effect is that every column of - coefficients in `c` is evaluated for every value in `x`. If False, - the `x` are broadcast over the columns of `c` in the usual way. - This gives some flexibility in evaluations in the multidimensional - case. The default value it ``True``. + If True, the shape of the coefficient array is extended with ones + on the right, one for each dimension of `x`. Scalars have dimension 0 + for this action. The result is that every column of coefficients in + `c` is evaluated for every element of `x`. If False, `x` is broadcast + over the columns of `c` for the evaluation. This keyword is useful + when `c` is multidimensional. The default value is True. + + .. versionadded:: 1.7.0 Returns ------- @@ -1173,29 +1170,39 @@ def chebval(x, c, tensor=True): def chebval2d(x, y, c): """ - Evaluate 2D Chebyshev series at points (x,y). + Evaluate a 2-D Chebyshev series at points (x, y). This function returns the values: - ``p(x,y) = \sum_{i,j} c[i,j] * T_i(x) * T_j(y)`` + .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * T_i(x) * T_j(y) + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars and they + must have the same shape after conversion. In either case, either `x` + and `y` or their elements must support multiplication and addition both + with themselves and with the elements of `c`. + + If `c` is a 1-D array a one is implicitly appended to its shape to make + it 2-D. The shape of the result will be c.shape[2:] + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Chebyshev seres is evaluated at the points - ``(x,y)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray. - Otherwise it is left unchanged and if it isn't an ndarray it is - treated as a scalar. See `chebval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points `(x, y)`, + where `x` and `y` must have the same shape. If `x` or `y` is a list + or tuple, it is first converted to an ndarray, otherwise it is left + unchanged and if it isn't an ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of - coefficients. + Array of coefficients ordered so that the coefficient of the term + of multi-degree i,j is contained in ``c[i,j]``. If `c` has + dimension greater than 2 the remaining indices enumerate multiple + sets of coefficients. Returns ------- - values : ndarray, algebra_like + values : ndarray, compatible object The values of the two dimensional Chebyshev series at points formed from pairs of corresponding values from `x` and `y`. @@ -1216,33 +1223,43 @@ def chebval2d(x, y, c): def chebgrid2d(x, y, c): """ - Evaluate 2D Chebyshev series on the Cartesion product of x,y. + Evaluate a 2-D Chebyshev series on the Cartesion product of x and y. This function returns the values: - ``p(a,b) = \sum_{i,j} c[i,j] * T_i(a) * T_j(b)`` + .. math:: p(a,b) = \sum_{i,j} c_{i,j} * T_i(a) * T_j(b) + + where the points `(a, b)` consist of all pairs formed by taking + `a` from `x` and `b` from `y`. The resulting points form a grid with + `x` in the first dimension and `y` in the second. - where the points ``(a,b)`` consist of all pairs of points formed by - taking ``a`` from `x` and ``b`` from `y`. The resulting points form a - grid with `x` in the first dimension and `y` in the second. + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars. In either + case, either `x` and `y` or their elements must support multiplication + and addition both with themselves and with the elements of `c`. + + If `c` has fewer than two dimensions, ones are implicitly appended to + its shape to make it 2-D. The shape of the result will be c.shape[2:] + + x.shape + y.shape. + + .. versionadded:: 1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Chebyshev series is evaluated at the points in - the Cartesian product of `x` and `y`. If `x` or `y` is a list or - tuple, it is first converted to an ndarray, Otherwise it is left - unchanged and if it isn't an ndarray it is treated as a scalar. See - `chebval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points in the + Cartesian product of `x` and `y`. If `x` or `y` is a list or + tuple, it is first converted to an ndarray, otherwise it is left + unchanged and, if it isn't an ndarray, it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + Array of coefficients ordered so that the coefficient of the term of + multi-degree i,j is contained in `c[i,j]`. If `c` has dimension + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like + values : ndarray, compatible object The values of the two dimensional Chebyshev series at points in the Cartesion product of `x` and `y`. @@ -1258,32 +1275,43 @@ def chebgrid2d(x, y, c): def chebval3d(x, y, z, c): """ - Evaluate 3D Chebyshev series at points (x,y,z). + Evaluate a 3-D Chebyshev series at points (x, y, z). This function returns the values: - ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * T_i(x) * T_j(y) * T_k(z)`` + .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * T_i(x) * T_j(y) * T_k(z) + + The parameters `x`, `y`, and `z` are converted to arrays only if + they are tuples or a lists, otherwise they are treated as a scalars and + they must have the same shape after conversion. In either case, either + `x`, `y`, and `z` or their elements must support multiplication and + addition both with themselves and with the elements of `c`. + + If `c` has fewer than 3 dimensions, ones are implicitly appended to its + shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Chebyshev seres is evaluated at the points - ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape. If + x, y, z : array_like, compatible object + The three dimensional series is evaluated at the points + `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If any of `x`, `y`, or `z` is a list or tuple, it is first converted - to an ndarray. Otherwise it is left unchanged and if it isn't an - ndarray it is treated as a scalar. See `chebval` for explanation of - algebra_like. + to an ndarray, otherwise it is left unchanged and if it isn't an + ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + Array of coefficients ordered so that the coefficient of the term of + multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension + greater than 3 the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Chebyshev series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the multidimension polynomial on points formed with + triples of corresponding values from `x`, `y`, and `z`. See Also -------- @@ -1303,37 +1331,48 @@ def chebval3d(x, y, z, c): def chebgrid3d(x, y, z, c): """ - Evaluate 3D Chebyshev series on the Cartesian product of x,y,z. + Evaluate a 3-D Chebyshev series on the Cartesian product of x, y, and z. This function returns the values: - ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * T_i(a) * T_j(b) * T_k(c)`` + .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * T_i(a) * T_j(b) * T_k(c) + + where the points `(a, b, c)` consist of all triples formed by taking + `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form + a grid with `x` in the first dimension, `y` in the second, and `z` in + the third. - where the points ``(a,b,c)`` consist of all triples formed by taking - ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting - points form a grid with `x` in the first dimension, `y` in the second, - and `z` in the third. + The parameters `x`, `y`, and `z` are converted to arrays only if they + are tuples or a lists, otherwise they are treated as a scalars. In + either case, either `x`, `y`, and `z` or their elements must support + multiplication and addition both with themselves and with the elements + of `c`. + + If `c` has fewer than three dimensions, ones are implicitly appended to + its shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape + yshape + z.shape. + + .. versionadded:: 1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Chebyshev seres is evaluated at the points - in the cartesian product of `x`, `y`, and `z` - ``(x,y,z)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray, - otherwise it is left unchanged and treated as a scalar. See - `chebval` for explanation of algebra_like. + x, y, z : array_like, compatible objects + The three dimensional series is evaluated at the points in the + Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a + list or tuple, it is first converted to an ndarray, otherwise it is + left unchanged and, if it isn't an ndarray, it is treated as a + scalar. c : array_like Array of coefficients ordered so that the coefficients for terms of degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Chebyshev series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points in the Cartesion + product of `x` and `y`. See Also -------- @@ -1390,7 +1429,7 @@ def chebvander(x, deg) : def chebvander2d(x, y, deg) : - """Psuedo Vandermonde matrix of given degree. + """Pseudo Vandermonde matrix of given degree. Returns the pseudo Vandermonde matrix for 2D Chebyshev series in `x` and `y`. The sample point coordinates must all have the same shape @@ -1435,7 +1474,7 @@ def chebvander2d(x, y, deg) : def chebvander3d(x, y, z, deg) : - """Psuedo Vandermonde matrix of given degree. + """Pseudo Vandermonde matrix of given degree. Returns the pseudo Vandermonde matrix for 3D Chebyshev series in `x`, `y`, or `z`. The sample point coordinates must all have the same shape @@ -1636,18 +1675,18 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None): return c -def chebcompanion(cs): - """Return the scaled companion matrix of cs. +def chebcompanion(c): + """Return the scaled companion matrix of c. The basis polynomials are scaled so that the companion matrix is - symmetric when `cs` represents a single Chebyshev polynomial. This + symmetric when `c` represents a single Chebyshev polynomial. This provides better eigenvalue estimates than the unscaled case and in the single polynomial case the eigenvalues are guaranteed to be real if np.eigvalsh is used to obtain them. Parameters ---------- - cs : array_like + c : array_like 1-d array of Legendre series coefficients ordered from low to high degree. @@ -1657,56 +1696,60 @@ def chebcompanion(cs): Scaled companion matrix of dimensions (deg, deg). """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) < 2: + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) < 2: raise ValueError('Series must have maximum degree of at least 1.') - if len(cs) == 2: - return np.array(-cs[0]/cs[1]) + if len(c) == 2: + return np.array(-c[0]/c[1]) - n = len(cs) - 1 - mat = np.zeros((n, n), dtype=cs.dtype) + n = len(c) - 1 + mat = np.zeros((n, n), dtype=c.dtype) scl = np.array([1.] + [np.sqrt(.5)]*(n-1)) top = mat.reshape(-1)[1::n+1] bot = mat.reshape(-1)[n::n+1] top[0] = np.sqrt(.5) top[1:] = 1/2 bot[...] = top - mat[:,-1] -= (cs[:-1]/cs[-1])*(scl/scl[-1])*.5 + mat[:,-1] -= (c[:-1]/c[-1])*(scl/scl[-1])*.5 return mat -def chebroots(cs): +def chebroots(c): """ Compute the roots of a Chebyshev series. - Return the roots (a.k.a "zeros") of the C-series represented by `cs`, - which is the sequence of the C-series' coefficients from lowest order - "term" to highest, e.g., [1,2,3] represents the C-series - ``T_0 + 2*T_1 + 3*T_2``. + Return the roots (a.k.a. "zeros") of the polynomial + + .. math:: p(x) = \\sum_i c[i] * T_i(x). Parameters ---------- - cs : array_like - 1-d array of C-series coefficients ordered from low to high. + c : 1-D array_like + 1-D array of coefficients. Returns ------- out : ndarray - Array of the roots. If all the roots are real, then so is the - dtype of ``out``; otherwise, ``out``'s dtype is complex. + Array of the roots of the series. If all the roots are real, + then `out` is also real, otherwise it is complex. See Also -------- - polyroots + polyroots, legroots, lagroots, hermroots, hermeroots Notes ----- - Algorithm(s) used: + The root estimates are obtained as the eigenvalues of the companion + matrix, Roots far from the origin of the complex plane may have large + errors due to the numerical instability of the series for such + values. Roots with multiplicity greater than 1 will also show larger + errors as the value of the series near such points is relatively + insensitive to errors in the roots. Isolated roots near the origin can + be improved by a few iterations of Newton's method. - Remember: because the C-series basis set is different from the - "standard" basis set, the results of this function *may* not be what - one is expecting. + The Chebyshev series basis polynomials aren't powers of `x` so the + results of this function may seem unintuitive. Examples -------- @@ -1715,14 +1758,14 @@ def chebroots(cs): array([ -5.00000000e-01, 2.60860684e-17, 1.00000000e+00]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) < 2: - return np.array([], dtype=cs.dtype) - if len(cs) == 2: - return np.array([-cs[0]/cs[1]]) - - m = chebcompanion(cs) + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) < 2: + return np.array([], dtype=c.dtype) + if len(c) == 2: + return np.array([-c[0]/c[1]]) + + m = chebcompanion(c) r = la.eigvals(m) r.sort() return r diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py index 58cc655d7..b9862ad5a 100644 --- a/numpy/polynomial/hermite.py +++ b/numpy/polynomial/hermite.py @@ -94,7 +94,7 @@ def poly2herm(pol) : Returns ------- - cs : ndarray + c : ndarray 1-d array containing the coefficients of the equivalent Hermite series. @@ -122,7 +122,7 @@ def poly2herm(pol) : return res -def herm2poly(cs) : +def herm2poly(c) : """ Convert a Hermite series to a polynomial. @@ -133,7 +133,7 @@ def herm2poly(cs) : Parameters ---------- - cs : array_like + c : array_like 1-d array containing the Hermite series coefficients, ordered from lowest order term to highest. @@ -162,20 +162,20 @@ def herm2poly(cs) : """ from polynomial import polyadd, polysub, polymulx - [cs] = pu.as_series([cs]) - n = len(cs) + [c] = pu.as_series([c]) + n = len(c) if n == 1: - return cs + return c if n == 2: - cs[1] *= 2 - return cs + c[1] *= 2 + return c else: - c0 = cs[-2] - c1 = cs[-1] + c0 = c[-2] + c1 = c[-1] # i is the current degree of c1 for i in range(n - 1, 1, -1) : tmp = c0 - c0 = polysub(cs[i - 2], c1*(2*(i - 1))) + c0 = polysub(c[i - 2], c1*(2*(i - 1))) c1 = polyadd(tmp, polymulx(c1)*2) return polyadd(c0, polymulx(c1)*2) @@ -235,14 +235,24 @@ def hermline(off, scl) : def hermfromroots(roots) : """ - Generate a Hermite series with the given roots. + Generate a Hermite series with given roots. + + The function returns the coefficients of the polynomial - Return the array of coefficients for the P-series whose roots (a.k.a. - "zeros") are given by *roots*. The returned array of coefficients is - ordered from lowest order "term" to highest, and zeros of multiplicity - greater than one must be included in *roots* a number of times equal - to their multiplicity (e.g., if `2` is a root of multiplicity three, - then [2,2,2] must be in *roots*). + .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), + + in Hermite form, where the `r_n` are the roots specified in `roots`. + If a zero has multiplicity n, then it must appear in `roots` n times. + For instance, if 2 is a root of multiplicity three and 3 is a root of + multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The + roots can appear in any order. + + If the returned coefficients are `c`, then + + .. math:: p(x) = c_0 + c_1 * H_1(x) + ... + c_n * H_n(x) + + The coefficient of the last term is not generally 1 for monic + polynomials in Hermite form. Parameters ---------- @@ -252,28 +262,15 @@ def hermfromroots(roots) : Returns ------- out : ndarray - 1-d array of the Hermite series coefficients, ordered from low to - high. If all roots are real, ``out.dtype`` is a float type; - otherwise, ``out.dtype`` is a complex type, even if all the - coefficients in the result are real (see Examples below). + 1-D array of coefficients. If all roots are real then `out` is a + real array, if some of the roots are complex, then `out` is complex + even if all the coefficients in the result are real (see Examples + below). See Also -------- - polyfromroots, chebfromroots - - Notes - ----- - What is returned are the :math:`c_i` such that: - - .. math:: - - \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i]) - - where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Hermite - (basis) polynomial over the domain `[-1,1]`. Note that, unlike - `polyfromroots`, due to the nature of the Hermite basis set, the - above identity *does not* imply :math:`c_n = 1` identically (see - Examples). + polyfromroots, legfromroots, lagfromroots, chebfromroots, + hermefromroots. Examples -------- @@ -393,16 +390,16 @@ def hermsub(c1, c2): return pu.trimseq(ret) -def hermmulx(cs): +def hermmulx(c): """Multiply a Hermite series by x. - Multiply the Hermite series `cs` by x, where x is the independent + Multiply the Hermite series `c` by x, where x is the independent variable. Parameters ---------- - cs : array_like + c : array_like 1-d array of Hermite series coefficients ordered from low to high. @@ -427,18 +424,18 @@ def hermmulx(cs): array([ 2. , 6.5, 1. , 1.5]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) # The zero series needs special treatment - if len(cs) == 1 and cs[0] == 0: - return cs - - prd = np.empty(len(cs) + 1, dtype=cs.dtype) - prd[0] = cs[0]*0 - prd[1] = cs[0]/2 - for i in range(1, len(cs)): - prd[i + 1] = cs[i]/2 - prd[i - 1] += cs[i]*i + if len(c) == 1 and c[0] == 0: + return c + + prd = np.empty(len(c) + 1, dtype=c.dtype) + prd[0] = c[0]*0 + prd[1] = c[0]/2 + for i in range(1, len(c)): + prd[i + 1] = c[i]/2 + prd[i - 1] += c[i]*i return prd @@ -484,26 +481,26 @@ def hermmul(c1, c2): [c1, c2] = pu.as_series([c1, c2]) if len(c1) > len(c2): - cs = c2 + c = c2 xs = c1 else: - cs = c1 + c = c1 xs = c2 - if len(cs) == 1: - c0 = cs[0]*xs + if len(c) == 1: + c0 = c[0]*xs c1 = 0 - elif len(cs) == 2: - c0 = cs[0]*xs - c1 = cs[1]*xs + elif len(c) == 2: + c0 = c[0]*xs + c1 = c[1]*xs else : - nd = len(cs) - c0 = cs[-2]*xs - c1 = cs[-1]*xs - for i in range(3, len(cs) + 1) : + nd = len(c) + c0 = c[-2]*xs + c1 = c[-1]*xs + for i in range(3, len(c) + 1) : tmp = c0 nd = nd - 1 - c0 = hermsub(cs[-i]*xs, c1*(2*(nd - 1))) + c0 = hermsub(c[-i]*xs, c1*(2*(nd - 1))) c1 = hermadd(tmp, hermmulx(c1)*2) return hermadd(c0, hermmulx(c1)*2) @@ -575,16 +572,16 @@ def hermdiv(c1, c2): return quo, pu.trimseq(rem) -def hermpow(cs, pow, maxpower=16) : +def hermpow(c, pow, maxpower=16) : """Raise a Hermite series to a power. - Returns the Hermite series `cs` raised to the power `pow`. The - arguement `cs` is a sequence of coefficients ordered from low to high. + Returns the Hermite series `c` raised to the power `pow`. The + arguement `c` is a sequence of coefficients ordered from low to high. i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.`` Parameters ---------- - cs : array_like + c : array_like 1d array of Hermite series coefficients ordered from low to high. pow : integer @@ -609,23 +606,23 @@ def hermpow(cs, pow, maxpower=16) : array([ 81., 52., 82., 12., 9.]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) power = int(pow) if power != pow or power < 0 : raise ValueError("Power must be a non-negative integer.") elif maxpower is not None and power > maxpower : raise ValueError("Power is too large") elif power == 0 : - return np.array([1], dtype=cs.dtype) + return np.array([1], dtype=c.dtype) elif power == 1 : - return cs + return c else : # This can be made more efficient by using powers of two # in the usual way. - prd = cs + prd = c for i in range(2, power + 1) : - prd = hermmul(prd, cs) + prd = hermmul(prd, c) return prd @@ -775,8 +772,8 @@ def hermint(c, m=1, k=[], lbnd=0, scl=1, axis=0): Note that the result of each integration is *multiplied* by `scl`. Why is this important to note? Say one is making a linear change of variable :math:`u = ax + b` in an integral relative to `x`. Then - :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a` - - perhaps not what one would have first thought. + .. math::`dx = du/a`, so one will need to set `scl` equal to + :math:`1/a` - perhaps not what one would have first thought. Also note that, in general, the result of integrating a C-series needs to be "re-projected" onto the C-series basis set. Thus, typically, @@ -842,48 +839,48 @@ def hermint(c, m=1, k=[], lbnd=0, scl=1, axis=0): def hermval(x, c, tensor=True): """ - Evaluate a Hermite series. + Evaluate an Hermite series at points x. - If `c` is of length ``n + 1``, this function returns the value: + If `c` is of length `n + 1`, this function returns the value: - ``p(x) = c[0]*H_0(x) + c[1]*H_1(x) + ... + c[n]*H_n(x)`` + .. math:: p(x) = c_0 * H_0(x) + c_1 * H_1(x) + ... + c_n * H_n(x) - If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)`` - will have the same shape as `x`. If `x` is a algebra_like object that - supports multiplication and addition with itself and the values in `c`, - then an object of the same type is returned. + The parameter `x` is converted to an array only if it is a tuple or a + list, otherwise it is treated as a scalar. In either case, either `x` + or its elements must support multiplication and addition both with + themselves and with the elements of `c`. - In the case where c is multidimensional, the shape of the result - depends on the value of `tensor`. If tensor is true the shape of the - return will be ``c.shape[1:] + x.shape``, where the shape of a scalar - is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if - `x` is broadcast compatible with that. + If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If + `c` is multidimensional, then the shape of the result depends on the + value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + + x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that + scalars have shape (,). - If there are trailing zeros in the coefficients they still take part in - the evaluation, so they should be avoided if efficiency is a concern. + Trailing zeros in the coefficients will be used in the evaluation, so + they should be avoided if efficiency is a concern. Parameters ---------- - x : array_like, algebra_like - If x is a list or tuple, it is converted to an ndarray. Otherwise - it is left unchanged and if it isn't an ndarray it is treated as a - scalar. In either case, `x` or any element of an ndarray must - support addition and multiplication with itself and the elements of - `c`. + x : array_like, compatible object + If `x` is a list or tuple, it is converted to an ndarray, otherwise + it is left unchanged and treated as a scalar. In either case, `x` + or its elements must support addition and multiplication with + with themselves and with the elements of `c`. c : array_like Array of coefficients ordered so that the coefficients for terms of - degree n are contained in ``c[n]``. If `c` is multidimesional the - remaining indices enumerate multiple Hermite series. In the two + degree n are contained in c[n]. If `c` is multidimesional the + remaining indices enumerate multiple polynomials. In the two dimensional case the coefficients may be thought of as stored in the columns of `c`. tensor : boolean, optional - If true, the coefficient array shape is extended with ones on the - right, one for each dimension of `x`. Scalars are treated as having - dimension 0 for this action. The effect is that every column of - coefficients in `c` is evaluated for every value in `x`. If False, - the `x` are broadcast over the columns of `c` in the usual way. - This gives some flexibility in evaluations in the multidimensional - case. The default value it ``True``. + If True, the shape of the coefficient array is extended with ones + on the right, one for each dimension of `x`. Scalars have dimension 0 + for this action. The result is that every column of coefficients in + `c` is evaluated for every element of `x`. If False, `x` is broadcast + over the columns of `c` for the evaluation. This keyword is useful + when `c` is multidimensional. The default value is True. + + .. versionadded:: 1.7.0 Returns ------- @@ -938,31 +935,41 @@ def hermval(x, c, tensor=True): def hermval2d(x, y, c): """ - Evaluate 2D Hermite series at points (x,y). + Evaluate a 2-D Hermite series at points (x, y). This function returns the values: - ``p(x,y) = \sum_{i,j} c[i,j] * H_i(x) * H_j(y)`` + .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * H_i(x) * H_j(y) + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars and they + must have the same shape after conversion. In either case, either `x` + and `y` or their elements must support multiplication and addition both + with themselves and with the elements of `c`. + + If `c` is a 1-D array a one is implicitly appended to its shape to make + it 2-D. The shape of the result will be c.shape[2:] + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Hermite seres is evaluated at the points - ``(x,y)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray. - Otherwise it is left unchanged and if it isn't an ndarray it is - treated as a scalar. See `hermval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points `(x, y)`, + where `x` and `y` must have the same shape. If `x` or `y` is a list + or tuple, it is first converted to an ndarray, otherwise it is left + unchanged and if it isn't an ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of - coefficients. + Array of coefficients ordered so that the coefficient of the term + of multi-degree i,j is contained in ``c[i,j]``. If `c` has + dimension greater than two the remaining indices enumerate multiple + sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the two dimensional Hermite series at points formed - from pairs of corresponding values from `x` and `y`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points formed with + pairs of corresponding values from `x` and `y`. See Also -------- @@ -981,35 +988,45 @@ def hermval2d(x, y, c): def hermgrid2d(x, y, c): """ - Evaluate 2D Hermite series on the Cartesion product of x,y. + Evaluate a 2-D Hermite series on the Cartesion product of x and y. This function returns the values: - ``p(a,b) = \sum_{i,j} c[i,j] * H_i(a) * H_j(b)`` + .. math:: p(a,b) = \sum_{i,j} c_{i,j} * H_i(a) * H_j(b) - where the points ``(a,b)`` consist of all pairs of points formed by - taking ``a`` from `x` and ``b`` from `y`. The resulting points form a - grid with `x` in the first dimension and `y` in the second. + where the points `(a, b)` consist of all pairs formed by taking + `a` from `x` and `b` from `y`. The resulting points form a grid with + `x` in the first dimension and `y` in the second. + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars. In either + case, either `x` and `y` or their elements must support multiplication + and addition both with themselves and with the elements of `c`. + + If `c` has fewer than two dimensions, ones are implicitly appended to + its shape to make it 2-D. The shape of the result will be c.shape[2:] + + x.shape. + + .. versionadded:: 1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Hermite series is evaluated at the points in - the Cartesian product of `x` and `y`. If `x` or `y` is a list or - tuple, it is first converted to an ndarray, Otherwise it is left - unchanged and if it isn't an ndarray it is treated as a scalar. See - `hermval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points in the + Cartesian product of `x` and `y`. If `x` or `y` is a list or + tuple, it is first converted to an ndarray, otherwise it is left + unchanged and, if it isn't an ndarray, it is treated as a scalar. c : array_like Array of coefficients ordered so that the coefficients for terms of degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the two dimensional Hermite series at points in the - Cartesion product of `x` and `y`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points in the Cartesion + product of `x` and `y`. See Also -------- @@ -1023,32 +1040,43 @@ def hermgrid2d(x, y, c): def hermval3d(x, y, z, c): """ - Evaluate 3D Hermite series at points (x,y,z). + Evaluate a 3-D Hermite series at points (x, y, z). This function returns the values: - ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * H_i(x) * H_j(y) * H_k(z)`` + .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * H_i(x) * H_j(y) * H_k(z) + + The parameters `x`, `y`, and `z` are converted to arrays only if + they are tuples or a lists, otherwise they are treated as a scalars and + they must have the same shape after conversion. In either case, either + `x`, `y`, and `z` or their elements must support multiplication and + addition both with themselves and with the elements of `c`. + + If `c` has fewer than 3 dimensions, ones are implicitly appended to its + shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Hermite seres is evaluated at the points - ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape. If + x, y, z : array_like, compatible object + The three dimensional series is evaluated at the points + `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If any of `x`, `y`, or `z` is a list or tuple, it is first converted - to an ndarray. Otherwise it is left unchanged and if it isn't an - ndarray it is treated as a scalar. See `hermval` for explanation of - algebra_like. + to an ndarray, otherwise it is left unchanged and if it isn't an + ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + Array of coefficients ordered so that the coefficient of the term of + multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension + greater than 3 the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Hermite series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the multidimension polynomial on points formed with + triples of corresponding values from `x`, `y`, and `z`. See Also -------- @@ -1068,37 +1096,48 @@ def hermval3d(x, y, z, c): def hermgrid3d(x, y, z, c): """ - Evaluate 3D Hermite series on the Cartesian product of x,y,z. + Evaluate a 3-D Hermite series on the Cartesian product of x, y, and z. This function returns the values: - ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * H_i(a) * H_j(b) * H_k(c)`` + .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * H_i(a) * H_j(b) * H_k(c) + + where the points `(a, b, c)` consist of all triples formed by taking + `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form + a grid with `x` in the first dimension, `y` in the second, and `z` in + the third. + + The parameters `x`, `y`, and `z` are converted to arrays only if they + are tuples or a lists, otherwise they are treated as a scalars. In + either case, either `x`, `y`, and `z` or their elements must support + multiplication and addition both with themselves and with the elements + of `c`. + + If `c` has fewer than three dimensions, ones are implicitly appended to + its shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape + yshape + z.shape. - where the points ``(a,b,c)`` consist of all triples formed by taking - ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting - points form a grid with `x` in the first dimension, `y` in the second, - and `z` in the third. + .. versionadded:: 1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Hermite seres is evaluated at the points - in the cartesian product of `x`, `y`, and `z` - ``(x,y,z)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray, - otherwise it is left unchanged and treated as a scalar. See - `hermval` for explanation of algebra_like. + x, y, z : array_like, compatible objects + The three dimensional series is evaluated at the points in the + Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a + list or tuple, it is first converted to an ndarray, otherwise it is + left unchanged and, if it isn't an ndarray, it is treated as a + scalar. c : array_like Array of coefficients ordered so that the coefficients for terms of degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Hermite series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points in the Cartesion + product of `x` and `y`. See Also -------- @@ -1208,7 +1247,7 @@ def hermvander2d(x, y, deg) : def hermvander3d(x, y, z, deg) : - """Psuedo Vandermonde matrix of given degree. + """Pseudo Vandermonde matrix of given degree. Returns the pseudo Vandermonde matrix for 3D Hermite series in `x`, `y`, or `z`. The sample point coordinates must all have the same shape @@ -1415,18 +1454,18 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None): return c -def hermcompanion(cs): - """Return the scaled companion matrix of cs. +def hermcompanion(c): + """Return the scaled companion matrix of c. The basis polynomials are scaled so that the companion matrix is - symmetric when `cs` represents a single Hermite polynomial. This + symmetric when `c` represents a single Hermite polynomial. This provides better eigenvalue estimates than the unscaled case and in the single polynomial case the eigenvalues are guaranteed to be real if `numpy.linalg.eigvalsh` is used to obtain them. Parameters ---------- - cs : array_like + c : array_like 1-d array of Legendre series coefficients ordered from low to high degree. @@ -1437,56 +1476,60 @@ def hermcompanion(cs): """ accprod = np.multiply.accumulate - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) < 2: + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) < 2: raise ValueError('Series must have maximum degree of at least 1.') - if len(cs) == 2: - return np.array(-.5*cs[0]/cs[1]) + if len(c) == 2: + return np.array(-.5*c[0]/c[1]) - n = len(cs) - 1 - mat = np.zeros((n, n), dtype=cs.dtype) + n = len(c) - 1 + mat = np.zeros((n, n), dtype=c.dtype) scl = np.hstack((1., np.sqrt(2.*np.arange(1,n)))) scl = np.multiply.accumulate(scl) top = mat.reshape(-1)[1::n+1] bot = mat.reshape(-1)[n::n+1] top[...] = np.sqrt(.5*np.arange(1,n)) bot[...] = top - mat[:,-1] -= (cs[:-1]/cs[-1])*(scl/scl[-1])*.5 + mat[:,-1] -= (c[:-1]/c[-1])*(scl/scl[-1])*.5 return mat -def hermroots(cs): +def hermroots(c): """ Compute the roots of a Hermite series. - Return the roots (a.k.a "zeros") of the Hermite series represented by - `cs`, which is the sequence of coefficients from lowest order "term" - to highest, e.g., [1,2,3] is the series ``L_0 + 2*L_1 + 3*L_2``. + Return the roots (a.k.a. "zeros") of the polynomial + + .. math:: p(x) = \\sum_i c[i] * H_i(x). Parameters ---------- - cs : array_like - 1-d array of Hermite series coefficients ordered from low to high. + c : 1-D array_like + 1-D array of coefficients. Returns ------- out : ndarray - Array of the roots. If all the roots are real, then so is the - dtype of ``out``; otherwise, ``out``'s dtype is complex. + Array of the roots of the series. If all the roots are real, + then `out` is also real, otherwise it is complex. See Also -------- - polyroots - chebroots + polyroots, legroots, lagroots, chebroots, hermeroots Notes ----- - Algorithm(s) used: + The root estimates are obtained as the eigenvalues of the companion + matrix, Roots far from the origin of the complex plane may have large + errors due to the numerical instability of the series for such + values. Roots with multiplicity greater than 1 will also show larger + errors as the value of the series near such points is relatively + insensitive to errors in the roots. Isolated roots near the origin can + be improved by a few iterations of Newton's method. - Remember: because the Hermite series basis set is different from the - "standard" basis set, the results of this function *may* not be what - one is expecting. + The Hermite series basis polynomials aren't powers of `x` so the + results of this function may seem unintuitive. Examples -------- @@ -1498,14 +1541,14 @@ def hermroots(cs): array([ -1.00000000e+00, -1.38777878e-17, 1.00000000e+00]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) <= 1 : - return np.array([], dtype=cs.dtype) - if len(cs) == 2 : - return np.array([-.5*cs[0]/cs[1]]) - - m = hermcompanion(cs) + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) <= 1 : + return np.array([], dtype=c.dtype) + if len(c) == 2 : + return np.array([-.5*c[0]/c[1]]) + + m = hermcompanion(c) r = la.eigvals(m) r.sort() return r diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py index c8d641734..4f39827f9 100644 --- a/numpy/polynomial/hermite_e.py +++ b/numpy/polynomial/hermite_e.py @@ -94,7 +94,7 @@ def poly2herme(pol) : Returns ------- - cs : ndarray + c : ndarray 1-d array containing the coefficients of the equivalent Hermite series. @@ -122,7 +122,7 @@ def poly2herme(pol) : return res -def herme2poly(cs) : +def herme2poly(c) : """ Convert a Hermite series to a polynomial. @@ -133,7 +133,7 @@ def herme2poly(cs) : Parameters ---------- - cs : array_like + c : array_like 1-d array containing the Hermite series coefficients, ordered from lowest order term to highest. @@ -162,19 +162,19 @@ def herme2poly(cs) : """ from polynomial import polyadd, polysub, polymulx - [cs] = pu.as_series([cs]) - n = len(cs) + [c] = pu.as_series([c]) + n = len(c) if n == 1: - return cs + return c if n == 2: - return cs + return c else: - c0 = cs[-2] - c1 = cs[-1] + c0 = c[-2] + c1 = c[-1] # i is the current degree of c1 for i in range(n - 1, 1, -1) : tmp = c0 - c0 = polysub(cs[i - 2], c1*(i - 1)) + c0 = polysub(c[i - 2], c1*(i - 1)) c1 = polyadd(tmp, polymulx(c1)) return polyadd(c0, polymulx(c1)) @@ -235,14 +235,24 @@ def hermeline(off, scl) : def hermefromroots(roots) : """ - Generate a Hermite series with the given roots. + Generate a HermiteE series with given roots. + + The function returns the coefficients of the polynomial - Return the array of coefficients for the P-series whose roots (a.k.a. - "zeros") are given by *roots*. The returned array of coefficients is - ordered from lowest order "term" to highest, and zeros of multiplicity - greater than one must be included in *roots* a number of times equal - to their multiplicity (e.g., if `2` is a root of multiplicity three, - then [2,2,2] must be in *roots*). + .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), + + in HermiteE form, where the `r_n` are the roots specified in `roots`. + If a zero has multiplicity n, then it must appear in `roots` n times. + For instance, if 2 is a root of multiplicity three and 3 is a root of + multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The + roots can appear in any order. + + If the returned coefficients are `c`, then + + .. math:: p(x) = c_0 + c_1 * He_1(x) + ... + c_n * He_n(x) + + The coefficient of the last term is not generally 1 for monic + polynomials in HermiteE form. Parameters ---------- @@ -252,28 +262,15 @@ def hermefromroots(roots) : Returns ------- out : ndarray - 1-d array of the Hermite series coefficients, ordered from low to - high. If all roots are real, ``out.dtype`` is a float type; - otherwise, ``out.dtype`` is a complex type, even if all the - coefficients in the result are real (see Examples below). + 1-D array of coefficients. If all roots are real then `out` is a + real array, if some of the roots are complex, then `out` is complex + even if all the coefficients in the result are real (see Examples + below). See Also -------- - polyfromroots, chebfromroots - - Notes - ----- - What is returned are the :math:`c_i` such that: - - .. math:: - - \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i]) - - where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Hermite - (basis) polynomial over the domain `[-1,1]`. Note that, unlike - `polyfromroots`, due to the nature of the Hermite basis set, the - above identity *does not* imply :math:`c_n = 1` identically (see - Examples). + polyfromroots, legfromroots, lagfromroots, hermfromroots, + chebfromroots. Examples -------- @@ -393,16 +390,16 @@ def hermesub(c1, c2): return pu.trimseq(ret) -def hermemulx(cs): +def hermemulx(c): """Multiply a Hermite series by x. - Multiply the Hermite series `cs` by x, where x is the independent + Multiply the Hermite series `c` by x, where x is the independent variable. Parameters ---------- - cs : array_like + c : array_like 1-d array of Hermite series coefficients ordered from low to high. @@ -427,18 +424,18 @@ def hermemulx(cs): array([ 2., 7., 2., 3.]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) # The zero series needs special treatment - if len(cs) == 1 and cs[0] == 0: - return cs - - prd = np.empty(len(cs) + 1, dtype=cs.dtype) - prd[0] = cs[0]*0 - prd[1] = cs[0] - for i in range(1, len(cs)): - prd[i + 1] = cs[i] - prd[i - 1] += cs[i]*i + if len(c) == 1 and c[0] == 0: + return c + + prd = np.empty(len(c) + 1, dtype=c.dtype) + prd[0] = c[0]*0 + prd[1] = c[0] + for i in range(1, len(c)): + prd[i + 1] = c[i] + prd[i - 1] += c[i]*i return prd @@ -484,26 +481,26 @@ def hermemul(c1, c2): [c1, c2] = pu.as_series([c1, c2]) if len(c1) > len(c2): - cs = c2 + c = c2 xs = c1 else: - cs = c1 + c = c1 xs = c2 - if len(cs) == 1: - c0 = cs[0]*xs + if len(c) == 1: + c0 = c[0]*xs c1 = 0 - elif len(cs) == 2: - c0 = cs[0]*xs - c1 = cs[1]*xs + elif len(c) == 2: + c0 = c[0]*xs + c1 = c[1]*xs else : - nd = len(cs) - c0 = cs[-2]*xs - c1 = cs[-1]*xs - for i in range(3, len(cs) + 1) : + nd = len(c) + c0 = c[-2]*xs + c1 = c[-1]*xs + for i in range(3, len(c) + 1) : tmp = c0 nd = nd - 1 - c0 = hermesub(cs[-i]*xs, c1*(nd - 1)) + c0 = hermesub(c[-i]*xs, c1*(nd - 1)) c1 = hermeadd(tmp, hermemulx(c1)) return hermeadd(c0, hermemulx(c1)) @@ -573,16 +570,16 @@ def hermediv(c1, c2): return quo, pu.trimseq(rem) -def hermepow(cs, pow, maxpower=16) : +def hermepow(c, pow, maxpower=16) : """Raise a Hermite series to a power. - Returns the Hermite series `cs` raised to the power `pow`. The - arguement `cs` is a sequence of coefficients ordered from low to high. + Returns the Hermite series `c` raised to the power `pow`. The + arguement `c` is a sequence of coefficients ordered from low to high. i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.`` Parameters ---------- - cs : array_like + c : array_like 1d array of Hermite series coefficients ordered from low to high. pow : integer @@ -607,23 +604,23 @@ def hermepow(cs, pow, maxpower=16) : array([ 23., 28., 46., 12., 9.]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) power = int(pow) if power != pow or power < 0 : raise ValueError("Power must be a non-negative integer.") elif maxpower is not None and power > maxpower : raise ValueError("Power is too large") elif power == 0 : - return np.array([1], dtype=cs.dtype) + return np.array([1], dtype=c.dtype) elif power == 1 : - return cs + return c else : # This can be made more efficient by using powers of two # in the usual way. - prd = cs + prd = c for i in range(2, power + 1) : - prd = hermemul(prd, cs) + prd = hermemul(prd, c) return prd @@ -773,8 +770,8 @@ def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0): Note that the result of each integration is *multiplied* by `scl`. Why is this important to note? Say one is making a linear change of variable :math:`u = ax + b` in an integral relative to `x`. Then - :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a` - - perhaps not what one would have first thought. + .. math::`dx = du/a`, so one will need to set `scl` equal to + :math:`1/a` - perhaps not what one would have first thought. Also note that, in general, the result of integrating a C-series needs to be "re-projected" onto the C-series basis set. Thus, typically, @@ -839,48 +836,49 @@ def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0): def hermeval(x, c, tensor=True): - """ Evaluate a Hermite_e series. + """ + Evaluate an HermiteE series at points x. - If `c` is of length ``n + 1``, this function returns the value: + If `c` is of length `n + 1`, this function returns the value: - ``p(x) = c[0]*He_0(x) + c[1]*He_1(x) + ... + c[n]*He_n(x)`` + .. math:: p(x) = c_0 * He_0(x) + c_1 * He_1(x) + ... + c_n * He_n(x) - If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)`` - will have the same shape as `x`. If `x` is a algebra_like object that - supports multiplication and addition with itself and the values in `c`, - then an object of the same type is returned. + The parameter `x` is converted to an array only if it is a tuple or a + list, otherwise it is treated as a scalar. In either case, either `x` + or its elements must support multiplication and addition both with + themselves and with the elements of `c`. - In the case where c is multidimensional, the shape of the result - depends on the value of `tensor`. If tensor is true the shape of the - return will be ``c.shape[1:] + x.shape``, where the shape of a scalar - is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if - `x` is broadcast compatible with that. + If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If + `c` is multidimensional, then the shape of the result depends on the + value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + + x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that + scalars have shape (,). - If there are trailing zeros in the coefficients they still take part in - the evaluation, so they should be avoided if efficiency is a concern. + Trailing zeros in the coefficients will be used in the evaluation, so + they should be avoided if efficiency is a concern. Parameters ---------- - x : array_like, algebra_like - If x is a list or tuple, it is converted to an ndarray. Otherwise - it is left unchanged and if it isn't an ndarray it is treated as a - scalar. In either case, `x` or any element of an ndarray must - support addition and multiplication with itself and the elements of - `c`. + x : array_like, compatible object + If `x` is a list or tuple, it is converted to an ndarray, otherwise + it is left unchanged and treated as a scalar. In either case, `x` + or its elements must support addition and multiplication with + with themselves and with the elements of `c`. c : array_like Array of coefficients ordered so that the coefficients for terms of - degree n are contained in ``c[n]``. If `c` is multidimesional the - remaining indices enumerate multiple Hermite_e series. In the two + degree n are contained in c[n]. If `c` is multidimesional the + remaining indices enumerate multiple polynomials. In the two dimensional case the coefficients may be thought of as stored in the columns of `c`. tensor : boolean, optional - If true, the coefficient array shape is extended with ones on the - right, one for each dimension of `x`. Scalars are treated as having - dimension 0 for this action. The effect is that every column of - coefficients in `c` is evaluated for every value in `x`. If False, - the `x` are broadcast over the columns of `c` in the usual way. - This gives some flexibility in evaluations in the multidimensional - case. The default value it ``True``. + If True, the shape of the coefficient array is extended with ones + on the right, one for each dimension of `x`. Scalars have dimension 0 + for this action. The result is that every column of coefficients in + `c` is evaluated for every element of `x`. If False, `x` is broadcast + over the columns of `c` for the evaluation. This keyword is useful + when `c` is multidimensional. The default value is True. + + .. versionadded:: 1.7.0 Returns ------- @@ -934,31 +932,41 @@ def hermeval(x, c, tensor=True): def hermeval2d(x, y, c): """ - Evaluate 2D Hermite_e series at points (x,y). + Evaluate a 2-D HermiteE series at points (x, y). This function returns the values: - ``p(x,y) = \sum_{i,j} c[i,j] * He_i(x) * He_j(y)`` + .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * He_i(x) * He_j(y) + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars and they + must have the same shape after conversion. In either case, either `x` + and `y` or their elements must support multiplication and addition both + with themselves and with the elements of `c`. + + If `c` is a 1-D array a one is implicitly appended to its shape to make + it 2-D. The shape of the result will be c.shape[2:] + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Hermite_e seres is evaluated at the points - ``(x,y)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray. - Otherwise it is left unchanged and if it isn't an ndarray it is - treated as a scalar. See `hermeval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points `(x, y)`, + where `x` and `y` must have the same shape. If `x` or `y` is a list + or tuple, it is first converted to an ndarray, otherwise it is left + unchanged and if it isn't an ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of - coefficients. + Array of coefficients ordered so that the coefficient of the term + of multi-degree i,j is contained in ``c[i,j]``. If `c` has + dimension greater than two the remaining indices enumerate multiple + sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the two dimensional Hermite_e series at points formed - from pairs of corresponding values from `x` and `y`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points formed with + pairs of corresponding values from `x` and `y`. See Also -------- @@ -977,35 +985,45 @@ def hermeval2d(x, y, c): def hermegrid2d(x, y, c): """ - Evaluate 2D Hermite_e series on the Cartesion product of x,y. + Evaluate a 2-D HermiteE series on the Cartesion product of x and y. This function returns the values: - ``p(a,b) = \sum_{i,j} c[i,j] * He_i(a) * He_j(b)`` + .. math:: p(a,b) = \sum_{i,j} c_{i,j} * H_i(a) * H_j(b) - where the points ``(a,b)`` consist of all pairs of points formed by - taking ``a`` from `x` and ``b`` from `y`. The resulting points form a - grid with `x` in the first dimension and `y` in the second. + where the points `(a, b)` consist of all pairs formed by taking + `a` from `x` and `b` from `y`. The resulting points form a grid with + `x` in the first dimension and `y` in the second. + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars. In either + case, either `x` and `y` or their elements must support multiplication + and addition both with themselves and with the elements of `c`. + + If `c` has fewer than two dimensions, ones are implicitly appended to + its shape to make it 2-D. The shape of the result will be c.shape[2:] + + x.shape. + + .. versionadded:: 1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Hermite_e series is evaluated at the points in - the Cartesian product of `x` and `y`. If `x` or `y` is a list or - tuple, it is first converted to an ndarray, Otherwise it is left - unchanged and if it isn't an ndarray it is treated as a scalar. See - `hermeval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points in the + Cartesian product of `x` and `y`. If `x` or `y` is a list or + tuple, it is first converted to an ndarray, otherwise it is left + unchanged and, if it isn't an ndarray, it is treated as a scalar. c : array_like Array of coefficients ordered so that the coefficients for terms of degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the two dimensional Hermite_e series at points in the - Cartesion product of `x` and `y`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points in the Cartesion + product of `x` and `y`. See Also -------- @@ -1019,32 +1037,43 @@ def hermegrid2d(x, y, c): def hermeval3d(x, y, z, c): """ - Evaluate 3D Hermite_e series at points (x,y,z). + Evaluate a 3-D Hermite_e series at points (x, y, z). This function returns the values: - ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * He_i(x) * He_j(y) * He_k(z)`` + .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * He_i(x) * He_j(y) * He_k(z) + + The parameters `x`, `y`, and `z` are converted to arrays only if + they are tuples or a lists, otherwise they are treated as a scalars and + they must have the same shape after conversion. In either case, either + `x`, `y`, and `z` or their elements must support multiplication and + addition both with themselves and with the elements of `c`. + + If `c` has fewer than 3 dimensions, ones are implicitly appended to its + shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Hermite_e seres is evaluated at the points - ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape. If + x, y, z : array_like, compatible object + The three dimensional series is evaluated at the points + `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If any of `x`, `y`, or `z` is a list or tuple, it is first converted - to an ndarray. Otherwise it is left unchanged and if it isn't an - ndarray it is treated as a scalar. See `hermeval` for explanation of - algebra_like. + to an ndarray, otherwise it is left unchanged and if it isn't an + ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + Array of coefficients ordered so that the coefficient of the term of + multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension + greater than 3 the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Hermite_e series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the multidimension polynomial on points formed with + triples of corresponding values from `x`, `y`, and `z`. See Also -------- @@ -1064,37 +1093,48 @@ def hermeval3d(x, y, z, c): def hermegrid3d(x, y, z, c): """ - Evaluate 3D Hermite_e series on the Cartesian product of x,y,z. + Evaluate a 3-D HermiteE series on the Cartesian product of x, y, and z. This function returns the values: - ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * He_i(a) * He_j(b) * He_k(c)`` + .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * He_i(a) * He_j(b) * He_k(c) + + where the points `(a, b, c)` consist of all triples formed by taking + `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form + a grid with `x` in the first dimension, `y` in the second, and `z` in + the third. + + The parameters `x`, `y`, and `z` are converted to arrays only if they + are tuples or a lists, otherwise they are treated as a scalars. In + either case, either `x`, `y`, and `z` or their elements must support + multiplication and addition both with themselves and with the elements + of `c`. + + If `c` has fewer than three dimensions, ones are implicitly appended to + its shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape + yshape + z.shape. - where the points ``(a,b,c)`` consist of all triples formed by taking - ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting - points form a grid with `x` in the first dimension, `y` in the second, - and `z` in the third. + .. versionadded:: 1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Hermite_e seres is evaluated at the points - in the cartesian product of `x`, `y`, and `z` - ``(x,y,z)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray, - otherwise it is left unchanged and treated as a scalar. See - `hermeval` for explanation of algebra_like. + x, y, z : array_like, compatible objects + The three dimensional series is evaluated at the points in the + Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a + list or tuple, it is first converted to an ndarray, otherwise it is + left unchanged and, if it isn't an ndarray, it is treated as a + scalar. c : array_like Array of coefficients ordered so that the coefficients for terms of degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Hermite_e series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points in the Cartesion + product of `x` and `y`. See Also -------- @@ -1203,7 +1243,7 @@ def hermevander2d(x, y, deg) : def hermevander3d(x, y, z, deg) : - """Psuedo Vandermonde matrix of given degree. + """Pseudo Vandermonde matrix of given degree. Returns the pseudo Vandermonde matrix for 3D Hermite_e series in `x`, `y`, or `z`. The sample point coordinates must all have the same shape @@ -1410,18 +1450,18 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None): return c -def hermecompanion(cs): - """Return the scaled companion matrix of cs. +def hermecompanion(c): + """Return the scaled companion matrix of c. The basis polynomials are scaled so that the companion matrix is - symmetric when `cs` represents a single HermiteE polynomial. This + symmetric when `c` represents a single HermiteE polynomial. This provides better eigenvalue estimates than the unscaled case and in the single polynomial case the eigenvalues are guaranteed to be real if `numpy.linalg.eigvalsh` is used to obtain them. Parameters ---------- - cs : array_like + c : array_like 1-d array of Legendre series coefficients ordered from low to high degree. @@ -1432,56 +1472,60 @@ def hermecompanion(cs): """ accprod = np.multiply.accumulate - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) < 2: + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) < 2: raise ValueError('Series must have maximum degree of at least 1.') - if len(cs) == 2: - return np.array(-cs[0]/cs[1]) + if len(c) == 2: + return np.array(-c[0]/c[1]) - n = len(cs) - 1 - mat = np.zeros((n, n), dtype=cs.dtype) + n = len(c) - 1 + mat = np.zeros((n, n), dtype=c.dtype) scl = np.hstack((1., np.sqrt(np.arange(1,n)))) scl = np.multiply.accumulate(scl) top = mat.reshape(-1)[1::n+1] bot = mat.reshape(-1)[n::n+1] top[...] = np.sqrt(np.arange(1,n)) bot[...] = top - mat[:,-1] -= (cs[:-1]/cs[-1])*(scl/scl[-1]) + mat[:,-1] -= (c[:-1]/c[-1])*(scl/scl[-1]) return mat -def hermeroots(cs): +def hermeroots(c): """ - Compute the roots of a Hermite series. + Compute the roots of a HermiteE series. - Return the roots (a.k.a "zeros") of the HermiteE series represented by - `cs`, which is the sequence of coefficients from lowest order "term" - to highest, e.g., [1,2,3] is the series ``L_0 + 2*L_1 + 3*L_2``. + Return the roots (a.k.a. "zeros") of the polynomial + + .. math:: p(x) = \\sum_i c[i] * He_i(x). Parameters ---------- - cs : array_like - 1-d array of HermiteE series coefficients ordered from low to high. + c : 1-D array_like + 1-D array of coefficients. Returns ------- out : ndarray - Array of the roots. If all the roots are real, then so is the - dtype of ``out``; otherwise, ``out``'s dtype is complex. + Array of the roots of the series. If all the roots are real, + then `out` is also real, otherwise it is complex. See Also -------- - polyroots - chebroots + polyroots, legroots, lagroots, hermroots, chebroots Notes ----- - Algorithm(s) used: + The root estimates are obtained as the eigenvalues of the companion + matrix, Roots far from the origin of the complex plane may have large + errors due to the numerical instability of the series for such + values. Roots with multiplicity greater than 1 will also show larger + errors as the value of the series near such points is relatively + insensitive to errors in the roots. Isolated roots near the origin can + be improved by a few iterations of Newton's method. - Remember: because the Hermite series basis set is different from the - "standard" basis set, the results of this function *may* not be what - one is expecting. + The HermiteE series basis polynomials aren't powers of `x` so the + results of this function may seem unintuitive. Examples -------- @@ -1493,14 +1537,14 @@ def hermeroots(cs): array([-1., 0., 1.]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) <= 1 : - return np.array([], dtype=cs.dtype) - if len(cs) == 2 : - return np.array([-cs[0]/cs[1]]) - - m = hermecompanion(cs) + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) <= 1 : + return np.array([], dtype=c.dtype) + if len(c) == 2 : + return np.array([-c[0]/c[1]]) + + m = hermecompanion(c) r = la.eigvals(m) r.sort() return r diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py index 872829fbc..15ea8d870 100644 --- a/numpy/polynomial/laguerre.py +++ b/numpy/polynomial/laguerre.py @@ -93,7 +93,7 @@ def poly2lag(pol) : Returns ------- - cs : ndarray + c : ndarray 1-d array containing the coefficients of the equivalent Laguerre series. @@ -121,7 +121,7 @@ def poly2lag(pol) : return res -def lag2poly(cs) : +def lag2poly(c) : """ Convert a Laguerre series to a polynomial. @@ -132,7 +132,7 @@ def lag2poly(cs) : Parameters ---------- - cs : array_like + c : array_like 1-d array containing the Laguerre series coefficients, ordered from lowest order term to highest. @@ -161,17 +161,17 @@ def lag2poly(cs) : """ from polynomial import polyadd, polysub, polymulx - [cs] = pu.as_series([cs]) - n = len(cs) + [c] = pu.as_series([c]) + n = len(c) if n == 1: - return cs + return c else: - c0 = cs[-2] - c1 = cs[-1] + c0 = c[-2] + c1 = c[-1] # i is the current degree of c1 for i in range(n - 1, 1, -1): tmp = c0 - c0 = polysub(cs[i - 2], (c1*(i - 1))/i) + c0 = polysub(c[i - 2], (c1*(i - 1))/i) c1 = polyadd(tmp, polysub((2*i - 1)*c1, polymulx(c1))/i) return polyadd(c0, polysub(c1, polymulx(c1))) @@ -231,14 +231,24 @@ def lagline(off, scl) : def lagfromroots(roots) : """ - Generate a Laguerre series with the given roots. + Generate a Laguerre series with given roots. + + The function returns the coefficients of the polynomial + + .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), - Return the array of coefficients for the P-series whose roots (a.k.a. - "zeros") are given by *roots*. The returned array of coefficients is - ordered from lowest order "term" to highest, and zeros of multiplicity - greater than one must be included in *roots* a number of times equal - to their multiplicity (e.g., if `2` is a root of multiplicity three, - then [2,2,2] must be in *roots*). + in Laguerre form, where the `r_n` are the roots specified in `roots`. + If a zero has multiplicity n, then it must appear in `roots` n times. + For instance, if 2 is a root of multiplicity three and 3 is a root of + multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The + roots can appear in any order. + + If the returned coefficients are `c`, then + + .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x) + + The coefficient of the last term is not generally 1 for monic + polynomials in Laguerre form. Parameters ---------- @@ -248,28 +258,15 @@ def lagfromroots(roots) : Returns ------- out : ndarray - 1-d array of the Laguerre series coefficients, ordered from low to - high. If all roots are real, ``out.dtype`` is a float type; - otherwise, ``out.dtype`` is a complex type, even if all the - coefficients in the result are real (see Examples below). + 1-D array of coefficients. If all roots are real then `out` is a + real array, if some of the roots are complex, then `out` is complex + even if all the coefficients in the result are real (see Examples + below). See Also -------- - polyfromroots, chebfromroots - - Notes - ----- - What is returned are the :math:`c_i` such that: - - .. math:: - - \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i]) - - where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Laguerre - (basis) polynomial over the domain `[-1,1]`. Note that, unlike - `polyfromroots`, due to the nature of the Laguerre basis set, the - above identity *does not* imply :math:`c_n = 1` identically (see - Examples). + polyfromroots, legfromroots, chebfromroots, hermfromroots, + hermefromroots. Examples -------- @@ -390,16 +387,16 @@ def lagsub(c1, c2): return pu.trimseq(ret) -def lagmulx(cs): +def lagmulx(c): """Multiply a Laguerre series by x. - Multiply the Laguerre series `cs` by x, where x is the independent + Multiply the Laguerre series `c` by x, where x is the independent variable. Parameters ---------- - cs : array_like + c : array_like 1-d array of Laguerre series coefficients ordered from low to high. @@ -424,19 +421,19 @@ def lagmulx(cs): array([ -1., -1., 11., -9.]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) # The zero series needs special treatment - if len(cs) == 1 and cs[0] == 0: - return cs - - prd = np.empty(len(cs) + 1, dtype=cs.dtype) - prd[0] = cs[0] - prd[1] = -cs[0] - for i in range(1, len(cs)): - prd[i + 1] = -cs[i]*(i + 1) - prd[i] += cs[i]*(2*i + 1) - prd[i - 1] -= cs[i]*i + if len(c) == 1 and c[0] == 0: + return c + + prd = np.empty(len(c) + 1, dtype=c.dtype) + prd[0] = c[0] + prd[1] = -c[0] + for i in range(1, len(c)): + prd[i + 1] = -c[i]*(i + 1) + prd[i] += c[i]*(2*i + 1) + prd[i - 1] -= c[i]*i return prd @@ -482,26 +479,26 @@ def lagmul(c1, c2): [c1, c2] = pu.as_series([c1, c2]) if len(c1) > len(c2): - cs = c2 + c = c2 xs = c1 else: - cs = c1 + c = c1 xs = c2 - if len(cs) == 1: - c0 = cs[0]*xs + if len(c) == 1: + c0 = c[0]*xs c1 = 0 - elif len(cs) == 2: - c0 = cs[0]*xs - c1 = cs[1]*xs + elif len(c) == 2: + c0 = c[0]*xs + c1 = c[1]*xs else : - nd = len(cs) - c0 = cs[-2]*xs - c1 = cs[-1]*xs - for i in range(3, len(cs) + 1) : + nd = len(c) + c0 = c[-2]*xs + c1 = c[-1]*xs + for i in range(3, len(c) + 1) : tmp = c0 nd = nd - 1 - c0 = lagsub(cs[-i]*xs, (c1*(nd - 1))/nd) + c0 = lagsub(c[-i]*xs, (c1*(nd - 1))/nd) c1 = lagadd(tmp, lagsub((2*nd - 1)*c1, lagmulx(c1))/nd) return lagadd(c0, lagsub(c1, lagmulx(c1))) @@ -571,16 +568,16 @@ def lagdiv(c1, c2): return quo, pu.trimseq(rem) -def lagpow(cs, pow, maxpower=16) : +def lagpow(c, pow, maxpower=16) : """Raise a Laguerre series to a power. - Returns the Laguerre series `cs` raised to the power `pow`. The - arguement `cs` is a sequence of coefficients ordered from low to high. + Returns the Laguerre series `c` raised to the power `pow`. The + arguement `c` is a sequence of coefficients ordered from low to high. i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.`` Parameters ---------- - cs : array_like + c : array_like 1d array of Laguerre series coefficients ordered from low to high. pow : integer @@ -605,23 +602,23 @@ def lagpow(cs, pow, maxpower=16) : array([ 14., -16., 56., -72., 54.]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) power = int(pow) if power != pow or power < 0 : raise ValueError("Power must be a non-negative integer.") elif maxpower is not None and power > maxpower : raise ValueError("Power is too large") elif power == 0 : - return np.array([1], dtype=cs.dtype) + return np.array([1], dtype=c.dtype) elif power == 1 : - return cs + return c else : # This can be made more efficient by using powers of two # in the usual way. - prd = cs + prd = c for i in range(2, power + 1) : - prd = lagmul(prd, cs) + prd = lagmul(prd, c) return prd @@ -774,8 +771,8 @@ def lagint(c, m=1, k=[], lbnd=0, scl=1, axis=0): Note that the result of each integration is *multiplied* by `scl`. Why is this important to note? Say one is making a linear change of variable :math:`u = ax + b` in an integral relative to `x`. Then - :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a` - - perhaps not what one would have first thought. + .. math::`dx = du/a`, so one will need to set `scl` equal to + :math:`1/a` - perhaps not what one would have first thought. Also note that, in general, the result of integrating a C-series needs to be "re-projected" onto the C-series basis set. Thus, typically, @@ -842,48 +839,48 @@ def lagint(c, m=1, k=[], lbnd=0, scl=1, axis=0): def lagval(x, c, tensor=True): """ - Evaluate a Laguerre series. + Evaluate a Laguerre series at points x. - If `c` is of length ``n + 1``, this function returns the value: + If `c` is of length `n + 1`, this function returns the value: - ``p(x) = c[0]*L_0(x) + c[1]*L_1(x) + ... + c[n]*L_n(x)`` + .. math:: p(x) = c_0 * L_0(x) + c_1 * L_1(x) + ... + c_n * L_n(x) - If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)`` - will have the same shape as `x`. If `x` is a algebra_like object that - supports multiplication and addition with itself and the values in `c`, - then an object of the same type is returned. + The parameter `x` is converted to an array only if it is a tuple or a + list, otherwise it is treated as a scalar. In either case, either `x` + or its elements must support multiplication and addition both with + themselves and with the elements of `c`. - In the case where c is multidimensional, the shape of the result - depends on the value of `tensor`. If tensor is true the shape of the - return will be ``c.shape[1:] + x.shape``, where the shape of a scalar - is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if - `x` is broadcast compatible with that. + If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If + `c` is multidimensional, then the shape of the result depends on the + value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + + x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that + scalars have shape (,). - If there are trailing zeros in the coefficients they still take part in - the evaluation, so they should be avoided if efficiency is a concern. + Trailing zeros in the coefficients will be used in the evaluation, so + they should be avoided if efficiency is a concern. Parameters ---------- - x : array_like, algebra_like - If x is a list or tuple, it is converted to an ndarray. Otherwise - it is left unchanged and if it isn't an ndarray it is treated as a - scalar. In either case, `x` or any element of an ndarray must - support addition and multiplication with itself and the elements of - `c`. + x : array_like, compatible object + If `x` is a list or tuple, it is converted to an ndarray, otherwise + it is left unchanged and treated as a scalar. In either case, `x` + or its elements must support addition and multiplication with + with themselves and with the elements of `c`. c : array_like Array of coefficients ordered so that the coefficients for terms of - degree n are contained in ``c[n]``. If `c` is multidimesional the - remaining indices enumerate multiple Laguerre series. In the two + degree n are contained in c[n]. If `c` is multidimesional the + remaining indices enumerate multiple polynomials. In the two dimensional case the coefficients may be thought of as stored in the columns of `c`. tensor : boolean, optional - If true, the coefficient array shape is extended with ones on the - right, one for each dimension of `x`. Scalars are treated as having - dimension 0 for this action. The effect is that every column of - coefficients in `c` is evaluated for every value in `x`. If False, - the `x` are broadcast over the columns of `c` in the usual way. - This gives some flexibility in evaluations in the multidimensional - case. The default value it ``True``. + If True, the shape of the coefficient array is extended with ones + on the right, one for each dimension of `x`. Scalars have dimension 0 + for this action. The result is that every column of coefficients in + `c` is evaluated for every element of `x`. If False, `x` is broadcast + over the columns of `c` for the evaluation. This keyword is useful + when `c` is multidimensional. The default value is True. + + .. versionadded:: 1.7.0 Returns ------- @@ -938,31 +935,41 @@ def lagval(x, c, tensor=True): def lagval2d(x, y, c): """ - Evaluate 2D Laguerre series at points (x,y). + Evaluate a 2-D Laguerre series at points (x, y). This function returns the values: - ``p(x,y) = \sum_{i,j} c[i,j] * L_i(x) * L_j(y)`` + .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * L_i(x) * L_j(y) + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars and they + must have the same shape after conversion. In either case, either `x` + and `y` or their elements must support multiplication and addition both + with themselves and with the elements of `c`. + + If `c` is a 1-D array a one is implicitly appended to its shape to make + it 2-D. The shape of the result will be c.shape[2:] + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Laguerre seres is evaluated at the points - ``(x,y)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray. - Otherwise it is left unchanged and if it isn't an ndarray it is - treated as a scalar. See `lagval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points `(x, y)`, + where `x` and `y` must have the same shape. If `x` or `y` is a list + or tuple, it is first converted to an ndarray, otherwise it is left + unchanged and if it isn't an ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of - coefficients. + Array of coefficients ordered so that the coefficient of the term + of multi-degree i,j is contained in ``c[i,j]``. If `c` has + dimension greater than two the remaining indices enumerate multiple + sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the two dimensional Laguerre series at points formed - from pairs of corresponding values from `x` and `y`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points formed with + pairs of corresponding values from `x` and `y`. See Also -------- @@ -981,34 +988,44 @@ def lagval2d(x, y, c): def laggrid2d(x, y, c): """ - Evaluate 2D Laguerre series on the Cartesion product of x,y. + Evaluate a 2-D Laguerre series on the Cartesion product of x and y. This function returns the values: - ``p(a,b) = \sum_{i,j} c[i,j] * L_i(a) * L_j(b)`` + .. math:: p(a,b) = \sum_{i,j} c_{i,j} * L_i(a) * L_j(b) + + where the points `(a, b)` consist of all pairs formed by taking + `a` from `x` and `b` from `y`. The resulting points form a grid with + `x` in the first dimension and `y` in the second. + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars. In either + case, either `x` and `y` or their elements must support multiplication + and addition both with themselves and with the elements of `c`. - where the points ``(a,b)`` consist of all pairs of points formed by - taking ``a`` from `x` and ``b`` from `y`. The resulting points form a - grid with `x` in the first dimension and `y` in the second. + If `c` has fewer than two dimensions, ones are implicitly appended to + its shape to make it 2-D. The shape of the result will be c.shape[2:] + + x.shape + y.shape. + + .. versionadded:: 1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Laguerre series is evaluated at the points in - the Cartesian product of `x` and `y`. If `x` or `y` is a list or - tuple, it is first converted to an ndarray, Otherwise it is left - unchanged and if it isn't an ndarray it is treated as a scalar. See - `lagval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points in the + Cartesian product of `x` and `y`. If `x` or `y` is a list or + tuple, it is first converted to an ndarray, otherwise it is left + unchanged and, if it isn't an ndarray, it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + Array of coefficients ordered so that the coefficient of the term of + multi-degree i,j is contained in `c[i,j]`. If `c` has dimension + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the two dimensional Laguerre series at points in the + values : ndarray, compatible object + The values of the two dimensional Chebyshev series at points in the Cartesion product of `x` and `y`. See Also @@ -1023,32 +1040,43 @@ def laggrid2d(x, y, c): def lagval3d(x, y, z, c): """ - Evaluate 3D Laguerre series at points (x,y,z). + Evaluate a 3-D Laguerre series at points (x, y, z). This function returns the values: - ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * L_i(x) * L_j(y) * L_k(z)`` + .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * L_i(x) * L_j(y) * L_k(z) + + The parameters `x`, `y`, and `z` are converted to arrays only if + they are tuples or a lists, otherwise they are treated as a scalars and + they must have the same shape after conversion. In either case, either + `x`, `y`, and `z` or their elements must support multiplication and + addition both with themselves and with the elements of `c`. + + If `c` has fewer than 3 dimensions, ones are implicitly appended to its + shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Laguerre seres is evaluated at the points - ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape. If + x, y, z : array_like, compatible object + The three dimensional series is evaluated at the points + `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If any of `x`, `y`, or `z` is a list or tuple, it is first converted - to an ndarray. Otherwise it is left unchanged and if it isn't an - ndarray it is treated as a scalar. See `lagval` for explanation of - algebra_like. + to an ndarray, otherwise it is left unchanged and if it isn't an + ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + Array of coefficients ordered so that the coefficient of the term of + multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension + greater than 3 the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Laguerre series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the multidimension polynomial on points formed with + triples of corresponding values from `x`, `y`, and `z`. See Also -------- @@ -1068,37 +1096,48 @@ def lagval3d(x, y, z, c): def laggrid3d(x, y, z, c): """ - Evaluate 3D Laguerre series on the Cartesian product of x,y,z. + Evaluate a 3-D Laguerre series on the Cartesian product of x, y, and z. This function returns the values: - ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * L_i(a) * L_j(b) * L_k(c)`` + .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * L_i(a) * L_j(b) * L_k(c) + + where the points `(a, b, c)` consist of all triples formed by taking + `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form + a grid with `x` in the first dimension, `y` in the second, and `z` in + the third. + + The parameters `x`, `y`, and `z` are converted to arrays only if they + are tuples or a lists, otherwise they are treated as a scalars. In + either case, either `x`, `y`, and `z` or their elements must support + multiplication and addition both with themselves and with the elements + of `c`. + + If `c` has fewer than three dimensions, ones are implicitly appended to + its shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape + yshape + z.shape. - where the points ``(a,b,c)`` consist of all triples formed by taking - ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting - points form a grid with `x` in the first dimension, `y` in the second, - and `z` in the third. + .. versionadded:: 1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Laguerre seres is evaluated at the points - in the cartesian product of `x`, `y`, and `z` - ``(x,y,z)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray, - otherwise it is left unchanged and treated as a scalar. See - `lagval` for explanation of algebra_like. + x, y, z : array_like, compatible objects + The three dimensional series is evaluated at the points in the + Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a + list or tuple, it is first converted to an ndarray, otherwise it is + left unchanged and, if it isn't an ndarray, it is treated as a + scalar. c : array_like Array of coefficients ordered so that the coefficients for terms of degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Laguerre series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points in the Cartesion + product of `x` and `y`. See Also -------- @@ -1207,7 +1246,7 @@ def lagvander2d(x, y, deg) : def lagvander3d(x, y, z, deg) : - """Psuedo Vandermonde matrix of given degree. + """Pseudo Vandermonde matrix of given degree. Returns the pseudo Vandermonde matrix for 3D Laguerre series in `x`, `y`, or `z`. The sample point coordinates must all have the same shape @@ -1414,16 +1453,16 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None): return c -def lagcompanion(cs): - """Return the companion matrix of cs. +def lagcompanion(c): + """Return the companion matrix of c. The unscaled companion matrix of the Laguerre polynomials is already - symmetric when `cs` represents a single Laguerre polynomial, so no + symmetric when `c` represents a single Laguerre polynomial, so no further scaling is needed. Parameters ---------- - cs : array_like + c : array_like 1-d array of Laguerre series coefficients ordered from low to high degree. @@ -1434,56 +1473,60 @@ def lagcompanion(cs): """ accprod = np.multiply.accumulate - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) < 2: + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) < 2: raise ValueError('Series must have maximum degree of at least 1.') - if len(cs) == 2: - return np.array(1 + cs[0]/cs[1]) + if len(c) == 2: + return np.array(1 + c[0]/c[1]) - n = len(cs) - 1 - mat = np.zeros((n, n), dtype=cs.dtype) + n = len(c) - 1 + mat = np.zeros((n, n), dtype=c.dtype) top = mat.reshape(-1)[1::n+1] mid = mat.reshape(-1)[0::n+1] bot = mat.reshape(-1)[n::n+1] top[...] = -np.arange(1,n) mid[...] = 2.*np.arange(n) + 1. bot[...] = top - mat[:,-1] += (cs[:-1]/cs[-1])*n + mat[:,-1] += (c[:-1]/c[-1])*n return mat -def lagroots(cs): +def lagroots(c): """ Compute the roots of a Laguerre series. - Return the roots (a.k.a "zeros") of the Laguerre series represented by - `cs`, which is the sequence of coefficients from lowest order "term" - to highest, e.g., [1,2,3] is the series ``L_0 + 2*L_1 + 3*L_2``. + Return the roots (a.k.a. "zeros") of the polynomial + + .. math:: p(x) = \\sum_i c[i] * L_i(x). Parameters ---------- - cs : array_like - 1-d array of Laguerre series coefficients ordered from low to high. + c : 1-D array_like + 1-D array of coefficients. Returns ------- out : ndarray - Array of the roots. If all the roots are real, then so is the - dtype of ``out``; otherwise, ``out``'s dtype is complex. + Array of the roots of the series. If all the roots are real, + then `out` is also real, otherwise it is complex. See Also -------- - polyroots - chebroots + polyroots, legroots, chebroots, hermroots, hermeroots Notes ----- - Algorithm(s) used: + The root estimates are obtained as the eigenvalues of the companion + matrix, Roots far from the origin of the complex plane may have large + errors due to the numerical instability of the series for such + values. Roots with multiplicity greater than 1 will also show larger + errors as the value of the series near such points is relatively + insensitive to errors in the roots. Isolated roots near the origin can + be improved by a few iterations of Newton's method. - Remember: because the Laguerre series basis set is different from the - "standard" basis set, the results of this function *may* not be what - one is expecting. + The Laguerre series basis polynomials aren't powers of `x` so the + results of this function may seem unintuitive. Examples -------- @@ -1495,14 +1538,14 @@ def lagroots(cs): array([ -4.44089210e-16, 1.00000000e+00, 2.00000000e+00]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) <= 1 : - return np.array([], dtype=cs.dtype) - if len(cs) == 2 : - return np.array([1 + cs[0]/cs[1]]) - - m = lagcompanion(cs) + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) <= 1 : + return np.array([], dtype=c.dtype) + if len(c) == 2 : + return np.array([1 + c[0]/c[1]]) + + m = lagcompanion(c) r = la.eigvals(m) r.sort() return r diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py index f084617b6..319fb505b 100644 --- a/numpy/polynomial/legendre.py +++ b/numpy/polynomial/legendre.py @@ -115,7 +115,7 @@ def poly2leg(pol) : Returns ------- - cs : ndarray + c : ndarray 1-d array containing the coefficients of the equivalent Legendre series. @@ -147,7 +147,7 @@ def poly2leg(pol) : return res -def leg2poly(cs) : +def leg2poly(c) : """ Convert a Legendre series to a polynomial. @@ -158,7 +158,7 @@ def leg2poly(cs) : Parameters ---------- - cs : array_like + c : array_like 1-d array containing the Legendre series coefficients, ordered from lowest order term to highest. @@ -193,17 +193,17 @@ def leg2poly(cs) : """ from polynomial import polyadd, polysub, polymulx - [cs] = pu.as_series([cs]) - n = len(cs) + [c] = pu.as_series([c]) + n = len(c) if n < 3: - return cs + return c else: - c0 = cs[-2] - c1 = cs[-1] + c0 = c[-2] + c1 = c[-1] # i is the current degree of c1 for i in range(n - 1, 1, -1) : tmp = c0 - c0 = polysub(cs[i - 2], (c1*(i - 1))/i) + c0 = polysub(c[i - 2], (c1*(i - 1))/i) c1 = polyadd(tmp, (polymulx(c1)*(2*i - 1))/i) return polyadd(c0, polymulx(c1)) @@ -263,14 +263,24 @@ def legline(off, scl) : def legfromroots(roots) : """ - Generate a Legendre series with the given roots. + Generate a Legendre series with given roots. + + The function returns the coefficients of the polynomial + + .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), - Return the array of coefficients for the P-series whose roots (a.k.a. - "zeros") are given by *roots*. The returned array of coefficients is - ordered from lowest order "term" to highest, and zeros of multiplicity - greater than one must be included in *roots* a number of times equal - to their multiplicity (e.g., if `2` is a root of multiplicity three, - then [2,2,2] must be in *roots*). + in Legendre form, where the `r_n` are the roots specified in `roots`. + If a zero has multiplicity n, then it must appear in `roots` n times. + For instance, if 2 is a root of multiplicity three and 3 is a root of + multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The + roots can appear in any order. + + If the returned coefficients are `c`, then + + .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x) + + The coefficient of the last term is not generally 1 for monic + polynomials in Legendre form. Parameters ---------- @@ -280,28 +290,15 @@ def legfromroots(roots) : Returns ------- out : ndarray - 1-d array of the Legendre series coefficients, ordered from low to - high. If all roots are real, ``out.dtype`` is a float type; - otherwise, ``out.dtype`` is a complex type, even if all the - coefficients in the result are real (see Examples below). + 1-D array of coefficients. If all roots are real then `out` is a + real array, if some of the roots are complex, then `out` is complex + even if all the coefficients in the result are real (see Examples + below). See Also -------- - polyfromroots, chebfromroots - - Notes - ----- - What is returned are the :math:`c_i` such that: - - .. math:: - - \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i]) - - where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Legendre - (basis) polynomial over the domain `[-1,1]`. Note that, unlike - `polyfromroots`, due to the nature of the Legendre basis set, the - above identity *does not* imply :math:`c_n = 1` identically (see - Examples). + polyfromroots, chebfromroots, lagfromroots, hermfromroots, + hermefromroots. Examples -------- @@ -426,16 +423,16 @@ def legsub(c1, c2): return pu.trimseq(ret) -def legmulx(cs): +def legmulx(c): """Multiply a Legendre series by x. - Multiply the Legendre series `cs` by x, where x is the independent + Multiply the Legendre series `c` by x, where x is the independent variable. Parameters ---------- - cs : array_like + c : array_like 1-d array of Legendre series coefficients ordered from low to high. @@ -454,21 +451,21 @@ def legmulx(cs): xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) # The zero series needs special treatment - if len(cs) == 1 and cs[0] == 0: - return cs + if len(c) == 1 and c[0] == 0: + return c - prd = np.empty(len(cs) + 1, dtype=cs.dtype) - prd[0] = cs[0]*0 - prd[1] = cs[0] - for i in range(1, len(cs)): + prd = np.empty(len(c) + 1, dtype=c.dtype) + prd[0] = c[0]*0 + prd[1] = c[0] + for i in range(1, len(c)): j = i + 1 k = i - 1 s = i + j - prd[j] = (cs[i]*j)/s - prd[k] += (cs[i]*i)/s + prd[j] = (c[i]*j)/s + prd[k] += (c[i]*i)/s return prd @@ -516,26 +513,26 @@ def legmul(c1, c2): [c1, c2] = pu.as_series([c1, c2]) if len(c1) > len(c2): - cs = c2 + c = c2 xs = c1 else: - cs = c1 + c = c1 xs = c2 - if len(cs) == 1: - c0 = cs[0]*xs + if len(c) == 1: + c0 = c[0]*xs c1 = 0 - elif len(cs) == 2: - c0 = cs[0]*xs - c1 = cs[1]*xs + elif len(c) == 2: + c0 = c[0]*xs + c1 = c[1]*xs else : - nd = len(cs) - c0 = cs[-2]*xs - c1 = cs[-1]*xs - for i in range(3, len(cs) + 1) : + nd = len(c) + c0 = c[-2]*xs + c1 = c[-1]*xs + for i in range(3, len(c) + 1) : tmp = c0 nd = nd - 1 - c0 = legsub(cs[-i]*xs, (c1*(nd - 1))/nd) + c0 = legsub(c[-i]*xs, (c1*(nd - 1))/nd) c1 = legadd(tmp, (legmulx(c1)*(2*nd - 1))/nd) return legadd(c0, legmulx(c1)) @@ -608,16 +605,16 @@ def legdiv(c1, c2): return quo, pu.trimseq(rem) -def legpow(cs, pow, maxpower=16) : +def legpow(c, pow, maxpower=16) : """Raise a Legendre series to a power. - Returns the Legendre series `cs` raised to the power `pow`. The - arguement `cs` is a sequence of coefficients ordered from low to high. + Returns the Legendre series `c` raised to the power `pow`. The + arguement `c` is a sequence of coefficients ordered from low to high. i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.`` Parameters ---------- - cs : array_like + c : array_like 1d array of Legendre series coefficients ordered from low to high. pow : integer @@ -639,23 +636,23 @@ def legpow(cs, pow, maxpower=16) : -------- """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) power = int(pow) if power != pow or power < 0 : raise ValueError("Power must be a non-negative integer.") elif maxpower is not None and power > maxpower : raise ValueError("Power is too large") elif power == 0 : - return np.array([1], dtype=cs.dtype) + return np.array([1], dtype=c.dtype) elif power == 1 : - return cs + return c else : # This can be made more efficient by using powers of two # in the usual way. - prd = cs + prd = c for i in range(2, power + 1) : - prd = legmul(prd, cs) + prd = legmul(prd, c) return prd @@ -814,8 +811,8 @@ def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0): Note that the result of each integration is *multiplied* by `scl`. Why is this important to note? Say one is making a linear change of variable :math:`u = ax + b` in an integral relative to `x`. Then - :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a` - - perhaps not what one would have first thought. + .. math::`dx = du/a`, so one will need to set `scl` equal to + :math:`1/a` - perhaps not what one would have first thought. Also note that, in general, the result of integrating a C-series needs to be "re-projected" onto the C-series basis set. Thus, typically, @@ -826,7 +823,7 @@ def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0): -------- >>> from numpy.polynomial import legendre as L >>> c = (1,2,3) - >>> L.legint(cs) + >>> L.legint(c) array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) >>> L.legint(c, 3) array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, @@ -886,48 +883,49 @@ def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0): def legval(x, c, tensor=True): - """ Evaluate a Legendre series. + """ + Evaluate a Legendre series at points x. - If `c` is of length ``n + 1``, this function returns the value: + If `c` is of length `n + 1`, this function returns the value: - ``p(x) = c[0]*L_0(x) + c[1]*L_1(x) + ... + c[n]*L_n(x)`` + .. math:: p(x) = c_0 * L_0(x) + c_1 * L_1(x) + ... + c_n * L_n(x) - If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)`` - will have the same shape as `x`. If `x` is a algebra_like object that - supports multiplication and addition with itself and the values in `c`, - then an object of the same type is returned. + The parameter `x` is converted to an array only if it is a tuple or a + list, otherwise it is treated as a scalar. In either case, either `x` + or its elements must support multiplication and addition both with + themselves and with the elements of `c`. - In the case where c is multidimensional, the shape of the result - depends on the value of `tensor`. If tensor is true the shape of the - return will be ``c.shape[1:] + x.shape``, where the shape of a scalar - is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if - `x` is broadcast compatible with that. + If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If + `c` is multidimensional, then the shape of the result depends on the + value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + + x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that + scalars have shape (,). - If there are trailing zeros in the coefficients they still take part in - the evaluation, so they should be avoided if efficiency is a concern. + Trailing zeros in the coefficients will be used in the evaluation, so + they should be avoided if efficiency is a concern. Parameters ---------- - x : array_like, algebra_like - If x is a list or tuple, it is converted to an ndarray. Otherwise - it is left unchanged and if it isn't an ndarray it is treated as a - scalar. In either case, `x` or any element of an ndarray must - support addition and multiplication with itself and the elements of - `c`. + x : array_like, compatible object + If `x` is a list or tuple, it is converted to an ndarray, otherwise + it is left unchanged and treated as a scalar. In either case, `x` + or its elements must support addition and multiplication with + with themselves and with the elements of `c`. c : array_like Array of coefficients ordered so that the coefficients for terms of - degree n are contained in ``c[n]``. If `c` is multidimesional the - remaining indices enumerate multiple Legendre series. In the two + degree n are contained in c[n]. If `c` is multidimesional the + remaining indices enumerate multiple polynomials. In the two dimensional case the coefficients may be thought of as stored in the columns of `c`. tensor : boolean, optional - If true, the coefficient array shape is extended with ones on the - right, one for each dimension of `x`. Scalars are treated as having - dimension 0 for this action. The effect is that every column of - coefficients in `c` is evaluated for every value in `x`. If False, - the `x` are broadcast over the columns of `c` in the usual way. - This gives some flexibility in evaluations in the multidimensional - case. The default value it ``True``. + If True, the shape of the coefficient array is extended with ones + on the right, one for each dimension of `x`. Scalars have dimension 0 + for this action. The result is that every column of coefficients in + `c` is evaluated for every element of `x`. If False, `x` is broadcast + over the columns of `c` for the evaluation. This keyword is useful + when `c` is multidimensional. The default value is True. + + .. versionadded:: 1.7.0 Returns ------- @@ -974,29 +972,39 @@ def legval(x, c, tensor=True): def legval2d(x, y, c): """ - Evaluate 2D Legendre series at points (x,y). + Evaluate a 2-D Legendre series at points (x, y). This function returns the values: - ``p(x,y) = \sum_{i,j} c[i,j] * L_i(x) * L_j(y)`` + .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * L_i(x) * L_j(y) + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars and they + must have the same shape after conversion. In either case, either `x` + and `y` or their elements must support multiplication and addition both + with themselves and with the elements of `c`. + + If `c` is a 1-D array a one is implicitly appended to its shape to make + it 2-D. The shape of the result will be c.shape[2:] + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Legendre seres is evaluated at the points - ``(x,y)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray. - Otherwise it is left unchanged and if it isn't an ndarray it is - treated as a scalar. See `legval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points `(x, y)`, + where `x` and `y` must have the same shape. If `x` or `y` is a list + or tuple, it is first converted to an ndarray, otherwise it is left + unchanged and if it isn't an ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of - coefficients. + Array of coefficients ordered so that the coefficient of the term + of multi-degree i,j is contained in ``c[i,j]``. If `c` has + dimension greater than two the remaining indices enumerate multiple + sets of coefficients. Returns ------- - values : ndarray, algebra_like + values : ndarray, compatible object The values of the two dimensional Legendre series at points formed from pairs of corresponding values from `x` and `y`. @@ -1017,34 +1025,44 @@ def legval2d(x, y, c): def leggrid2d(x, y, c): """ - Evaluate 2D Legendre series on the Cartesion product of x,y. + Evaluate a 2-D Legendre series on the Cartesion product of x and y. This function returns the values: - ``p(a,b) = \sum_{i,j} c[i,j] * L_i(a) * L_j(b)`` + .. math:: p(a,b) = \sum_{i,j} c_{i,j} * L_i(a) * L_j(b) + + where the points `(a, b)` consist of all pairs formed by taking + `a` from `x` and `b` from `y`. The resulting points form a grid with + `x` in the first dimension and `y` in the second. + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars. In either + case, either `x` and `y` or their elements must support multiplication + and addition both with themselves and with the elements of `c`. - where the points ``(a,b)`` consist of all pairs of points formed by - taking ``a`` from `x` and ``b`` from `y`. The resulting points form a - grid with `x` in the first dimension and `y` in the second. + If `c` has fewer than two dimensions, ones are implicitly appended to + its shape to make it 2-D. The shape of the result will be c.shape[2:] + + x.shape + y.shape. + + .. versionadded:: 1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional Legendre series is evaluated at the points in - the Cartesian product of `x` and `y`. If `x` or `y` is a list or - tuple, it is first converted to an ndarray, Otherwise it is left - unchanged and if it isn't an ndarray it is treated as a scalar. See - `legval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points in the + Cartesian product of `x` and `y`. If `x` or `y` is a list or + tuple, it is first converted to an ndarray, otherwise it is left + unchanged and, if it isn't an ndarray, it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + Array of coefficients ordered so that the coefficient of the term of + multi-degree i,j is contained in `c[i,j]`. If `c` has dimension + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the two dimensional Legendre series at points in the + values : ndarray, compatible object + The values of the two dimensional Chebyshev series at points in the Cartesion product of `x` and `y`. See Also @@ -1059,32 +1077,43 @@ def leggrid2d(x, y, c): def legval3d(x, y, z, c): """ - Evaluate 3D Legendre series at points (x,y,z). + Evaluate a 3-D Legendre series at points (x, y, z). This function returns the values: - ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * L_i(x) * L_j(y) * L_k(z)`` + .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * L_i(x) * L_j(y) * L_k(z) + + The parameters `x`, `y`, and `z` are converted to arrays only if + they are tuples or a lists, otherwise they are treated as a scalars and + they must have the same shape after conversion. In either case, either + `x`, `y`, and `z` or their elements must support multiplication and + addition both with themselves and with the elements of `c`. + + If `c` has fewer than 3 dimensions, ones are implicitly appended to its + shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Legendre seres is evaluated at the points - ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape. If + x, y, z : array_like, compatible object + The three dimensional series is evaluated at the points + `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If any of `x`, `y`, or `z` is a list or tuple, it is first converted - to an ndarray. Otherwise it is left unchanged and if it isn't an - ndarray it is treated as a scalar. See `legval` for explanation of - algebra_like. + to an ndarray, otherwise it is left unchanged and if it isn't an + ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + Array of coefficients ordered so that the coefficient of the term of + multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension + greater than 3 the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Legendre series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the multidimension polynomial on points formed with + triples of corresponding values from `x`, `y`, and `z`. See Also -------- @@ -1104,37 +1133,48 @@ def legval3d(x, y, z, c): def leggrid3d(x, y, z, c): """ - Evaluate 3D Legendre series on the Cartesian product of x,y,z. + Evaluate a 3-D Legendre series on the Cartesian product of x, y, and z. This function returns the values: - ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * L_i(a) * L_j(b) * L_k(c)`` + .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * L_i(a) * L_j(b) * L_k(c) + + where the points `(a, b, c)` consist of all triples formed by taking + `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form + a grid with `x` in the first dimension, `y` in the second, and `z` in + the third. + + The parameters `x`, `y`, and `z` are converted to arrays only if they + are tuples or a lists, otherwise they are treated as a scalars. In + either case, either `x`, `y`, and `z` or their elements must support + multiplication and addition both with themselves and with the elements + of `c`. + + If `c` has fewer than three dimensions, ones are implicitly appended to + its shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape + yshape + z.shape. - where the points ``(a,b,c)`` consist of all triples formed by taking - ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting - points form a grid with `x` in the first dimension, `y` in the second, - and `z` in the third. + .. versionadded:: 1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional Legendre seres is evaluated at the points - in the cartesian product of `x`, `y`, and `z` - ``(x,y,z)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray, - otherwise it is left unchanged and treated as a scalar. See - `legval` for explanation of algebra_like. + x, y, z : array_like, compatible objects + The three dimensional series is evaluated at the points in the + Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a + list or tuple, it is first converted to an ndarray, otherwise it is + left unchanged and, if it isn't an ndarray, it is treated as a + scalar. c : array_like Array of coefficients ordered so that the coefficients for terms of degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the three dimensional Legendre series at points formed - from triples of corresponding values from `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points in the Cartesion + product of `x` and `y`. See Also -------- @@ -1236,7 +1276,7 @@ def legvander2d(x, y, deg) : def legvander3d(x, y, z, deg) : - """Psuedo Vandermonde matrix of given degree. + """Pseudo Vandermonde matrix of given degree. Returns the pseudo Vandermonde matrix for 3D Legendre series in `x`, `y`, or `z`. The sample point coordinates must all have the same shape after @@ -1437,18 +1477,18 @@ def legfit(x, y, deg, rcond=None, full=False, w=None): return c -def legcompanion(cs): - """Return the scaled companion matrix of cs. +def legcompanion(c): + """Return the scaled companion matrix of c. The basis polynomials are scaled so that the companion matrix is - symmetric when `cs` represents a single Legendre polynomial. This + symmetric when `c` represents a single Legendre polynomial. This provides better eigenvalue estimates than the unscaled case and in the single polynomial case the eigenvalues are guaranteed to be real if `numpy.linalg.eigvalsh` is used to obtain them. Parameters ---------- - cs : array_like + c : array_like 1-d array of Legendre series coefficients ordered from low to high degree. @@ -1458,61 +1498,56 @@ def legcompanion(cs): Scaled companion matrix of dimensions (deg, deg). """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) < 2: + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) < 2: raise ValueError('Series must have maximum degree of at least 1.') - if len(cs) == 2: - return np.array(-cs[0]/cs[1]) + if len(c) == 2: + return np.array(-c[0]/c[1]) - n = len(cs) - 1 - mat = np.zeros((n, n), dtype=cs.dtype) + n = len(c) - 1 + mat = np.zeros((n, n), dtype=c.dtype) scl = 1./np.sqrt(2*np.arange(n) + 1) top = mat.reshape(-1)[1::n+1] bot = mat.reshape(-1)[n::n+1] top[...] = np.arange(1, n)*scl[:n-1]*scl[1:n] bot[...] = top - mat[:,-1] -= (cs[:-1]/cs[-1])*(scl/scl[-1])*(n/(2*n - 1)) + mat[:,-1] -= (c[:-1]/c[-1])*(scl/scl[-1])*(n/(2*n - 1)) return mat -def legroots(cs): +def legroots(c): """ Compute the roots of a Legendre series. - Returns the roots (a.k.a "zeros") of the Legendre series represented by - `cs`, which is the sequence of coefficients from lowest order "term" - to highest, e.g., [1,2,3] is the series ``L_0 + 2*L_1 + 3*L_2``. + Return the roots (a.k.a. "zeros") of the polynomial + + .. math:: p(x) = \\sum_i c[i] * L_i(x). Parameters ---------- - cs : array_like - 1-d array of Legendre series coefficients ordered from low to high. - maxiter : int, optional - Maximum number of iterations of Newton to use in refining the - roots. + c : 1-D array_like + 1-D array of coefficients. Returns ------- out : ndarray - Sorted array of the roots. If all the roots are real, then so is - the dtype of ``out``; otherwise, ``out``'s dtype is complex. + Array of the roots of the series. If all the roots are real, + then `out` is also real, otherwise it is complex. See Also -------- - polyroots - chebroots + polyroots, chebroots, lagroots, hermroots, hermeroots Notes ----- The root estimates are obtained as the eigenvalues of the companion - matrix, Roots far from the real interval [-1, 1] in the complex plane - may have large errors due to the numerical instability of the Lengendre - series for such values. Roots with multiplicity greater than 1 will - also show larger errors as the value of the series near such points is - relatively insensitive to errors in the roots. Isolated roots near the - interval [-1, 1] can be improved by a few iterations of Newton's - method. + matrix, Roots far from the origin of the complex plane may have large + errors due to the numerical instability of the series for such + values. Roots with multiplicity greater than 1 will also show larger + errors as the value of the series near such points is relatively + insensitive to errors in the roots. Isolated roots near the origin can + be improved by a few iterations of Newton's method. The Legendre series basis polynomials aren't powers of ``x`` so the results of this function may seem unintuitive. @@ -1524,14 +1559,14 @@ def legroots(cs): array([-0.85099543, -0.11407192, 0.51506735]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) < 2: - return np.array([], dtype=cs.dtype) - if len(cs) == 2: - return np.array([-cs[0]/cs[1]]) - - m = legcompanion(cs) + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) < 2: + return np.array([], dtype=c.dtype) + if len(c) == 2: + return np.array([-c[0]/c[1]]) + + m = legcompanion(c) r = la.eigvals(m) r.sort() return r diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py index 1d8c1dbd3..73090244c 100644 --- a/numpy/polynomial/polynomial.py +++ b/numpy/polynomial/polynomial.py @@ -126,15 +126,24 @@ def polyline(off, scl) : def polyfromroots(roots) : """ - Generate a polynomial with the given roots. + Generate a monic polynomial with given roots. - Return the array of coefficients for the polynomial whose leading - coefficient (i.e., that of the highest order term) is `1` and whose - roots (a.k.a. "zeros") are given by *roots*. The returned array of - coefficients is ordered from lowest order term to highest, and zeros - of multiplicity greater than one must be included in *roots* a number - of times equal to their multiplicity (e.g., if `2` is a root of - multiplicity three, then [2,2,2] must be in *roots*). + Return the coefficients of the polynomial + + .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), + + where the `r_n` are the roots specified in `roots`. If a zero has + multiplicity n, then it must appear in `roots` n times. For instance, + if 2 is a root of multiplicity three and 3 is a root of multiplicity 2, + then `roots` looks something like [2, 2, 2, 3, 3]. The roots can appear + in any order. + + If the returned coefficients are `c`, then + + .. math:: p(x) = c_0 + c_1 * x + ... + x^n + + The coefficient of the last term is 1 for monic polynomials in this + form. Parameters ---------- @@ -144,24 +153,23 @@ def polyfromroots(roots) : Returns ------- out : ndarray - 1-d array of the polynomial's coefficients, ordered from low to - high. If all roots are real, ``out.dtype`` is a float type; - otherwise, ``out.dtype`` is a complex type, even if all the - coefficients in the result are real (see Examples below). + 1-D array of the polynomial's coefficients If all the roots are + real, then `out` is also real, otherwise it is complex. (see + Examples below). See Also -------- - chebfromroots + chebfromroots, legfromroots, lagfromroots, hermfromroots + hermefromroots Notes ----- - What is returned are the :math:`a_i` such that: + The coefficients are determined by multiplying together linear factors + of the form `(x - r_i)`, i.e. - .. math:: + .. math:: p(x) = (x - r_0) (x - r_1) ... (x - r_n) - \\sum_{i=0}^{n} a_ix^i = \\prod_{i=0}^{n} (x - roots[i]) - - where ``n == len(roots)``; note that this implies that `1` is always + where ``n == len(roots) - 1``; note that this implies that `1` is always returned for :math:`a_n`. Examples @@ -274,16 +282,16 @@ def polysub(c1, c2): return pu.trimseq(ret) -def polymulx(cs): +def polymulx(c): """Multiply a polynomial by x. - Multiply the polynomial `cs` by x, where x is the independent + Multiply the polynomial `c` by x, where x is the independent variable. Parameters ---------- - cs : array_like + c : array_like 1-d array of polynomial coefficients ordered from low to high. @@ -297,15 +305,15 @@ def polymulx(cs): .. versionadded:: 1.5.0 """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) # The zero series needs special treatment - if len(cs) == 1 and cs[0] == 0: - return cs + if len(c) == 1 and c[0] == 0: + return c - prd = np.empty(len(cs) + 1, dtype=cs.dtype) - prd[0] = cs[0]*0 - prd[1:] = cs + prd = np.empty(len(c) + 1, dtype=c.dtype) + prd[0] = c[0]*0 + prd[1:] = c return prd @@ -404,16 +412,16 @@ def polydiv(c1, c2): return c1[j+1:]/scl, pu.trimseq(c1[:j+1]) -def polypow(cs, pow, maxpower=None) : +def polypow(c, pow, maxpower=None) : """Raise a polynomial to a power. - Returns the polynomial `cs` raised to the power `pow`. The argument - `cs` is a sequence of coefficients ordered from low to high. i.e., + Returns the polynomial `c` raised to the power `pow`. The argument + `c` is a sequence of coefficients ordered from low to high. i.e., [1,2,3] is the series ``1 + 2*x + 3*x**2.`` Parameters ---------- - cs : array_like + c : array_like 1d array of array of series coefficients ordered from low to high degree. pow : integer @@ -435,23 +443,23 @@ def polypow(cs, pow, maxpower=None) : -------- """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) + # c is a trimmed copy + [c] = pu.as_series([c]) power = int(pow) if power != pow or power < 0 : raise ValueError("Power must be a non-negative integer.") elif maxpower is not None and power > maxpower : raise ValueError("Power is too large") elif power == 0 : - return np.array([1], dtype=cs.dtype) + return np.array([1], dtype=c.dtype) elif power == 1 : - return cs + return c else : # This can be made more efficient by using powers of two # in the usual way. - prd = cs + prd = c for i in range(2, power + 1) : - prd = np.convolve(prd, cs) + prd = np.convolve(prd, c) return prd @@ -592,11 +600,11 @@ def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0): Notes ----- - Note that the result of each integration is *multiplied* by `scl`. - Why is this important to note? Say one is making a linear change of - variable :math:`u = ax + b` in an integral relative to `x`. Then - :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a` - - perhaps not what one would have first thought. + Note that the result of each integration is *multiplied* by `scl`. Why + is this important to note? Say one is making a linear change of + variable :math:`u = ax + b` in an integral relative to `x`. Then + .. math::`dx = du/a`, so one will need to set `scl` equal to + :math:`1/a` - perhaps not what one would have first thought. Examples -------- @@ -660,53 +668,53 @@ def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0): def polyval(x, c, tensor=True): """ - Evaluate a polynomial. + Evaluate a polynomial at points x. - If `c` is of length ``n + 1``, this function returns the value: + If `c` is of length `n + 1`, this function returns the value - ``p(x) = c[0] + c[1]*x + ... + c[n]*x**(n)`` + .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n - If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)`` - will have the same shape as `x`. If `x` is a algebra_like object that - supports multiplication and addition with itself and the values in `c`, - then an object of the same type is returned. + The parameter `x` is converted to an array only if it is a tuple or a + list, otherwise it is treated as a scalar. In either case, either `x` + or its elements must support multiplication and addition both with + themselves and with the elements of `c`. - In the case where c is multidimensional, the shape of the result - depends on the value of `tensor`. If tensor is true the shape of the - return will be ``c.shape[1:] + x.shape``, where the shape of a scalar - is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if - `x` is broadcast compatible with that. + If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If + `c` is multidimensional, then the shape of the result depends on the + value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + + x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that + scalars have shape (,). - If there are trailing zeros in the coefficients they still take part in - the evaluation, so they should be avoided if efficiency is a concern. + Trailing zeros in the coefficients will be used in the evaluation, so + they should be avoided if efficiency is a concern. Parameters ---------- - x : array_like, algebra_like - If x is a list or tuple, it is converted to an ndarray. Otherwise - it is left unchanged and if it isn't an ndarray it is treated as a - scalar. In either case, `x` or any element of an ndarray must - support addition and multiplication with itself and the elements of - `c`. + x : array_like, compatible object + If `x` is a list or tuple, it is converted to an ndarray, otherwise + it is left unchanged and treated as a scalar. In either case, `x` + or its elements must support addition and multiplication with + with themselves and with the elements of `c`. c : array_like Array of coefficients ordered so that the coefficients for terms of - degree n are contained in ``c[n]``. If `c` is multidimesional the + degree n are contained in c[n]. If `c` is multidimesional the remaining indices enumerate multiple polynomials. In the two dimensional case the coefficients may be thought of as stored in the columns of `c`. tensor : boolean, optional - If true, the coefficient array shape is extended with ones on the - right, one for each dimension of `x`. Scalars are treated as having - dimension 0 for this action. The effect is that every column of - coefficients in `c` is evaluated for every value in `x`. If False, - the `x` are broadcast over the columns of `c` in the usual way. - This gives some flexibility in evaluations in the multidimensional - case. The default value it ``True``. + If True, the shape of the coefficient array is extended with ones + on the right, one for each dimension of `x`. Scalars have dimension 0 + for this action. The result is that every column of coefficients in + `c` is evaluated for every element of `x`. If False, `x` is broadcast + over the columns of `c` for the evaluation. This keyword is useful + when `c` is multidimensional. The default value is True. + + .. versionadded:: 1.7.0 Returns ------- - values : ndarray, algebra_like - The shape of the return value is described above. + values : ndarray, compatible object + The shape of the returned array is described above. See Also -------- @@ -755,29 +763,40 @@ def polyval(x, c, tensor=True): def polyval2d(x, y, c): """ - Evaluate 2D polynomials at points (x,y). + Evaluate a 2-D polynomial at points (x, y). - This function returns the values: + This function returns the value + + .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * x^i * y^j + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars and they + must have the same shape after conversion. In either case, either `x` + and `y` or their elements must support multiplication and addition both + with themselves and with the elements of `c`. - ``p(x,y) = \sum_{i,j} c[i,j] * x^i * y^j`` + If `c` has fewer than two dimensions, ones are implicitly appended to + its shape to make it 2-D. The shape of the result will be c.shape[2:] + + x.shape. + + .. versionadded:: 1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional polynomial is evaluated at the points - ``(x,y)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray, - otherwise it is left unchanged and if it isn't an ndarray it is - treated as a scalar. See `polyval` for explanation of algebra_like. + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points `(x, y)`, + where `x` and `y` must have the same shape. If `x` or `y` is a list + or tuple, it is first converted to an ndarray, otherwise it is left + unchanged and, if it isn't an ndarray, it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of - coefficients. + Array of coefficients ordered so that the coefficient of the term + of multi-degree i,j is contained in `c[i,j]`. If `c` has + dimension greater than two the remaining indices enumerate multiple + sets of coefficients. Returns ------- - values : ndarray, algebra_like + values : ndarray, compatible object The values of the two dimensional polynomial at points formed with pairs of corresponding values from `x` and `y`. @@ -798,33 +817,43 @@ def polyval2d(x, y, c): def polygrid2d(x, y, c): """ - Evaluate 2D polynomials on the Cartesion product of x,y. + Evaluate a 2-D polynomial on the Cartesion product of x and y. This function returns the values: - ``p(a,b) = \sum_{i,j} c[i,j] * a^i * b^j`` + .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * a^i * b^j - where the points ``(a,b)`` consist of all pairs of points formed by - taking ``a`` from `x` and ``b`` from `y`. The resulting points form a - grid with `x` in the first dimension and `y` in the second. + where the points `(a, b)` consist of all pairs formed by taking + `a` from `x` and `b` from `y`. The resulting points form a grid with + `x` in the first dimension and `y` in the second. + + The parameters `x` and `y` are converted to arrays only if they are + tuples or a lists, otherwise they are treated as a scalars. In either + case, either `x` and `y` or their elements must support multiplication + and addition both with themselves and with the elements of `c`. + + If `c` has fewer than two dimensions, ones are implicitly appended to + its shape to make it 2-D. The shape of the result will be c.shape[2:] + + x.shape + y.shape. + + .. versionadded:: 1.7.0 Parameters ---------- - x,y : array_like, algebra_like - The two dimensional polynomial is evaluated at the points in the + x, y : array_like, compatible objects + The two dimensional series is evaluated at the points in the Cartesian product of `x` and `y`. If `x` or `y` is a list or tuple, it is first converted to an ndarray, otherwise it is left - unchanged and if it isn't an ndarray it is treated as a scalar. See - `polyval` for explanation of algebra_like. + unchanged and, if it isn't an ndarray, it is treated as a scalar. c : array_like Array of coefficients ordered so that the coefficients for terms of degree i,j are contained in ``c[i,j]``. If `c` has dimension - greater than 2 the remaining indices enumerate multiple sets of + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like + values : ndarray, compatible object The values of the two dimensional polynomial at points in the Cartesion product of `x` and `y`. @@ -840,30 +869,41 @@ def polygrid2d(x, y, c): def polyval3d(x, y, z, c): """ - Evaluate 3D polynomials at points (x,y,z). + Evaluate a 3-D polynomial at points (x, y, z). This function returns the values: - ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * x^i * y^j * z^k`` + .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * x^i * y^j * z^k + + The parameters `x`, `y`, and `z` are converted to arrays only if + they are tuples or a lists, otherwise they are treated as a scalars and + they must have the same shape after conversion. In either case, either + `x`, `y`, and `z` or their elements must support multiplication and + addition both with themselves and with the elements of `c`. + + If `c` has fewer than 3 dimensions, ones are implicitly appended to its + shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape. + + .. versionadded::1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional polynomial is evaluated at the points - ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape. If + x, y, z : array_like, compatible object + The three dimensional series is evaluated at the points + `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If any of `x`, `y`, or `z` is a list or tuple, it is first converted to an ndarray, otherwise it is left unchanged and if it isn't an - ndarray it is treated as a scalar. See `polyval` for an - explanation of algebra_like. + ndarray it is treated as a scalar. c : array_like - Array of coefficients ordered so that the coefficients for terms of - degree i, j are contained in ``c[i,j,k]``. If `c` has dimension + Array of coefficients ordered so that the coefficient of the term of + multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension greater than 3 the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like + values : ndarray, compatible object The values of the multidimension polynomial on points formed with triples of corresponding values from `x`, `y`, and `z`. @@ -885,36 +925,48 @@ def polyval3d(x, y, z, c): def polygrid3d(x, y, z, c): """ - Evaluate 3D polynomials on the Cartesion product of x,y,z. + Evaluate a 3-D polynomial on the Cartesion product of x, y and z. This function returns the values: - ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * a^i * b^j * c^k`` + .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * a^i * b^j * c^k + + where the points `(a, b, c)` consist of all triples formed by taking + `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form + a grid with `x` in the first dimension, `y` in the second, and `z` in + the third. + + The parameters `x`, `y`, and `z` are converted to arrays only if they + are tuples or a lists, otherwise they are treated as a scalars. In + either case, either `x`, `y`, and `z` or their elements must support + multiplication and addition both with themselves and with the elements + of `c`. - where the points ``(a,b,c)`` consist of all triples formed by taking - ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting - points form a grid with `x` in the first dimension, `y` in the second, - and `z` in the third. + If `c` has fewer than three dimensions, ones are implicitly appended to + its shape to make it 3-D. The shape of the result will be c.shape[3:] + + x.shape + yshape + z.shape. + + .. versionadded:: 1.7.0 Parameters ---------- - x,y,z : array_like, algebra_like - The three dimensional polynomial is evaluated at the points - ``(x,y,z)``, where `x` and `y` must have the same shape. If `x` or - `y` is a list or tuple, it is first converted to an ndarray, - otherwise it is left unchanged and treated as a scalar. See - `polyval` for explanation of algebra_like. + x, y, z : array_like, compatible objects + The three dimensional series is evaluated at the points in the + Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a + list or tuple, it is first converted to an ndarray, otherwise it is + left unchanged and, if it isn't an ndarray, it is treated as a + scalar. c : array_like Array of coefficients ordered so that the coefficients for terms of - degree i,j are contained in ``c[i,j,k]``. If `c` has dimension - greater than 3 the remaining indices enumerate multiple sets of + degree i,j are contained in ``c[i,j]``. If `c` has dimension + greater than two the remaining indices enumerate multiple sets of coefficients. Returns ------- - values : ndarray, algebra_like - The values of the multidimensional polynomial on the cartesion - product of `x`, `y`, and `z`. + values : ndarray, compatible object + The values of the two dimensional polynomial at points in the Cartesion + product of `x` and `y`. See Also -------- @@ -1018,7 +1070,7 @@ def polyvander2d(x, y, deg) : def polyvander3d(x, y, z, deg) : - """Psuedo Vandermonde matrix of given degree. + """Pseudo Vandermonde matrix of given degree. Returns the pseudo Vandermonde matrix for 3D polynomials in `x`, `y`, or `z`. The sample point coordinates must all have the same shape after @@ -1246,13 +1298,13 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None): return c -def polycompanion(cs): - """Return the companion matrix of cs. +def polycompanion(c): + """Return the companion matrix of c. Parameters ---------- - cs : array_like + c : array_like 1-d array of series coefficients ordered from low to high degree. Returns @@ -1261,40 +1313,39 @@ def polycompanion(cs): Scaled companion matrix of dimensions (deg, deg). """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) < 2 : + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) < 2 : raise ValueError('Series must have maximum degree of at least 1.') - if len(cs) == 2: - return np.array(-cs[0]/cs[1]) + if len(c) == 2: + return np.array(-c[0]/c[1]) - n = len(cs) - 1 - mat = np.zeros((n, n), dtype=cs.dtype) + n = len(c) - 1 + mat = np.zeros((n, n), dtype=c.dtype) bot = mat.reshape(-1)[n::n+1] bot[...] = 1 - mat[:,-1] -= cs[:-1]/cs[-1] + mat[:,-1] -= c[:-1]/c[-1] return mat -def polyroots(cs): +def polyroots(c): """ Compute the roots of a polynomial. - Return the roots (a.k.a. "zeros") of the "polynomial" `cs`, the - polynomial's coefficients from lowest order term to highest - (e.g., [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``). + Return the roots (a.k.a. "zeros") of the polynomial + + .. math:: p(x) = \\sum_i c[i] * x^i. Parameters ---------- - cs : array_like of shape (M,) - 1-d array of polynomial coefficients ordered from low to high. + c : 1-D array_like + 1-D array of polynomial coefficients. Returns ------- out : ndarray - Array of the roots of the polynomial. If all the roots are real, - then so is the dtype of ``out``; otherwise, ``out``'s dtype is - complex. + Array of the roots of the polynomial. If all the roots are real, + then `out` is also real, otherwise it is complex. See Also -------- @@ -1322,14 +1373,14 @@ def polyroots(cs): array([ 0.00000000e+00+0.j, 0.00000000e+00+1.j, 2.77555756e-17-1.j]) """ - # cs is a trimmed copy - [cs] = pu.as_series([cs]) - if len(cs) < 2: - return np.array([], dtype=cs.dtype) - if len(cs) == 2: - return np.array([-cs[0]/cs[1]]) - - m = polycompanion(cs) + # c is a trimmed copy + [c] = pu.as_series([c]) + if len(c) < 2: + return np.array([], dtype=c.dtype) + if len(c) == 2: + return np.array([-c[0]/c[1]]) + + m = polycompanion(c) r = la.eigvals(m) r.sort() return r |