summaryrefslogtreecommitdiff
path: root/numpy/polynomial
diff options
context:
space:
mode:
Diffstat (limited to 'numpy/polynomial')
-rw-r--r--numpy/polynomial/chebyshev.py68
-rw-r--r--numpy/polynomial/hermite.py58
-rw-r--r--numpy/polynomial/hermite_e.py64
-rw-r--r--numpy/polynomial/laguerre.py58
-rw-r--r--numpy/polynomial/legendre.py93
-rw-r--r--numpy/polynomial/polynomial.py63
-rw-r--r--numpy/polynomial/tests/test_hermite_e.py2
7 files changed, 299 insertions, 107 deletions
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index bbf437525..c80123bd3 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -1632,6 +1632,46 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
return c
+def chebcompanion(cs):
+ """Return the scaled companion matrix of cs.
+
+ The basis polynomials are scaled so that the companion matrix is
+ symmetric when `cs` 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
+ 1-d array of Legendre series coefficients ordered from low to high
+ degree.
+
+ Returns
+ -------
+ mat : ndarray
+ Scaled companion matrix of dimensions (deg, deg).
+
+ """
+ # cs is a trimmed copy
+ [cs] = pu.as_series([cs])
+ if len(cs) < 2:
+ raise ValueError('Series must have maximum degree of at least 1.')
+ if len(cs) == 2:
+ return np.array(-cs[0]/cs[1])
+
+ n = len(cs) - 1
+ mat = np.zeros((n, n), dtype=cs.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
+ return mat
+
+
def chebroots(cs):
"""
Compute the roots of a Chebyshev series.
@@ -1666,34 +1706,22 @@ def chebroots(cs):
Examples
--------
- >>> import numpy.polynomial as P
- >>> import numpy.polynomial.chebyshev as C
- >>> P.polyroots((-1,1,-1,1)) # x^3 - x^2 + x - 1 has two complex roots
- array([ -4.99600361e-16-1.j, -4.99600361e-16+1.j, 1.00000e+00+0.j])
- >>> C.chebroots((-1,1,-1,1)) # T3 - T2 + T1 - T0 has only real roots
+ >>> import numpy.polynomial.chebyshev as cheb
+ >>> cheb.chebroots((-1, 1,-1, 1)) # T3 - T2 + T1 - T0 has real roots
array([ -5.00000000e-01, 2.60860684e-17, 1.00000000e+00])
"""
# cs is a trimmed copy
[cs] = pu.as_series([cs])
- if len(cs) <= 1 :
+ if len(cs) < 2:
return np.array([], dtype=cs.dtype)
- if len(cs) == 2 :
+ if len(cs) == 2:
return np.array([-cs[0]/cs[1]])
- n = len(cs) - 1
- cs /= cs[-1]
- cmat = np.zeros((n,n), dtype=cs.dtype)
- cmat[1, 0] = 1
- for i in range(1, n):
- cmat[i - 1, i] = .5
- if i != n - 1:
- cmat[i + 1, i] = .5
- else:
- cmat[:, i] -= cs[:-1]*.5
- roots = la.eigvals(cmat)
- roots.sort()
- return roots
+ m = chebcompanion(cs)
+ r = la.eigvals(m)
+ r.sort()
+ return r
def chebpts1(npts):
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index fadd866ee..5babfb7b1 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -1411,6 +1411,47 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
return c
+def hermcompanion(cs):
+ """Return the scaled companion matrix of cs.
+
+ The basis polynomials are scaled so that the companion matrix is
+ symmetric when `cs` 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
+ 1-d array of Legendre series coefficients ordered from low to high
+ degree.
+
+ Returns
+ -------
+ mat : ndarray
+ Scaled companion matrix of dimensions (deg, deg).
+
+ """
+ accprod = np.multiply.accumulate
+ # cs is a trimmed copy
+ [cs] = pu.as_series([cs])
+ if len(cs) < 2:
+ raise ValueError('Series must have maximum degree of at least 1.')
+ if len(cs) == 2:
+ return np.array(-.5*cs[0]/cs[1])
+
+ n = len(cs) - 1
+ mat = np.zeros((n, n), dtype=cs.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
+ return mat
+
+
def hermroots(cs):
"""
Compute the roots of a Hermite series.
@@ -1460,19 +1501,10 @@ def hermroots(cs):
if len(cs) == 2 :
return np.array([-.5*cs[0]/cs[1]])
- n = len(cs) - 1
- cs /= cs[-1]
- cmat = np.zeros((n,n), dtype=cs.dtype)
- cmat[1, 0] = .5
- for i in range(1, n):
- cmat[i - 1, i] = i
- if i != n - 1:
- cmat[i + 1, i] = .5
- else:
- cmat[:, i] -= cs[:-1]*.5
- roots = la.eigvals(cmat)
- roots.sort()
- return roots
+ m = hermcompanion(cs)
+ r = la.eigvals(m)
+ r.sort()
+ return r
#
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index 97749cca4..97ef366bc 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -1407,18 +1407,59 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
return c
+def hermecompanion(cs):
+ """Return the scaled companion matrix of cs.
+
+ The basis polynomials are scaled so that the companion matrix is
+ symmetric when `cs` 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
+ 1-d array of Legendre series coefficients ordered from low to high
+ degree.
+
+ Returns
+ -------
+ mat : ndarray
+ Scaled companion matrix of dimensions (deg, deg).
+
+ """
+ accprod = np.multiply.accumulate
+ # cs is a trimmed copy
+ [cs] = pu.as_series([cs])
+ if len(cs) < 2:
+ raise ValueError('Series must have maximum degree of at least 1.')
+ if len(cs) == 2:
+ return np.array(-cs[0]/cs[1])
+
+ n = len(cs) - 1
+ mat = np.zeros((n, n), dtype=cs.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])
+ return mat
+
+
def hermeroots(cs):
"""
Compute the roots of a Hermite series.
- Return the roots (a.k.a "zeros") of the Hermite series represented by
+ 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``.
Parameters
----------
cs : array_like
- 1-d array of Hermite series coefficients ordered from low to high.
+ 1-d array of HermiteE series coefficients ordered from low to high.
Returns
-------
@@ -1454,21 +1495,12 @@ def hermeroots(cs):
if len(cs) <= 1 :
return np.array([], dtype=cs.dtype)
if len(cs) == 2 :
- return np.array([-.5*cs[0]/cs[1]])
+ return np.array([-cs[0]/cs[1]])
- n = len(cs) - 1
- cs /= cs[-1]
- cmat = np.zeros((n,n), dtype=cs.dtype)
- cmat[1, 0] = 1
- for i in range(1, n):
- cmat[i - 1, i] = i
- if i != n - 1:
- cmat[i + 1, i] = 1
- else:
- cmat[:, i] -= cs[:-1]
- roots = la.eigvals(cmat)
- roots.sort()
- return roots
+ m = hermecompanion(cs)
+ r = la.eigvals(m)
+ r.sort()
+ return r
#
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index 5f9d2d214..5ab2d23ae 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -1410,6 +1410,45 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
return c
+def lagcompanion(cs):
+ """Return the companion matrix of cs.
+
+ The unscaled companion matrix of the Laguerre polynomials is already
+ symmetric when `cs` represents a single Laguerre polynomial, so no
+ further scaling is needed.
+
+ Parameters
+ ----------
+ cs : array_like
+ 1-d array of Laguerre series coefficients ordered from low to high
+ degree.
+
+ Returns
+ -------
+ mat : ndarray
+ Companion matrix of dimensions (deg, deg).
+
+ """
+ accprod = np.multiply.accumulate
+ # cs is a trimmed copy
+ [cs] = pu.as_series([cs])
+ if len(cs) < 2:
+ raise ValueError('Series must have maximum degree of at least 1.')
+ if len(cs) == 2:
+ return np.array(1 + cs[0]/cs[1])
+
+ n = len(cs) - 1
+ mat = np.zeros((n, n), dtype=cs.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
+ return mat
+
+
def lagroots(cs):
"""
Compute the roots of a Laguerre series.
@@ -1459,21 +1498,10 @@ def lagroots(cs):
if len(cs) == 2 :
return np.array([1 + cs[0]/cs[1]])
- n = len(cs) - 1
- cs /= cs[-1]
- cmat = np.zeros((n,n), dtype=cs.dtype)
- cmat[0, 0] = 1
- cmat[1, 0] = -1
- for i in range(1, n):
- cmat[i - 1, i] = -i
- cmat[i, i] = 2*i + 1
- if i != n - 1:
- cmat[i + 1, i] = -(i + 1)
- else:
- cmat[:, i] += cs[:-1]*(i + 1)
- roots = la.eigvals(cmat)
- roots.sort()
- return roots
+ m = lagcompanion(cs)
+ r = la.eigvals(m)
+ r.sort()
+ return r
#
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 5f956b112..b9a74501c 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -1410,11 +1410,50 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
return c
+def legcompanion(cs):
+ """Return the scaled companion matrix of cs.
+
+ The basis polynomials are scaled so that the companion matrix is
+ symmetric when `cs` 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
+ 1-d array of Legendre series coefficients ordered from low to high
+ degree.
+
+ Returns
+ -------
+ mat : ndarray
+ Scaled companion matrix of dimensions (deg, deg).
+
+ """
+ # cs is a trimmed copy
+ [cs] = pu.as_series([cs])
+ if len(cs) < 2:
+ raise ValueError('Series must have maximum degree of at least 1.')
+ if len(cs) == 2:
+ return np.array(-cs[0]/cs[1])
+
+ n = len(cs) - 1
+ mat = np.zeros((n, n), dtype=cs.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))
+ return mat
+
+
def legroots(cs):
"""
Compute the roots of a Legendre series.
- Return the roots (a.k.a "zeros") of the Legendre series represented by
+ 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``.
@@ -1422,12 +1461,15 @@ def legroots(cs):
----------
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.
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.
+ Sorted array of the roots. If all the roots are real, then so is
+ the dtype of ``out``; otherwise, ``out``'s dtype is complex.
See Also
--------
@@ -1436,43 +1478,36 @@ def legroots(cs):
Notes
-----
- Algorithm(s) used:
-
- Remember: because the Legendre series basis set is different from the
- "standard" basis set, the results of this function *may* not be what
- one is expecting.
+ 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.
+
+ The Legendre series basis polynomials aren't powers of ``x`` so the
+ results of this function may seem unintuitive.
Examples
--------
- >>> import numpy.polynomial as P
- >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots
- array([-0.60582959+0.j , -0.07208521-0.63832674j,
- -0.07208521+0.63832674j])
- >>> P.legroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots
+ >>> import numpy.polynomial.legendre as leg
+ >>> leg.legroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots
array([-0.85099543, -0.11407192, 0.51506735])
"""
# cs is a trimmed copy
[cs] = pu.as_series([cs])
- if len(cs) <= 1 :
+ if len(cs) < 2:
return np.array([], dtype=cs.dtype)
- if len(cs) == 2 :
+ if len(cs) == 2:
return np.array([-cs[0]/cs[1]])
- n = len(cs) - 1
- cs /= cs[-1]
- cmat = np.zeros((n,n), dtype=cs.dtype)
- cmat[1, 0] = 1
- for i in range(1, n):
- tmp = 2*i + 1
- cmat[i - 1, i] = i/tmp
- if i != n - 1:
- cmat[i + 1, i] = (i + 1)/tmp
- else:
- cmat[:, i] -= cs[:-1]*(i + 1)/tmp
- roots = la.eigvals(cmat)
- roots.sort()
- return roots
+ m = legcompanion(cs)
+ r = la.eigvals(m)
+ r.sort()
+ return r
#
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 345d76233..1d8c1dbd3 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -1246,6 +1246,36 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
return c
+def polycompanion(cs):
+ """Return the companion matrix of cs.
+
+
+ Parameters
+ ----------
+ cs : array_like
+ 1-d array of series coefficients ordered from low to high degree.
+
+ Returns
+ -------
+ mat : ndarray
+ Scaled companion matrix of dimensions (deg, deg).
+
+ """
+ # cs is a trimmed copy
+ [cs] = pu.as_series([cs])
+ if len(cs) < 2 :
+ raise ValueError('Series must have maximum degree of at least 1.')
+ if len(cs) == 2:
+ return np.array(-cs[0]/cs[1])
+
+ n = len(cs) - 1
+ mat = np.zeros((n, n), dtype=cs.dtype)
+ bot = mat.reshape(-1)[n::n+1]
+ bot[...] = 1
+ mat[:,-1] -= cs[:-1]/cs[-1]
+ return mat
+
+
def polyroots(cs):
"""
Compute the roots of a polynomial.
@@ -1270,32 +1300,39 @@ def polyroots(cs):
--------
chebroots
+ Notes
+ -----
+ 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 power 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.
+
Examples
--------
- >>> import numpy.polynomial as P
- >>> P.polyroots(P.polyfromroots((-1,0,1)))
+ >>> import numpy.polynomial.polynomial as poly
+ >>> poly.polyroots(poly.polyfromroots((-1,0,1)))
array([-1., 0., 1.])
- >>> P.polyroots(P.polyfromroots((-1,0,1))).dtype
+ >>> poly.polyroots(poly.polyfromroots((-1,0,1))).dtype
dtype('float64')
>>> j = complex(0,1)
- >>> P.polyroots(P.polyfromroots((-j,0,j)))
+ >>> poly.polyroots(poly.polyfromroots((-j,0,j)))
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) <= 1 :
+ if len(cs) < 2:
return np.array([], dtype=cs.dtype)
- if len(cs) == 2 :
+ if len(cs) == 2:
return np.array([-cs[0]/cs[1]])
- n = len(cs) - 1
- cmat = np.zeros((n,n), dtype=cs.dtype)
- cmat.flat[n::n+1] = 1
- cmat[:,-1] -= cs[:-1]/cs[-1]
- roots = la.eigvals(cmat)
- roots.sort()
- return roots
+ m = polycompanion(cs)
+ r = la.eigvals(m)
+ r.sort()
+ return r
#
diff --git a/numpy/polynomial/tests/test_hermite_e.py b/numpy/polynomial/tests/test_hermite_e.py
index 19d35b3d4..a29dcd240 100644
--- a/numpy/polynomial/tests/test_hermite_e.py
+++ b/numpy/polynomial/tests/test_hermite_e.py
@@ -398,7 +398,7 @@ class TestMisc(TestCase) :
def test_hermeroots(self) :
assert_almost_equal(herme.hermeroots([1]), [])
- assert_almost_equal(herme.hermeroots([1, 1]), [-.5])
+ assert_almost_equal(herme.hermeroots([1, 1]), [-1])
for i in range(2,5) :
tgt = np.linspace(-1, 1, i)
res = herme.hermeroots(herme.hermefromroots(tgt))
View Lorry Controller Status