diff options
| -rw-r--r-- | doc/release/1.14.0-notes.rst | 9 | ||||
| -rw-r--r-- | numpy/polynomial/chebyshev.py | 115 | ||||
| -rw-r--r-- | numpy/polynomial/tests/test_chebyshev.py | 25 | ||||
| -rw-r--r-- | numpy/polynomial/tests/test_classes.py | 25 |
4 files changed, 171 insertions, 3 deletions
diff --git a/doc/release/1.14.0-notes.rst b/doc/release/1.14.0-notes.rst index 5db0df5da..5e12b9f4a 100644 --- a/doc/release/1.14.0-notes.rst +++ b/doc/release/1.14.0-notes.rst @@ -14,7 +14,8 @@ Highlights New functions ============= -* ``parametrize`` decorator added to numpy.testing +* ``parametrize``: decorator added to numpy.testing +* ``chebinterpolate``: Interpolate function at Chebyshev points. Deprecations @@ -86,6 +87,12 @@ one in pytest. It doesn't work for classes, doesn't support nesting, and does not substitute variable names. Even so, it should be adequate to rewrite the NumPy tests. +``chebinterpolate`` function added to ``numpy.polynomial.chebyshev`` +-------------------------------------------------------------------- +The new ``chebinterpolate`` function interpolates a given function at the +Chebyshev points of the first kind. A new ``Chebyshev.interpolate`` class +method adds support for interpolation over arbitrary intervals using the scaled +and shifted Chebyshev points of the first kind. Improvements ============ diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py index b983b2fec..dbf380991 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 @@ -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. @@ -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) diff --git a/numpy/polynomial/tests/test_chebyshev.py b/numpy/polynomial/tests/test_chebyshev.py index 2e4344731..1a34f42b0 100644 --- a/numpy/polynomial/tests/test_chebyshev.py +++ b/numpy/polynomial/tests/test_chebyshev.py @@ -471,6 +471,31 @@ class TestFitting(object): assert_almost_equal(coef1, coef2) +class TestInterpolate(object): + + def f(self, x): + return x * (x - 1) * (x - 2) + + def test_raises(self): + assert_raises(ValueError, cheb.chebinterpolate, self.f, -1) + assert_raises(TypeError, cheb.chebinterpolate, self.f, 10.) + + def test_dimensions(self): + for deg in range(1, 5): + assert_(cheb.chebinterpolate(self.f, deg).shape == (deg + 1,)) + + def test_approximation(self): + + def powx(x, p): + return x**p + + x = np.linspace(-1, 1, 10) + for deg in range(0, 10): + for p in range(0, deg + 1): + c = cheb.chebinterpolate(powx, deg, (p,)) + assert_almost_equal(cheb.chebval(x, c), powx(x, p), decimal=12) + + class TestCompanion(object): def test_raises(self): diff --git a/numpy/polynomial/tests/test_classes.py b/numpy/polynomial/tests/test_classes.py index 46d721df4..2ec8277ff 100644 --- a/numpy/polynomial/tests/test_classes.py +++ b/numpy/polynomial/tests/test_classes.py @@ -583,5 +583,30 @@ def check_ufunc_override(Poly): assert_raises(TypeError, np.add, x, p) +class TestInterpolate(object): + + def f(self, x): + return x * (x - 1) * (x - 2) + + def test_raises(self): + assert_raises(ValueError, Chebyshev.interpolate, self.f, -1) + assert_raises(TypeError, Chebyshev.interpolate, self.f, 10.) + + def test_dimensions(self): + for deg in range(1, 5): + assert_(Chebyshev.interpolate(self.f, deg).degree() == deg) + + def test_approximation(self): + + def powx(x, p): + return x**p + + x = np.linspace(0, 2, 10) + for deg in range(0, 10): + for t in range(0, deg + 1): + p = Chebyshev.interpolate(powx, deg, domain=[0, 2], args=(t,)) + assert_almost_equal(p(x), powx(x, t), decimal=12) + + if __name__ == "__main__": run_module_suite() |