diff options
Diffstat (limited to 'numpy/polynomial/legendre.py')
| -rw-r--r-- | numpy/polynomial/legendre.py | 93 |
1 files changed, 64 insertions, 29 deletions
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 # |