diff options
Diffstat (limited to 'numpy/polynomial/chebyshev.py')
| -rw-r--r-- | numpy/polynomial/chebyshev.py | 149 |
1 files changed, 130 insertions, 19 deletions
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py index 49d0302e0..fe2805a03 100644 --- a/numpy/polynomial/chebyshev.py +++ b/numpy/polynomial/chebyshev.py @@ -52,6 +52,7 @@ Misc Functions - `chebline` -- Chebyshev series representing given straight line. - `cheb2poly` -- convert a Chebyshev series to a polynomial. - `poly2cheb` -- convert a polynomial to a Chebyshev series. +- `chebinterpolate` -- interpolate a function at the Chebyshev points. Classes ------- @@ -87,6 +88,7 @@ References """ from __future__ import division, absolute_import, print_function +import numbers import warnings import numpy as np import numpy.linalg as la @@ -102,7 +104,7 @@ __all__ = [ 'chebvander', 'chebfit', 'chebtrim', 'chebroots', 'chebpts1', 'chebpts2', 'Chebyshev', 'chebval2d', 'chebval3d', 'chebgrid2d', 'chebgrid3d', 'chebvander2d', 'chebvander3d', 'chebcompanion', - 'chebgauss', 'chebweight'] + 'chebgauss', 'chebweight', 'chebinterpolate'] chebtrim = pu.trimcoef @@ -359,10 +361,10 @@ def poly2cheb(pol): >>> from numpy import polynomial as P >>> p = P.Polynomial(range(4)) >>> p - Polynomial([ 0., 1., 2., 3.], [-1., 1.]) + Polynomial([ 0., 1., 2., 3.], domain=[-1, 1], window=[-1, 1]) >>> c = p.convert(kind=P.Chebyshev) >>> c - Chebyshev([ 1. , 3.25, 1. , 0.75], [-1., 1.]) + Chebyshev([ 1. , 3.25, 1. , 0.75], domain=[-1, 1], window=[-1, 1]) >>> P.poly2cheb(range(4)) array([ 1. , 3.25, 1. , 0.75]) @@ -942,7 +944,7 @@ def chebder(c, m=1, scl=1, axis=0): if cnt == 0: return c - c = np.rollaxis(c, iaxis) + c = np.moveaxis(c, iaxis, 0) n = len(c) if cnt >= n: c = c[:1]*0 @@ -958,7 +960,7 @@ def chebder(c, m=1, scl=1, axis=0): der[1] = 4*c[2] der[0] = c[1] c = der - c = np.rollaxis(c, 0, iaxis + 1) + c = np.moveaxis(c, 0, iaxis) return c @@ -1022,7 +1024,7 @@ 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:`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 @@ -1067,7 +1069,7 @@ def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0): if cnt == 0: return c - c = np.rollaxis(c, iaxis) + c = np.moveaxis(c, iaxis, 0) k = list(k) + [0]*(cnt - len(k)) for i in range(cnt): n = len(c) @@ -1086,7 +1088,7 @@ def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0): tmp[j - 1] -= c[j]/(2*(j - 1)) tmp[0] += k[i] - chebval(lbnd, tmp) c = tmp - c = np.rollaxis(c, 0, iaxis + 1) + c = np.moveaxis(c, 0, iaxis) return c @@ -1220,12 +1222,12 @@ def chebval2d(x, y, c): Notes ----- - .. versionadded::1.7.0 + .. versionadded:: 1.7.0 """ try: x, y = np.array((x, y), copy=0) - except: + except Exception: raise ValueError('x, y are incompatible') c = chebval(x, c) @@ -1280,7 +1282,7 @@ def chebgrid2d(x, y, c): Notes ----- - .. versionadded::1.7.0 + .. versionadded:: 1.7.0 """ c = chebval(x, c) @@ -1333,12 +1335,12 @@ def chebval3d(x, y, z, c): Notes ----- - .. versionadded::1.7.0 + .. versionadded:: 1.7.0 """ try: x, y, z = np.array((x, y, z), copy=0) - except: + except Exception: raise ValueError('x, y, z are incompatible') c = chebval(x, c) @@ -1397,7 +1399,7 @@ def chebgrid3d(x, y, z, c): Notes ----- - .. versionadded::1.7.0 + .. versionadded:: 1.7.0 """ c = chebval(x, c) @@ -1458,7 +1460,7 @@ def chebvander(x, deg): v[1] = x for i in range(2, ideg + 1): v[i] = v[i-1]*x2 - v[i-2] - return np.rollaxis(v, 0, v.ndim) + return np.moveaxis(v, 0, -1) def chebvander2d(x, y, deg): @@ -1508,7 +1510,7 @@ def chebvander2d(x, y, deg): Notes ----- - .. versionadded::1.7.0 + .. versionadded:: 1.7.0 """ ideg = [int(d) for d in deg] @@ -1572,7 +1574,7 @@ def chebvander3d(x, y, z, deg): Notes ----- - .. versionadded::1.7.0 + .. versionadded:: 1.7.0 """ ideg = [int(d) for d in deg] @@ -1613,7 +1615,7 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None): points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column. deg : int or 1-D array_like - Degree(s) of the fitting polynomials. If `deg` is a single integer + Degree(s) of the fitting polynomials. If `deg` is a single integer, all terms up to and including the `deg`'th term are included in the fit. For NumPy versions >= 1.11.0 a list of integers specifying the degrees of the terms to include may be used instead. @@ -1808,7 +1810,7 @@ def chebcompanion(c): Notes ----- - .. versionadded::1.7.0 + .. versionadded:: 1.7.0 """ # c is a trimmed copy @@ -1886,6 +1888,73 @@ def chebroots(c): return r +def chebinterpolate(func, deg, args=()): + """Interpolate a function at the Chebyshev points of the first kind. + + Returns the Chebyshev series that interpolates `func` at the Chebyshev + points of the first kind in the interval [-1, 1]. The interpolating + series tends to a minmax approximation to `func` with increasing `deg` + if the function is continuous in the interval. + + .. versionadded:: 1.14.0 + + Parameters + ---------- + func : function + The function to be approximated. It must be a function of a single + variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are + extra arguments passed in the `args` parameter. + deg : int + Degree of the interpolating polynomial + args : tuple, optional + Extra arguments to be used in the function call. Default is no extra + arguments. + + Returns + ------- + coef : ndarray, shape (deg + 1,) + Chebyshev coefficients of the interpolating series ordered from low to + high. + + Examples + -------- + >>> import numpy.polynomial.chebyshev as C + >>> C.chebfromfunction(lambda x: np.tanh(x) + 0.5, 8) + array([ 5.00000000e-01, 8.11675684e-01, -9.86864911e-17, + -5.42457905e-02, -2.71387850e-16, 4.51658839e-03, + 2.46716228e-17, -3.79694221e-04, -3.26899002e-16]) + + Notes + ----- + + The Chebyshev polynomials used in the interpolation are orthogonal when + sampled at the Chebyshev points of the first kind. If it is desired to + constrain some of the coefficients they can simply be set to the desired + value after the interpolation, no new interpolation or fit is needed. This + is especially useful if it is known apriori that some of coefficients are + zero. For instance, if the function is even then the coefficients of the + terms of odd degree in the result can be set to zero. + + """ + deg = np.asarray(deg) + + # check arguments. + if deg.ndim > 0 or deg.dtype.kind not in 'iu' or deg.size == 0: + raise TypeError("deg must be an int") + if deg < 0: + raise ValueError("expected deg >= 0") + + order = deg + 1 + xcheb = chebpts1(order) + yfunc = func(xcheb, *args) + m = chebvander(xcheb, deg) + c = np.dot(m.T, yfunc) + c[0] /= order + c[1:] /= 0.5*order + + return c + + def chebgauss(deg): """ Gauss-Chebyshev quadrature. @@ -2069,6 +2138,48 @@ class Chebyshev(ABCPolyBase): _roots = staticmethod(chebroots) _fromroots = staticmethod(chebfromroots) + @classmethod + def interpolate(cls, func, deg, domain=None, args=()): + """Interpolate a function at the Chebyshev points of the first kind. + + Returns the series that interpolates `func` at the Chebyshev points of + the first kind scaled and shifted to the `domain`. The resulting series + tends to a minmax approximation of `func` when the function is + continuous in the domain. + + .. versionadded:: 1.14.0 + + Parameters + ---------- + func : function + The function to be interpolated. It must be a function of a single + variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are + extra arguments passed in the `args` parameter. + deg : int + Degree of the interpolating polynomial. + domain : {None, [beg, end]}, optional + Domain over which `func` is interpolated. The default is None, in + which case the domain is [-1, 1]. + args : tuple, optional + Extra arguments to be used in the function call. Default is no + extra arguments. + + Returns + ------- + polynomial : Chebyshev instance + Interpolating Chebyshev instance. + + Notes + ----- + See `numpy.polynomial.chebfromfunction` for more details. + + """ + if domain is None: + domain = cls.domain + xfunc = lambda x: func(pu.mapdomain(x, cls.window, domain), *args) + coef = chebinterpolate(xfunc, deg) + return cls(coef, domain=domain) + # Virtual properties nickname = 'cheb' domain = np.array(chebdomain) |