Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/numpy/polynomial/chebyshev.py : 18%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2====================================================
3Chebyshev Series (:mod:`numpy.polynomial.chebyshev`)
4====================================================
6This module provides a number of objects (mostly functions) useful for
7dealing with Chebyshev series, including a `Chebyshev` class that
8encapsulates the usual arithmetic operations. (General information
9on how this module represents and works with such polynomials is in the
10docstring for its "parent" sub-package, `numpy.polynomial`).
12Classes
13-------
15.. autosummary::
16 :toctree: generated/
18 Chebyshev
21Constants
22---------
24.. autosummary::
25 :toctree: generated/
27 chebdomain
28 chebzero
29 chebone
30 chebx
32Arithmetic
33----------
35.. autosummary::
36 :toctree: generated/
38 chebadd
39 chebsub
40 chebmulx
41 chebmul
42 chebdiv
43 chebpow
44 chebval
45 chebval2d
46 chebval3d
47 chebgrid2d
48 chebgrid3d
50Calculus
51--------
53.. autosummary::
54 :toctree: generated/
56 chebder
57 chebint
59Misc Functions
60--------------
62.. autosummary::
63 :toctree: generated/
65 chebfromroots
66 chebroots
67 chebvander
68 chebvander2d
69 chebvander3d
70 chebgauss
71 chebweight
72 chebcompanion
73 chebfit
74 chebpts1
75 chebpts2
76 chebtrim
77 chebline
78 cheb2poly
79 poly2cheb
80 chebinterpolate
82See also
83--------
84`numpy.polynomial`
86Notes
87-----
88The implementations of multiplication, division, integration, and
89differentiation use the algebraic identities [1]_:
91.. math ::
92 T_n(x) = \\frac{z^n + z^{-n}}{2} \\\\
93 z\\frac{dx}{dz} = \\frac{z - z^{-1}}{2}.
95where
97.. math :: x = \\frac{z + z^{-1}}{2}.
99These identities allow a Chebyshev series to be expressed as a finite,
100symmetric Laurent series. In this module, this sort of Laurent series
101is referred to as a "z-series."
103References
104----------
105.. [1] A. T. Benjamin, et al., "Combinatorial Trigonometry with Chebyshev
106 Polynomials," *Journal of Statistical Planning and Inference 14*, 2008
107 (preprint: https://www.math.hmc.edu/~benjamin/papers/CombTrig.pdf, pg. 4)
109"""
110import numpy as np
111import numpy.linalg as la
112from numpy.core.multiarray import normalize_axis_index
114from . import polyutils as pu
115from ._polybase import ABCPolyBase
117__all__ = [
118 'chebzero', 'chebone', 'chebx', 'chebdomain', 'chebline', 'chebadd',
119 'chebsub', 'chebmulx', 'chebmul', 'chebdiv', 'chebpow', 'chebval',
120 'chebder', 'chebint', 'cheb2poly', 'poly2cheb', 'chebfromroots',
121 'chebvander', 'chebfit', 'chebtrim', 'chebroots', 'chebpts1',
122 'chebpts2', 'Chebyshev', 'chebval2d', 'chebval3d', 'chebgrid2d',
123 'chebgrid3d', 'chebvander2d', 'chebvander3d', 'chebcompanion',
124 'chebgauss', 'chebweight', 'chebinterpolate']
126chebtrim = pu.trimcoef
128#
129# A collection of functions for manipulating z-series. These are private
130# functions and do minimal error checking.
131#
133def _cseries_to_zseries(c):
134 """Covert Chebyshev series to z-series.
136 Covert a Chebyshev series to the equivalent z-series. The result is
137 never an empty array. The dtype of the return is the same as that of
138 the input. No checks are run on the arguments as this routine is for
139 internal use.
141 Parameters
142 ----------
143 c : 1-D ndarray
144 Chebyshev coefficients, ordered from low to high
146 Returns
147 -------
148 zs : 1-D ndarray
149 Odd length symmetric z-series, ordered from low to high.
151 """
152 n = c.size
153 zs = np.zeros(2*n-1, dtype=c.dtype)
154 zs[n-1:] = c/2
155 return zs + zs[::-1]
158def _zseries_to_cseries(zs):
159 """Covert z-series to a Chebyshev series.
161 Covert a z series to the equivalent Chebyshev series. The result is
162 never an empty array. The dtype of the return is the same as that of
163 the input. No checks are run on the arguments as this routine is for
164 internal use.
166 Parameters
167 ----------
168 zs : 1-D ndarray
169 Odd length symmetric z-series, ordered from low to high.
171 Returns
172 -------
173 c : 1-D ndarray
174 Chebyshev coefficients, ordered from low to high.
176 """
177 n = (zs.size + 1)//2
178 c = zs[n-1:].copy()
179 c[1:n] *= 2
180 return c
183def _zseries_mul(z1, z2):
184 """Multiply two z-series.
186 Multiply two z-series to produce a z-series.
188 Parameters
189 ----------
190 z1, z2 : 1-D ndarray
191 The arrays must be 1-D but this is not checked.
193 Returns
194 -------
195 product : 1-D ndarray
196 The product z-series.
198 Notes
199 -----
200 This is simply convolution. If symmetric/anti-symmetric z-series are
201 denoted by S/A then the following rules apply:
203 S*S, A*A -> S
204 S*A, A*S -> A
206 """
207 return np.convolve(z1, z2)
210def _zseries_div(z1, z2):
211 """Divide the first z-series by the second.
213 Divide `z1` by `z2` and return the quotient and remainder as z-series.
214 Warning: this implementation only applies when both z1 and z2 have the
215 same symmetry, which is sufficient for present purposes.
217 Parameters
218 ----------
219 z1, z2 : 1-D ndarray
220 The arrays must be 1-D and have the same symmetry, but this is not
221 checked.
223 Returns
224 -------
226 (quotient, remainder) : 1-D ndarrays
227 Quotient and remainder as z-series.
229 Notes
230 -----
231 This is not the same as polynomial division on account of the desired form
232 of the remainder. If symmetric/anti-symmetric z-series are denoted by S/A
233 then the following rules apply:
235 S/S -> S,S
236 A/A -> S,A
238 The restriction to types of the same symmetry could be fixed but seems like
239 unneeded generality. There is no natural form for the remainder in the case
240 where there is no symmetry.
242 """
243 z1 = z1.copy()
244 z2 = z2.copy()
245 lc1 = len(z1)
246 lc2 = len(z2)
247 if lc2 == 1:
248 z1 /= z2
249 return z1, z1[:1]*0
250 elif lc1 < lc2:
251 return z1[:1]*0, z1
252 else:
253 dlen = lc1 - lc2
254 scl = z2[0]
255 z2 /= scl
256 quo = np.empty(dlen + 1, dtype=z1.dtype)
257 i = 0
258 j = dlen
259 while i < j:
260 r = z1[i]
261 quo[i] = z1[i]
262 quo[dlen - i] = r
263 tmp = r*z2
264 z1[i:i+lc2] -= tmp
265 z1[j:j+lc2] -= tmp
266 i += 1
267 j -= 1
268 r = z1[i]
269 quo[i] = r
270 tmp = r*z2
271 z1[i:i+lc2] -= tmp
272 quo /= scl
273 rem = z1[i+1:i-1+lc2].copy()
274 return quo, rem
277def _zseries_der(zs):
278 """Differentiate a z-series.
280 The derivative is with respect to x, not z. This is achieved using the
281 chain rule and the value of dx/dz given in the module notes.
283 Parameters
284 ----------
285 zs : z-series
286 The z-series to differentiate.
288 Returns
289 -------
290 derivative : z-series
291 The derivative
293 Notes
294 -----
295 The zseries for x (ns) has been multiplied by two in order to avoid
296 using floats that are incompatible with Decimal and likely other
297 specialized scalar types. This scaling has been compensated by
298 multiplying the value of zs by two also so that the two cancels in the
299 division.
301 """
302 n = len(zs)//2
303 ns = np.array([-1, 0, 1], dtype=zs.dtype)
304 zs *= np.arange(-n, n+1)*2
305 d, r = _zseries_div(zs, ns)
306 return d
309def _zseries_int(zs):
310 """Integrate a z-series.
312 The integral is with respect to x, not z. This is achieved by a change
313 of variable using dx/dz given in the module notes.
315 Parameters
316 ----------
317 zs : z-series
318 The z-series to integrate
320 Returns
321 -------
322 integral : z-series
323 The indefinite integral
325 Notes
326 -----
327 The zseries for x (ns) has been multiplied by two in order to avoid
328 using floats that are incompatible with Decimal and likely other
329 specialized scalar types. This scaling has been compensated by
330 dividing the resulting zs by two.
332 """
333 n = 1 + len(zs)//2
334 ns = np.array([-1, 0, 1], dtype=zs.dtype)
335 zs = _zseries_mul(zs, ns)
336 div = np.arange(-n, n+1)*2
337 zs[:n] /= div[:n]
338 zs[n+1:] /= div[n+1:]
339 zs[n] = 0
340 return zs
342#
343# Chebyshev series functions
344#
347def poly2cheb(pol):
348 """
349 Convert a polynomial to a Chebyshev series.
351 Convert an array representing the coefficients of a polynomial (relative
352 to the "standard" basis) ordered from lowest degree to highest, to an
353 array of the coefficients of the equivalent Chebyshev series, ordered
354 from lowest to highest degree.
356 Parameters
357 ----------
358 pol : array_like
359 1-D array containing the polynomial coefficients
361 Returns
362 -------
363 c : ndarray
364 1-D array containing the coefficients of the equivalent Chebyshev
365 series.
367 See Also
368 --------
369 cheb2poly
371 Notes
372 -----
373 The easy way to do conversions between polynomial basis sets
374 is to use the convert method of a class instance.
376 Examples
377 --------
378 >>> from numpy import polynomial as P
379 >>> p = P.Polynomial(range(4))
380 >>> p
381 Polynomial([0., 1., 2., 3.], domain=[-1, 1], window=[-1, 1])
382 >>> c = p.convert(kind=P.Chebyshev)
383 >>> c
384 Chebyshev([1. , 3.25, 1. , 0.75], domain=[-1., 1.], window=[-1., 1.])
385 >>> P.chebyshev.poly2cheb(range(4))
386 array([1. , 3.25, 1. , 0.75])
388 """
389 [pol] = pu.as_series([pol])
390 deg = len(pol) - 1
391 res = 0
392 for i in range(deg, -1, -1):
393 res = chebadd(chebmulx(res), pol[i])
394 return res
397def cheb2poly(c):
398 """
399 Convert a Chebyshev series to a polynomial.
401 Convert an array representing the coefficients of a Chebyshev series,
402 ordered from lowest degree to highest, to an array of the coefficients
403 of the equivalent polynomial (relative to the "standard" basis) ordered
404 from lowest to highest degree.
406 Parameters
407 ----------
408 c : array_like
409 1-D array containing the Chebyshev series coefficients, ordered
410 from lowest order term to highest.
412 Returns
413 -------
414 pol : ndarray
415 1-D array containing the coefficients of the equivalent polynomial
416 (relative to the "standard" basis) ordered from lowest order term
417 to highest.
419 See Also
420 --------
421 poly2cheb
423 Notes
424 -----
425 The easy way to do conversions between polynomial basis sets
426 is to use the convert method of a class instance.
428 Examples
429 --------
430 >>> from numpy import polynomial as P
431 >>> c = P.Chebyshev(range(4))
432 >>> c
433 Chebyshev([0., 1., 2., 3.], domain=[-1, 1], window=[-1, 1])
434 >>> p = c.convert(kind=P.Polynomial)
435 >>> p
436 Polynomial([-2., -8., 4., 12.], domain=[-1., 1.], window=[-1., 1.])
437 >>> P.chebyshev.cheb2poly(range(4))
438 array([-2., -8., 4., 12.])
440 """
441 from .polynomial import polyadd, polysub, polymulx
443 [c] = pu.as_series([c])
444 n = len(c)
445 if n < 3:
446 return c
447 else:
448 c0 = c[-2]
449 c1 = c[-1]
450 # i is the current degree of c1
451 for i in range(n - 1, 1, -1):
452 tmp = c0
453 c0 = polysub(c[i - 2], c1)
454 c1 = polyadd(tmp, polymulx(c1)*2)
455 return polyadd(c0, polymulx(c1))
458#
459# These are constant arrays are of integer type so as to be compatible
460# with the widest range of other types, such as Decimal.
461#
463# Chebyshev default domain.
464chebdomain = np.array([-1, 1])
466# Chebyshev coefficients representing zero.
467chebzero = np.array([0])
469# Chebyshev coefficients representing one.
470chebone = np.array([1])
472# Chebyshev coefficients representing the identity x.
473chebx = np.array([0, 1])
476def chebline(off, scl):
477 """
478 Chebyshev series whose graph is a straight line.
482 Parameters
483 ----------
484 off, scl : scalars
485 The specified line is given by ``off + scl*x``.
487 Returns
488 -------
489 y : ndarray
490 This module's representation of the Chebyshev series for
491 ``off + scl*x``.
493 See Also
494 --------
495 polyline
497 Examples
498 --------
499 >>> import numpy.polynomial.chebyshev as C
500 >>> C.chebline(3,2)
501 array([3, 2])
502 >>> C.chebval(-3, C.chebline(3,2)) # should be -3
503 -3.0
505 """
506 if scl != 0:
507 return np.array([off, scl])
508 else:
509 return np.array([off])
512def chebfromroots(roots):
513 """
514 Generate a Chebyshev series with given roots.
516 The function returns the coefficients of the polynomial
518 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
520 in Chebyshev form, where the `r_n` are the roots specified in `roots`.
521 If a zero has multiplicity n, then it must appear in `roots` n times.
522 For instance, if 2 is a root of multiplicity three and 3 is a root of
523 multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The
524 roots can appear in any order.
526 If the returned coefficients are `c`, then
528 .. math:: p(x) = c_0 + c_1 * T_1(x) + ... + c_n * T_n(x)
530 The coefficient of the last term is not generally 1 for monic
531 polynomials in Chebyshev form.
533 Parameters
534 ----------
535 roots : array_like
536 Sequence containing the roots.
538 Returns
539 -------
540 out : ndarray
541 1-D array of coefficients. If all roots are real then `out` is a
542 real array, if some of the roots are complex, then `out` is complex
543 even if all the coefficients in the result are real (see Examples
544 below).
546 See Also
547 --------
548 polyfromroots, legfromroots, lagfromroots, hermfromroots, hermefromroots
550 Examples
551 --------
552 >>> import numpy.polynomial.chebyshev as C
553 >>> C.chebfromroots((-1,0,1)) # x^3 - x relative to the standard basis
554 array([ 0. , -0.25, 0. , 0.25])
555 >>> j = complex(0,1)
556 >>> C.chebfromroots((-j,j)) # x^2 + 1 relative to the standard basis
557 array([1.5+0.j, 0. +0.j, 0.5+0.j])
559 """
560 return pu._fromroots(chebline, chebmul, roots)
563def chebadd(c1, c2):
564 """
565 Add one Chebyshev series to another.
567 Returns the sum of two Chebyshev series `c1` + `c2`. The arguments
568 are sequences of coefficients ordered from lowest order term to
569 highest, i.e., [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
571 Parameters
572 ----------
573 c1, c2 : array_like
574 1-D arrays of Chebyshev series coefficients ordered from low to
575 high.
577 Returns
578 -------
579 out : ndarray
580 Array representing the Chebyshev series of their sum.
582 See Also
583 --------
584 chebsub, chebmulx, chebmul, chebdiv, chebpow
586 Notes
587 -----
588 Unlike multiplication, division, etc., the sum of two Chebyshev series
589 is a Chebyshev series (without having to "reproject" the result onto
590 the basis set) so addition, just like that of "standard" polynomials,
591 is simply "component-wise."
593 Examples
594 --------
595 >>> from numpy.polynomial import chebyshev as C
596 >>> c1 = (1,2,3)
597 >>> c2 = (3,2,1)
598 >>> C.chebadd(c1,c2)
599 array([4., 4., 4.])
601 """
602 return pu._add(c1, c2)
605def chebsub(c1, c2):
606 """
607 Subtract one Chebyshev series from another.
609 Returns the difference of two Chebyshev series `c1` - `c2`. The
610 sequences of coefficients are from lowest order term to highest, i.e.,
611 [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
613 Parameters
614 ----------
615 c1, c2 : array_like
616 1-D arrays of Chebyshev series coefficients ordered from low to
617 high.
619 Returns
620 -------
621 out : ndarray
622 Of Chebyshev series coefficients representing their difference.
624 See Also
625 --------
626 chebadd, chebmulx, chebmul, chebdiv, chebpow
628 Notes
629 -----
630 Unlike multiplication, division, etc., the difference of two Chebyshev
631 series is a Chebyshev series (without having to "reproject" the result
632 onto the basis set) so subtraction, just like that of "standard"
633 polynomials, is simply "component-wise."
635 Examples
636 --------
637 >>> from numpy.polynomial import chebyshev as C
638 >>> c1 = (1,2,3)
639 >>> c2 = (3,2,1)
640 >>> C.chebsub(c1,c2)
641 array([-2., 0., 2.])
642 >>> C.chebsub(c2,c1) # -C.chebsub(c1,c2)
643 array([ 2., 0., -2.])
645 """
646 return pu._sub(c1, c2)
649def chebmulx(c):
650 """Multiply a Chebyshev series by x.
652 Multiply the polynomial `c` by x, where x is the independent
653 variable.
656 Parameters
657 ----------
658 c : array_like
659 1-D array of Chebyshev series coefficients ordered from low to
660 high.
662 Returns
663 -------
664 out : ndarray
665 Array representing the result of the multiplication.
667 Notes
668 -----
670 .. versionadded:: 1.5.0
672 Examples
673 --------
674 >>> from numpy.polynomial import chebyshev as C
675 >>> C.chebmulx([1,2,3])
676 array([1. , 2.5, 1. , 1.5])
678 """
679 # c is a trimmed copy
680 [c] = pu.as_series([c])
681 # The zero series needs special treatment
682 if len(c) == 1 and c[0] == 0:
683 return c
685 prd = np.empty(len(c) + 1, dtype=c.dtype)
686 prd[0] = c[0]*0
687 prd[1] = c[0]
688 if len(c) > 1:
689 tmp = c[1:]/2
690 prd[2:] = tmp
691 prd[0:-2] += tmp
692 return prd
695def chebmul(c1, c2):
696 """
697 Multiply one Chebyshev series by another.
699 Returns the product of two Chebyshev series `c1` * `c2`. The arguments
700 are sequences of coefficients, from lowest order "term" to highest,
701 e.g., [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
703 Parameters
704 ----------
705 c1, c2 : array_like
706 1-D arrays of Chebyshev series coefficients ordered from low to
707 high.
709 Returns
710 -------
711 out : ndarray
712 Of Chebyshev series coefficients representing their product.
714 See Also
715 --------
716 chebadd, chebsub, chebmulx, chebdiv, chebpow
718 Notes
719 -----
720 In general, the (polynomial) product of two C-series results in terms
721 that are not in the Chebyshev polynomial basis set. Thus, to express
722 the product as a C-series, it is typically necessary to "reproject"
723 the product onto said basis set, which typically produces
724 "unintuitive live" (but correct) results; see Examples section below.
726 Examples
727 --------
728 >>> from numpy.polynomial import chebyshev as C
729 >>> c1 = (1,2,3)
730 >>> c2 = (3,2,1)
731 >>> C.chebmul(c1,c2) # multiplication requires "reprojection"
732 array([ 6.5, 12. , 12. , 4. , 1.5])
734 """
735 # c1, c2 are trimmed copies
736 [c1, c2] = pu.as_series([c1, c2])
737 z1 = _cseries_to_zseries(c1)
738 z2 = _cseries_to_zseries(c2)
739 prd = _zseries_mul(z1, z2)
740 ret = _zseries_to_cseries(prd)
741 return pu.trimseq(ret)
744def chebdiv(c1, c2):
745 """
746 Divide one Chebyshev series by another.
748 Returns the quotient-with-remainder of two Chebyshev series
749 `c1` / `c2`. The arguments are sequences of coefficients from lowest
750 order "term" to highest, e.g., [1,2,3] represents the series
751 ``T_0 + 2*T_1 + 3*T_2``.
753 Parameters
754 ----------
755 c1, c2 : array_like
756 1-D arrays of Chebyshev series coefficients ordered from low to
757 high.
759 Returns
760 -------
761 [quo, rem] : ndarrays
762 Of Chebyshev series coefficients representing the quotient and
763 remainder.
765 See Also
766 --------
767 chebadd, chebsub, chemulx, chebmul, chebpow
769 Notes
770 -----
771 In general, the (polynomial) division of one C-series by another
772 results in quotient and remainder terms that are not in the Chebyshev
773 polynomial basis set. Thus, to express these results as C-series, it
774 is typically necessary to "reproject" the results onto said basis
775 set, which typically produces "unintuitive" (but correct) results;
776 see Examples section below.
778 Examples
779 --------
780 >>> from numpy.polynomial import chebyshev as C
781 >>> c1 = (1,2,3)
782 >>> c2 = (3,2,1)
783 >>> C.chebdiv(c1,c2) # quotient "intuitive," remainder not
784 (array([3.]), array([-8., -4.]))
785 >>> c2 = (0,1,2,3)
786 >>> C.chebdiv(c2,c1) # neither "intuitive"
787 (array([0., 2.]), array([-2., -4.]))
789 """
790 # c1, c2 are trimmed copies
791 [c1, c2] = pu.as_series([c1, c2])
792 if c2[-1] == 0:
793 raise ZeroDivisionError()
795 # note: this is more efficient than `pu._div(chebmul, c1, c2)`
796 lc1 = len(c1)
797 lc2 = len(c2)
798 if lc1 < lc2:
799 return c1[:1]*0, c1
800 elif lc2 == 1:
801 return c1/c2[-1], c1[:1]*0
802 else:
803 z1 = _cseries_to_zseries(c1)
804 z2 = _cseries_to_zseries(c2)
805 quo, rem = _zseries_div(z1, z2)
806 quo = pu.trimseq(_zseries_to_cseries(quo))
807 rem = pu.trimseq(_zseries_to_cseries(rem))
808 return quo, rem
811def chebpow(c, pow, maxpower=16):
812 """Raise a Chebyshev series to a power.
814 Returns the Chebyshev series `c` raised to the power `pow`. The
815 argument `c` is a sequence of coefficients ordered from low to high.
816 i.e., [1,2,3] is the series ``T_0 + 2*T_1 + 3*T_2.``
818 Parameters
819 ----------
820 c : array_like
821 1-D array of Chebyshev series coefficients ordered from low to
822 high.
823 pow : integer
824 Power to which the series will be raised
825 maxpower : integer, optional
826 Maximum power allowed. This is mainly to limit growth of the series
827 to unmanageable size. Default is 16
829 Returns
830 -------
831 coef : ndarray
832 Chebyshev series of power.
834 See Also
835 --------
836 chebadd, chebsub, chebmulx, chebmul, chebdiv
838 Examples
839 --------
840 >>> from numpy.polynomial import chebyshev as C
841 >>> C.chebpow([1, 2, 3, 4], 2)
842 array([15.5, 22. , 16. , ..., 12.5, 12. , 8. ])
844 """
845 # note: this is more efficient than `pu._pow(chebmul, c1, c2)`, as it
846 # avoids converting between z and c series repeatedly
848 # c is a trimmed copy
849 [c] = pu.as_series([c])
850 power = int(pow)
851 if power != pow or power < 0:
852 raise ValueError("Power must be a non-negative integer.")
853 elif maxpower is not None and power > maxpower:
854 raise ValueError("Power is too large")
855 elif power == 0:
856 return np.array([1], dtype=c.dtype)
857 elif power == 1:
858 return c
859 else:
860 # This can be made more efficient by using powers of two
861 # in the usual way.
862 zs = _cseries_to_zseries(c)
863 prd = zs
864 for i in range(2, power + 1):
865 prd = np.convolve(prd, zs)
866 return _zseries_to_cseries(prd)
869def chebder(c, m=1, scl=1, axis=0):
870 """
871 Differentiate a Chebyshev series.
873 Returns the Chebyshev series coefficients `c` differentiated `m` times
874 along `axis`. At each iteration the result is multiplied by `scl` (the
875 scaling factor is for use in a linear change of variable). The argument
876 `c` is an array of coefficients from low to high degree along each
877 axis, e.g., [1,2,3] represents the series ``1*T_0 + 2*T_1 + 3*T_2``
878 while [[1,2],[1,2]] represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) +
879 2*T_0(x)*T_1(y) + 2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is
880 ``y``.
882 Parameters
883 ----------
884 c : array_like
885 Array of Chebyshev series coefficients. If c is multidimensional
886 the different axis correspond to different variables with the
887 degree in each axis given by the corresponding index.
888 m : int, optional
889 Number of derivatives taken, must be non-negative. (Default: 1)
890 scl : scalar, optional
891 Each differentiation is multiplied by `scl`. The end result is
892 multiplication by ``scl**m``. This is for use in a linear change of
893 variable. (Default: 1)
894 axis : int, optional
895 Axis over which the derivative is taken. (Default: 0).
897 .. versionadded:: 1.7.0
899 Returns
900 -------
901 der : ndarray
902 Chebyshev series of the derivative.
904 See Also
905 --------
906 chebint
908 Notes
909 -----
910 In general, the result of differentiating a C-series needs to be
911 "reprojected" onto the C-series basis set. Thus, typically, the
912 result of this function is "unintuitive," albeit correct; see Examples
913 section below.
915 Examples
916 --------
917 >>> from numpy.polynomial import chebyshev as C
918 >>> c = (1,2,3,4)
919 >>> C.chebder(c)
920 array([14., 12., 24.])
921 >>> C.chebder(c,3)
922 array([96.])
923 >>> C.chebder(c,scl=-1)
924 array([-14., -12., -24.])
925 >>> C.chebder(c,2,-1)
926 array([12., 96.])
928 """
929 c = np.array(c, ndmin=1, copy=True)
930 if c.dtype.char in '?bBhHiIlLqQpP':
931 c = c.astype(np.double)
932 cnt = pu._deprecate_as_int(m, "the order of derivation")
933 iaxis = pu._deprecate_as_int(axis, "the axis")
934 if cnt < 0:
935 raise ValueError("The order of derivation must be non-negative")
936 iaxis = normalize_axis_index(iaxis, c.ndim)
938 if cnt == 0:
939 return c
941 c = np.moveaxis(c, iaxis, 0)
942 n = len(c)
943 if cnt >= n:
944 c = c[:1]*0
945 else:
946 for i in range(cnt):
947 n = n - 1
948 c *= scl
949 der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
950 for j in range(n, 2, -1):
951 der[j - 1] = (2*j)*c[j]
952 c[j - 2] += (j*c[j])/(j - 2)
953 if n > 1:
954 der[1] = 4*c[2]
955 der[0] = c[1]
956 c = der
957 c = np.moveaxis(c, 0, iaxis)
958 return c
961def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
962 """
963 Integrate a Chebyshev series.
965 Returns the Chebyshev series coefficients `c` integrated `m` times from
966 `lbnd` along `axis`. At each iteration the resulting series is
967 **multiplied** by `scl` and an integration constant, `k`, is added.
968 The scaling factor is for use in a linear change of variable. ("Buyer
969 beware": note that, depending on what one is doing, one may want `scl`
970 to be the reciprocal of what one might expect; for more information,
971 see the Notes section below.) The argument `c` is an array of
972 coefficients from low to high degree along each axis, e.g., [1,2,3]
973 represents the series ``T_0 + 2*T_1 + 3*T_2`` while [[1,2],[1,2]]
974 represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) + 2*T_0(x)*T_1(y) +
975 2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
977 Parameters
978 ----------
979 c : array_like
980 Array of Chebyshev series coefficients. If c is multidimensional
981 the different axis correspond to different variables with the
982 degree in each axis given by the corresponding index.
983 m : int, optional
984 Order of integration, must be positive. (Default: 1)
985 k : {[], list, scalar}, optional
986 Integration constant(s). The value of the first integral at zero
987 is the first value in the list, the value of the second integral
988 at zero is the second value, etc. If ``k == []`` (the default),
989 all constants are set to zero. If ``m == 1``, a single scalar can
990 be given instead of a list.
991 lbnd : scalar, optional
992 The lower bound of the integral. (Default: 0)
993 scl : scalar, optional
994 Following each integration the result is *multiplied* by `scl`
995 before the integration constant is added. (Default: 1)
996 axis : int, optional
997 Axis over which the integral is taken. (Default: 0).
999 .. versionadded:: 1.7.0
1001 Returns
1002 -------
1003 S : ndarray
1004 C-series coefficients of the integral.
1006 Raises
1007 ------
1008 ValueError
1009 If ``m < 1``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
1010 ``np.ndim(scl) != 0``.
1012 See Also
1013 --------
1014 chebder
1016 Notes
1017 -----
1018 Note that the result of each integration is *multiplied* by `scl`.
1019 Why is this important to note? Say one is making a linear change of
1020 variable :math:`u = ax + b` in an integral relative to `x`. Then
1021 :math:`dx = du/a`, so one will need to set `scl` equal to
1022 :math:`1/a`- perhaps not what one would have first thought.
1024 Also note that, in general, the result of integrating a C-series needs
1025 to be "reprojected" onto the C-series basis set. Thus, typically,
1026 the result of this function is "unintuitive," albeit correct; see
1027 Examples section below.
1029 Examples
1030 --------
1031 >>> from numpy.polynomial import chebyshev as C
1032 >>> c = (1,2,3)
1033 >>> C.chebint(c)
1034 array([ 0.5, -0.5, 0.5, 0.5])
1035 >>> C.chebint(c,3)
1036 array([ 0.03125 , -0.1875 , 0.04166667, -0.05208333, 0.01041667, # may vary
1037 0.00625 ])
1038 >>> C.chebint(c, k=3)
1039 array([ 3.5, -0.5, 0.5, 0.5])
1040 >>> C.chebint(c,lbnd=-2)
1041 array([ 8.5, -0.5, 0.5, 0.5])
1042 >>> C.chebint(c,scl=-2)
1043 array([-1., 1., -1., -1.])
1045 """
1046 c = np.array(c, ndmin=1, copy=True)
1047 if c.dtype.char in '?bBhHiIlLqQpP':
1048 c = c.astype(np.double)
1049 if not np.iterable(k):
1050 k = [k]
1051 cnt = pu._deprecate_as_int(m, "the order of integration")
1052 iaxis = pu._deprecate_as_int(axis, "the axis")
1053 if cnt < 0:
1054 raise ValueError("The order of integration must be non-negative")
1055 if len(k) > cnt:
1056 raise ValueError("Too many integration constants")
1057 if np.ndim(lbnd) != 0:
1058 raise ValueError("lbnd must be a scalar.")
1059 if np.ndim(scl) != 0:
1060 raise ValueError("scl must be a scalar.")
1061 iaxis = normalize_axis_index(iaxis, c.ndim)
1063 if cnt == 0:
1064 return c
1066 c = np.moveaxis(c, iaxis, 0)
1067 k = list(k) + [0]*(cnt - len(k))
1068 for i in range(cnt):
1069 n = len(c)
1070 c *= scl
1071 if n == 1 and np.all(c[0] == 0):
1072 c[0] += k[i]
1073 else:
1074 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
1075 tmp[0] = c[0]*0
1076 tmp[1] = c[0]
1077 if n > 1:
1078 tmp[2] = c[1]/4
1079 for j in range(2, n):
1080 tmp[j + 1] = c[j]/(2*(j + 1))
1081 tmp[j - 1] -= c[j]/(2*(j - 1))
1082 tmp[0] += k[i] - chebval(lbnd, tmp)
1083 c = tmp
1084 c = np.moveaxis(c, 0, iaxis)
1085 return c
1088def chebval(x, c, tensor=True):
1089 """
1090 Evaluate a Chebyshev series at points x.
1092 If `c` is of length `n + 1`, this function returns the value:
1094 .. math:: p(x) = c_0 * T_0(x) + c_1 * T_1(x) + ... + c_n * T_n(x)
1096 The parameter `x` is converted to an array only if it is a tuple or a
1097 list, otherwise it is treated as a scalar. In either case, either `x`
1098 or its elements must support multiplication and addition both with
1099 themselves and with the elements of `c`.
1101 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
1102 `c` is multidimensional, then the shape of the result depends on the
1103 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
1104 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
1105 scalars have shape (,).
1107 Trailing zeros in the coefficients will be used in the evaluation, so
1108 they should be avoided if efficiency is a concern.
1110 Parameters
1111 ----------
1112 x : array_like, compatible object
1113 If `x` is a list or tuple, it is converted to an ndarray, otherwise
1114 it is left unchanged and treated as a scalar. In either case, `x`
1115 or its elements must support addition and multiplication with
1116 with themselves and with the elements of `c`.
1117 c : array_like
1118 Array of coefficients ordered so that the coefficients for terms of
1119 degree n are contained in c[n]. If `c` is multidimensional the
1120 remaining indices enumerate multiple polynomials. In the two
1121 dimensional case the coefficients may be thought of as stored in
1122 the columns of `c`.
1123 tensor : boolean, optional
1124 If True, the shape of the coefficient array is extended with ones
1125 on the right, one for each dimension of `x`. Scalars have dimension 0
1126 for this action. The result is that every column of coefficients in
1127 `c` is evaluated for every element of `x`. If False, `x` is broadcast
1128 over the columns of `c` for the evaluation. This keyword is useful
1129 when `c` is multidimensional. The default value is True.
1131 .. versionadded:: 1.7.0
1133 Returns
1134 -------
1135 values : ndarray, algebra_like
1136 The shape of the return value is described above.
1138 See Also
1139 --------
1140 chebval2d, chebgrid2d, chebval3d, chebgrid3d
1142 Notes
1143 -----
1144 The evaluation uses Clenshaw recursion, aka synthetic division.
1146 Examples
1147 --------
1149 """
1150 c = np.array(c, ndmin=1, copy=True)
1151 if c.dtype.char in '?bBhHiIlLqQpP':
1152 c = c.astype(np.double)
1153 if isinstance(x, (tuple, list)):
1154 x = np.asarray(x)
1155 if isinstance(x, np.ndarray) and tensor:
1156 c = c.reshape(c.shape + (1,)*x.ndim)
1158 if len(c) == 1:
1159 c0 = c[0]
1160 c1 = 0
1161 elif len(c) == 2:
1162 c0 = c[0]
1163 c1 = c[1]
1164 else:
1165 x2 = 2*x
1166 c0 = c[-2]
1167 c1 = c[-1]
1168 for i in range(3, len(c) + 1):
1169 tmp = c0
1170 c0 = c[-i] - c1
1171 c1 = tmp + c1*x2
1172 return c0 + c1*x
1175def chebval2d(x, y, c):
1176 """
1177 Evaluate a 2-D Chebyshev series at points (x, y).
1179 This function returns the values:
1181 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * T_i(x) * T_j(y)
1183 The parameters `x` and `y` are converted to arrays only if they are
1184 tuples or a lists, otherwise they are treated as a scalars and they
1185 must have the same shape after conversion. In either case, either `x`
1186 and `y` or their elements must support multiplication and addition both
1187 with themselves and with the elements of `c`.
1189 If `c` is a 1-D array a one is implicitly appended to its shape to make
1190 it 2-D. The shape of the result will be c.shape[2:] + x.shape.
1192 Parameters
1193 ----------
1194 x, y : array_like, compatible objects
1195 The two dimensional series is evaluated at the points `(x, y)`,
1196 where `x` and `y` must have the same shape. If `x` or `y` is a list
1197 or tuple, it is first converted to an ndarray, otherwise it is left
1198 unchanged and if it isn't an ndarray it is treated as a scalar.
1199 c : array_like
1200 Array of coefficients ordered so that the coefficient of the term
1201 of multi-degree i,j is contained in ``c[i,j]``. If `c` has
1202 dimension greater than 2 the remaining indices enumerate multiple
1203 sets of coefficients.
1205 Returns
1206 -------
1207 values : ndarray, compatible object
1208 The values of the two dimensional Chebyshev series at points formed
1209 from pairs of corresponding values from `x` and `y`.
1211 See Also
1212 --------
1213 chebval, chebgrid2d, chebval3d, chebgrid3d
1215 Notes
1216 -----
1218 .. versionadded:: 1.7.0
1220 """
1221 return pu._valnd(chebval, c, x, y)
1224def chebgrid2d(x, y, c):
1225 """
1226 Evaluate a 2-D Chebyshev series on the Cartesian product of x and y.
1228 This function returns the values:
1230 .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * T_i(a) * T_j(b),
1232 where the points `(a, b)` consist of all pairs formed by taking
1233 `a` from `x` and `b` from `y`. The resulting points form a grid with
1234 `x` in the first dimension and `y` in the second.
1236 The parameters `x` and `y` are converted to arrays only if they are
1237 tuples or a lists, otherwise they are treated as a scalars. In either
1238 case, either `x` and `y` or their elements must support multiplication
1239 and addition both with themselves and with the elements of `c`.
1241 If `c` has fewer than two dimensions, ones are implicitly appended to
1242 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
1243 x.shape + y.shape.
1245 Parameters
1246 ----------
1247 x, y : array_like, compatible objects
1248 The two dimensional series is evaluated at the points in the
1249 Cartesian product of `x` and `y`. If `x` or `y` is a list or
1250 tuple, it is first converted to an ndarray, otherwise it is left
1251 unchanged and, if it isn't an ndarray, it is treated as a scalar.
1252 c : array_like
1253 Array of coefficients ordered so that the coefficient of the term of
1254 multi-degree i,j is contained in `c[i,j]`. If `c` has dimension
1255 greater than two the remaining indices enumerate multiple sets of
1256 coefficients.
1258 Returns
1259 -------
1260 values : ndarray, compatible object
1261 The values of the two dimensional Chebyshev series at points in the
1262 Cartesian product of `x` and `y`.
1264 See Also
1265 --------
1266 chebval, chebval2d, chebval3d, chebgrid3d
1268 Notes
1269 -----
1271 .. versionadded:: 1.7.0
1273 """
1274 return pu._gridnd(chebval, c, x, y)
1277def chebval3d(x, y, z, c):
1278 """
1279 Evaluate a 3-D Chebyshev series at points (x, y, z).
1281 This function returns the values:
1283 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * T_i(x) * T_j(y) * T_k(z)
1285 The parameters `x`, `y`, and `z` are converted to arrays only if
1286 they are tuples or a lists, otherwise they are treated as a scalars and
1287 they must have the same shape after conversion. In either case, either
1288 `x`, `y`, and `z` or their elements must support multiplication and
1289 addition both with themselves and with the elements of `c`.
1291 If `c` has fewer than 3 dimensions, ones are implicitly appended to its
1292 shape to make it 3-D. The shape of the result will be c.shape[3:] +
1293 x.shape.
1295 Parameters
1296 ----------
1297 x, y, z : array_like, compatible object
1298 The three dimensional series is evaluated at the points
1299 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
1300 any of `x`, `y`, or `z` is a list or tuple, it is first converted
1301 to an ndarray, otherwise it is left unchanged and if it isn't an
1302 ndarray it is treated as a scalar.
1303 c : array_like
1304 Array of coefficients ordered so that the coefficient of the term of
1305 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
1306 greater than 3 the remaining indices enumerate multiple sets of
1307 coefficients.
1309 Returns
1310 -------
1311 values : ndarray, compatible object
1312 The values of the multidimensional polynomial on points formed with
1313 triples of corresponding values from `x`, `y`, and `z`.
1315 See Also
1316 --------
1317 chebval, chebval2d, chebgrid2d, chebgrid3d
1319 Notes
1320 -----
1322 .. versionadded:: 1.7.0
1324 """
1325 return pu._valnd(chebval, c, x, y, z)
1328def chebgrid3d(x, y, z, c):
1329 """
1330 Evaluate a 3-D Chebyshev series on the Cartesian product of x, y, and z.
1332 This function returns the values:
1334 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * T_i(a) * T_j(b) * T_k(c)
1336 where the points `(a, b, c)` consist of all triples formed by taking
1337 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
1338 a grid with `x` in the first dimension, `y` in the second, and `z` in
1339 the third.
1341 The parameters `x`, `y`, and `z` are converted to arrays only if they
1342 are tuples or a lists, otherwise they are treated as a scalars. In
1343 either case, either `x`, `y`, and `z` or their elements must support
1344 multiplication and addition both with themselves and with the elements
1345 of `c`.
1347 If `c` has fewer than three dimensions, ones are implicitly appended to
1348 its shape to make it 3-D. The shape of the result will be c.shape[3:] +
1349 x.shape + y.shape + z.shape.
1351 Parameters
1352 ----------
1353 x, y, z : array_like, compatible objects
1354 The three dimensional series is evaluated at the points in the
1355 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
1356 list or tuple, it is first converted to an ndarray, otherwise it is
1357 left unchanged and, if it isn't an ndarray, it is treated as a
1358 scalar.
1359 c : array_like
1360 Array of coefficients ordered so that the coefficients for terms of
1361 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1362 greater than two the remaining indices enumerate multiple sets of
1363 coefficients.
1365 Returns
1366 -------
1367 values : ndarray, compatible object
1368 The values of the two dimensional polynomial at points in the Cartesian
1369 product of `x` and `y`.
1371 See Also
1372 --------
1373 chebval, chebval2d, chebgrid2d, chebval3d
1375 Notes
1376 -----
1378 .. versionadded:: 1.7.0
1380 """
1381 return pu._gridnd(chebval, c, x, y, z)
1384def chebvander(x, deg):
1385 """Pseudo-Vandermonde matrix of given degree.
1387 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
1388 `x`. The pseudo-Vandermonde matrix is defined by
1390 .. math:: V[..., i] = T_i(x),
1392 where `0 <= i <= deg`. The leading indices of `V` index the elements of
1393 `x` and the last index is the degree of the Chebyshev polynomial.
1395 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
1396 matrix ``V = chebvander(x, n)``, then ``np.dot(V, c)`` and
1397 ``chebval(x, c)`` are the same up to roundoff. This equivalence is
1398 useful both for least squares fitting and for the evaluation of a large
1399 number of Chebyshev series of the same degree and sample points.
1401 Parameters
1402 ----------
1403 x : array_like
1404 Array of points. The dtype is converted to float64 or complex128
1405 depending on whether any of the elements are complex. If `x` is
1406 scalar it is converted to a 1-D array.
1407 deg : int
1408 Degree of the resulting matrix.
1410 Returns
1411 -------
1412 vander : ndarray
1413 The pseudo Vandermonde matrix. The shape of the returned matrix is
1414 ``x.shape + (deg + 1,)``, where The last index is the degree of the
1415 corresponding Chebyshev polynomial. The dtype will be the same as
1416 the converted `x`.
1418 """
1419 ideg = pu._deprecate_as_int(deg, "deg")
1420 if ideg < 0:
1421 raise ValueError("deg must be non-negative")
1423 x = np.array(x, copy=False, ndmin=1) + 0.0
1424 dims = (ideg + 1,) + x.shape
1425 dtyp = x.dtype
1426 v = np.empty(dims, dtype=dtyp)
1427 # Use forward recursion to generate the entries.
1428 v[0] = x*0 + 1
1429 if ideg > 0:
1430 x2 = 2*x
1431 v[1] = x
1432 for i in range(2, ideg + 1):
1433 v[i] = v[i-1]*x2 - v[i-2]
1434 return np.moveaxis(v, 0, -1)
1437def chebvander2d(x, y, deg):
1438 """Pseudo-Vandermonde matrix of given degrees.
1440 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1441 points `(x, y)`. The pseudo-Vandermonde matrix is defined by
1443 .. math:: V[..., (deg[1] + 1)*i + j] = T_i(x) * T_j(y),
1445 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
1446 `V` index the points `(x, y)` and the last index encodes the degrees of
1447 the Chebyshev polynomials.
1449 If ``V = chebvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
1450 correspond to the elements of a 2-D coefficient array `c` of shape
1451 (xdeg + 1, ydeg + 1) in the order
1453 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
1455 and ``np.dot(V, c.flat)`` and ``chebval2d(x, y, c)`` will be the same
1456 up to roundoff. This equivalence is useful both for least squares
1457 fitting and for the evaluation of a large number of 2-D Chebyshev
1458 series of the same degrees and sample points.
1460 Parameters
1461 ----------
1462 x, y : array_like
1463 Arrays of point coordinates, all of the same shape. The dtypes
1464 will be converted to either float64 or complex128 depending on
1465 whether any of the elements are complex. Scalars are converted to
1466 1-D arrays.
1467 deg : list of ints
1468 List of maximum degrees of the form [x_deg, y_deg].
1470 Returns
1471 -------
1472 vander2d : ndarray
1473 The shape of the returned matrix is ``x.shape + (order,)``, where
1474 :math:`order = (deg[0]+1)*(deg[1]+1)`. The dtype will be the same
1475 as the converted `x` and `y`.
1477 See Also
1478 --------
1479 chebvander, chebvander3d, chebval2d, chebval3d
1481 Notes
1482 -----
1484 .. versionadded:: 1.7.0
1486 """
1487 return pu._vander_nd_flat((chebvander, chebvander), (x, y), deg)
1490def chebvander3d(x, y, z, deg):
1491 """Pseudo-Vandermonde matrix of given degrees.
1493 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1494 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
1495 then The pseudo-Vandermonde matrix is defined by
1497 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = T_i(x)*T_j(y)*T_k(z),
1499 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
1500 indices of `V` index the points `(x, y, z)` and the last index encodes
1501 the degrees of the Chebyshev polynomials.
1503 If ``V = chebvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
1504 of `V` correspond to the elements of a 3-D coefficient array `c` of
1505 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
1507 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
1509 and ``np.dot(V, c.flat)`` and ``chebval3d(x, y, z, c)`` will be the
1510 same up to roundoff. This equivalence is useful both for least squares
1511 fitting and for the evaluation of a large number of 3-D Chebyshev
1512 series of the same degrees and sample points.
1514 Parameters
1515 ----------
1516 x, y, z : array_like
1517 Arrays of point coordinates, all of the same shape. The dtypes will
1518 be converted to either float64 or complex128 depending on whether
1519 any of the elements are complex. Scalars are converted to 1-D
1520 arrays.
1521 deg : list of ints
1522 List of maximum degrees of the form [x_deg, y_deg, z_deg].
1524 Returns
1525 -------
1526 vander3d : ndarray
1527 The shape of the returned matrix is ``x.shape + (order,)``, where
1528 :math:`order = (deg[0]+1)*(deg[1]+1)*(deg[2]+1)`. The dtype will
1529 be the same as the converted `x`, `y`, and `z`.
1531 See Also
1532 --------
1533 chebvander, chebvander3d, chebval2d, chebval3d
1535 Notes
1536 -----
1538 .. versionadded:: 1.7.0
1540 """
1541 return pu._vander_nd_flat((chebvander, chebvander, chebvander), (x, y, z), deg)
1544def chebfit(x, y, deg, rcond=None, full=False, w=None):
1545 """
1546 Least squares fit of Chebyshev series to data.
1548 Return the coefficients of a Chebyshev series of degree `deg` that is the
1549 least squares fit to the data values `y` given at points `x`. If `y` is
1550 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
1551 fits are done, one for each column of `y`, and the resulting
1552 coefficients are stored in the corresponding columns of a 2-D return.
1553 The fitted polynomial(s) are in the form
1555 .. math:: p(x) = c_0 + c_1 * T_1(x) + ... + c_n * T_n(x),
1557 where `n` is `deg`.
1559 Parameters
1560 ----------
1561 x : array_like, shape (M,)
1562 x-coordinates of the M sample points ``(x[i], y[i])``.
1563 y : array_like, shape (M,) or (M, K)
1564 y-coordinates of the sample points. Several data sets of sample
1565 points sharing the same x-coordinates can be fitted at once by
1566 passing in a 2D-array that contains one dataset per column.
1567 deg : int or 1-D array_like
1568 Degree(s) of the fitting polynomials. If `deg` is a single integer,
1569 all terms up to and including the `deg`'th term are included in the
1570 fit. For NumPy versions >= 1.11.0 a list of integers specifying the
1571 degrees of the terms to include may be used instead.
1572 rcond : float, optional
1573 Relative condition number of the fit. Singular values smaller than
1574 this relative to the largest singular value will be ignored. The
1575 default value is len(x)*eps, where eps is the relative precision of
1576 the float type, about 2e-16 in most cases.
1577 full : bool, optional
1578 Switch determining nature of return value. When it is False (the
1579 default) just the coefficients are returned, when True diagnostic
1580 information from the singular value decomposition is also returned.
1581 w : array_like, shape (`M`,), optional
1582 Weights. If not None, the contribution of each point
1583 ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the
1584 weights are chosen so that the errors of the products ``w[i]*y[i]``
1585 all have the same variance. The default value is None.
1587 .. versionadded:: 1.5.0
1589 Returns
1590 -------
1591 coef : ndarray, shape (M,) or (M, K)
1592 Chebyshev coefficients ordered from low to high. If `y` was 2-D,
1593 the coefficients for the data in column k of `y` are in column
1594 `k`.
1596 [residuals, rank, singular_values, rcond] : list
1597 These values are only returned if `full` = True
1599 resid -- sum of squared residuals of the least squares fit
1600 rank -- the numerical rank of the scaled Vandermonde matrix
1601 sv -- singular values of the scaled Vandermonde matrix
1602 rcond -- value of `rcond`.
1604 For more details, see `linalg.lstsq`.
1606 Warns
1607 -----
1608 RankWarning
1609 The rank of the coefficient matrix in the least-squares fit is
1610 deficient. The warning is only raised if `full` = False. The
1611 warnings can be turned off by
1613 >>> import warnings
1614 >>> warnings.simplefilter('ignore', np.RankWarning)
1616 See Also
1617 --------
1618 polyfit, legfit, lagfit, hermfit, hermefit
1619 chebval : Evaluates a Chebyshev series.
1620 chebvander : Vandermonde matrix of Chebyshev series.
1621 chebweight : Chebyshev weight function.
1622 linalg.lstsq : Computes a least-squares fit from the matrix.
1623 scipy.interpolate.UnivariateSpline : Computes spline fits.
1625 Notes
1626 -----
1627 The solution is the coefficients of the Chebyshev series `p` that
1628 minimizes the sum of the weighted squared errors
1630 .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
1632 where :math:`w_j` are the weights. This problem is solved by setting up
1633 as the (typically) overdetermined matrix equation
1635 .. math:: V(x) * c = w * y,
1637 where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
1638 coefficients to be solved for, `w` are the weights, and `y` are the
1639 observed values. This equation is then solved using the singular value
1640 decomposition of `V`.
1642 If some of the singular values of `V` are so small that they are
1643 neglected, then a `RankWarning` will be issued. This means that the
1644 coefficient values may be poorly determined. Using a lower order fit
1645 will usually get rid of the warning. The `rcond` parameter can also be
1646 set to a value smaller than its default, but the resulting fit may be
1647 spurious and have large contributions from roundoff error.
1649 Fits using Chebyshev series are usually better conditioned than fits
1650 using power series, but much can depend on the distribution of the
1651 sample points and the smoothness of the data. If the quality of the fit
1652 is inadequate splines may be a good alternative.
1654 References
1655 ----------
1656 .. [1] Wikipedia, "Curve fitting",
1657 https://en.wikipedia.org/wiki/Curve_fitting
1659 Examples
1660 --------
1662 """
1663 return pu._fit(chebvander, x, y, deg, rcond, full, w)
1666def chebcompanion(c):
1667 """Return the scaled companion matrix of c.
1669 The basis polynomials are scaled so that the companion matrix is
1670 symmetric when `c` is a Chebyshev basis polynomial. This provides
1671 better eigenvalue estimates than the unscaled case and for basis
1672 polynomials the eigenvalues are guaranteed to be real if
1673 `numpy.linalg.eigvalsh` is used to obtain them.
1675 Parameters
1676 ----------
1677 c : array_like
1678 1-D array of Chebyshev series coefficients ordered from low to high
1679 degree.
1681 Returns
1682 -------
1683 mat : ndarray
1684 Scaled companion matrix of dimensions (deg, deg).
1686 Notes
1687 -----
1689 .. versionadded:: 1.7.0
1691 """
1692 # c is a trimmed copy
1693 [c] = pu.as_series([c])
1694 if len(c) < 2:
1695 raise ValueError('Series must have maximum degree of at least 1.')
1696 if len(c) == 2:
1697 return np.array([[-c[0]/c[1]]])
1699 n = len(c) - 1
1700 mat = np.zeros((n, n), dtype=c.dtype)
1701 scl = np.array([1.] + [np.sqrt(.5)]*(n-1))
1702 top = mat.reshape(-1)[1::n+1]
1703 bot = mat.reshape(-1)[n::n+1]
1704 top[0] = np.sqrt(.5)
1705 top[1:] = 1/2
1706 bot[...] = top
1707 mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*.5
1708 return mat
1711def chebroots(c):
1712 """
1713 Compute the roots of a Chebyshev series.
1715 Return the roots (a.k.a. "zeros") of the polynomial
1717 .. math:: p(x) = \\sum_i c[i] * T_i(x).
1719 Parameters
1720 ----------
1721 c : 1-D array_like
1722 1-D array of coefficients.
1724 Returns
1725 -------
1726 out : ndarray
1727 Array of the roots of the series. If all the roots are real,
1728 then `out` is also real, otherwise it is complex.
1730 See Also
1731 --------
1732 polyroots, legroots, lagroots, hermroots, hermeroots
1734 Notes
1735 -----
1736 The root estimates are obtained as the eigenvalues of the companion
1737 matrix, Roots far from the origin of the complex plane may have large
1738 errors due to the numerical instability of the series for such
1739 values. Roots with multiplicity greater than 1 will also show larger
1740 errors as the value of the series near such points is relatively
1741 insensitive to errors in the roots. Isolated roots near the origin can
1742 be improved by a few iterations of Newton's method.
1744 The Chebyshev series basis polynomials aren't powers of `x` so the
1745 results of this function may seem unintuitive.
1747 Examples
1748 --------
1749 >>> import numpy.polynomial.chebyshev as cheb
1750 >>> cheb.chebroots((-1, 1,-1, 1)) # T3 - T2 + T1 - T0 has real roots
1751 array([ -5.00000000e-01, 2.60860684e-17, 1.00000000e+00]) # may vary
1753 """
1754 # c is a trimmed copy
1755 [c] = pu.as_series([c])
1756 if len(c) < 2:
1757 return np.array([], dtype=c.dtype)
1758 if len(c) == 2:
1759 return np.array([-c[0]/c[1]])
1761 # rotated companion matrix reduces error
1762 m = chebcompanion(c)[::-1,::-1]
1763 r = la.eigvals(m)
1764 r.sort()
1765 return r
1768def chebinterpolate(func, deg, args=()):
1769 """Interpolate a function at the Chebyshev points of the first kind.
1771 Returns the Chebyshev series that interpolates `func` at the Chebyshev
1772 points of the first kind in the interval [-1, 1]. The interpolating
1773 series tends to a minmax approximation to `func` with increasing `deg`
1774 if the function is continuous in the interval.
1776 .. versionadded:: 1.14.0
1778 Parameters
1779 ----------
1780 func : function
1781 The function to be approximated. It must be a function of a single
1782 variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are
1783 extra arguments passed in the `args` parameter.
1784 deg : int
1785 Degree of the interpolating polynomial
1786 args : tuple, optional
1787 Extra arguments to be used in the function call. Default is no extra
1788 arguments.
1790 Returns
1791 -------
1792 coef : ndarray, shape (deg + 1,)
1793 Chebyshev coefficients of the interpolating series ordered from low to
1794 high.
1796 Examples
1797 --------
1798 >>> import numpy.polynomial.chebyshev as C
1799 >>> C.chebfromfunction(lambda x: np.tanh(x) + 0.5, 8)
1800 array([ 5.00000000e-01, 8.11675684e-01, -9.86864911e-17,
1801 -5.42457905e-02, -2.71387850e-16, 4.51658839e-03,
1802 2.46716228e-17, -3.79694221e-04, -3.26899002e-16])
1804 Notes
1805 -----
1807 The Chebyshev polynomials used in the interpolation are orthogonal when
1808 sampled at the Chebyshev points of the first kind. If it is desired to
1809 constrain some of the coefficients they can simply be set to the desired
1810 value after the interpolation, no new interpolation or fit is needed. This
1811 is especially useful if it is known apriori that some of coefficients are
1812 zero. For instance, if the function is even then the coefficients of the
1813 terms of odd degree in the result can be set to zero.
1815 """
1816 deg = np.asarray(deg)
1818 # check arguments.
1819 if deg.ndim > 0 or deg.dtype.kind not in 'iu' or deg.size == 0:
1820 raise TypeError("deg must be an int")
1821 if deg < 0:
1822 raise ValueError("expected deg >= 0")
1824 order = deg + 1
1825 xcheb = chebpts1(order)
1826 yfunc = func(xcheb, *args)
1827 m = chebvander(xcheb, deg)
1828 c = np.dot(m.T, yfunc)
1829 c[0] /= order
1830 c[1:] /= 0.5*order
1832 return c
1835def chebgauss(deg):
1836 """
1837 Gauss-Chebyshev quadrature.
1839 Computes the sample points and weights for Gauss-Chebyshev quadrature.
1840 These sample points and weights will correctly integrate polynomials of
1841 degree :math:`2*deg - 1` or less over the interval :math:`[-1, 1]` with
1842 the weight function :math:`f(x) = 1/\\sqrt{1 - x^2}`.
1844 Parameters
1845 ----------
1846 deg : int
1847 Number of sample points and weights. It must be >= 1.
1849 Returns
1850 -------
1851 x : ndarray
1852 1-D ndarray containing the sample points.
1853 y : ndarray
1854 1-D ndarray containing the weights.
1856 Notes
1857 -----
1859 .. versionadded:: 1.7.0
1861 The results have only been tested up to degree 100, higher degrees may
1862 be problematic. For Gauss-Chebyshev there are closed form solutions for
1863 the sample points and weights. If n = `deg`, then
1865 .. math:: x_i = \\cos(\\pi (2 i - 1) / (2 n))
1867 .. math:: w_i = \\pi / n
1869 """
1870 ideg = pu._deprecate_as_int(deg, "deg")
1871 if ideg <= 0:
1872 raise ValueError("deg must be a positive integer")
1874 x = np.cos(np.pi * np.arange(1, 2*ideg, 2) / (2.0*ideg))
1875 w = np.ones(ideg)*(np.pi/ideg)
1877 return x, w
1880def chebweight(x):
1881 """
1882 The weight function of the Chebyshev polynomials.
1884 The weight function is :math:`1/\\sqrt{1 - x^2}` and the interval of
1885 integration is :math:`[-1, 1]`. The Chebyshev polynomials are
1886 orthogonal, but not normalized, with respect to this weight function.
1888 Parameters
1889 ----------
1890 x : array_like
1891 Values at which the weight function will be computed.
1893 Returns
1894 -------
1895 w : ndarray
1896 The weight function at `x`.
1898 Notes
1899 -----
1901 .. versionadded:: 1.7.0
1903 """
1904 w = 1./(np.sqrt(1. + x) * np.sqrt(1. - x))
1905 return w
1908def chebpts1(npts):
1909 """
1910 Chebyshev points of the first kind.
1912 The Chebyshev points of the first kind are the points ``cos(x)``,
1913 where ``x = [pi*(k + .5)/npts for k in range(npts)]``.
1915 Parameters
1916 ----------
1917 npts : int
1918 Number of sample points desired.
1920 Returns
1921 -------
1922 pts : ndarray
1923 The Chebyshev points of the first kind.
1925 See Also
1926 --------
1927 chebpts2
1929 Notes
1930 -----
1932 .. versionadded:: 1.5.0
1934 """
1935 _npts = int(npts)
1936 if _npts != npts:
1937 raise ValueError("npts must be integer")
1938 if _npts < 1:
1939 raise ValueError("npts must be >= 1")
1941 x = np.linspace(-np.pi, 0, _npts, endpoint=False) + np.pi/(2*_npts)
1942 return np.cos(x)
1945def chebpts2(npts):
1946 """
1947 Chebyshev points of the second kind.
1949 The Chebyshev points of the second kind are the points ``cos(x)``,
1950 where ``x = [pi*k/(npts - 1) for k in range(npts)]``.
1952 Parameters
1953 ----------
1954 npts : int
1955 Number of sample points desired.
1957 Returns
1958 -------
1959 pts : ndarray
1960 The Chebyshev points of the second kind.
1962 Notes
1963 -----
1965 .. versionadded:: 1.5.0
1967 """
1968 _npts = int(npts)
1969 if _npts != npts:
1970 raise ValueError("npts must be integer")
1971 if _npts < 2:
1972 raise ValueError("npts must be >= 2")
1974 x = np.linspace(-np.pi, 0, _npts)
1975 return np.cos(x)
1978#
1979# Chebyshev series class
1980#
1982class Chebyshev(ABCPolyBase):
1983 """A Chebyshev series class.
1985 The Chebyshev class provides the standard Python numerical methods
1986 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
1987 methods listed below.
1989 Parameters
1990 ----------
1991 coef : array_like
1992 Chebyshev coefficients in order of increasing degree, i.e.,
1993 ``(1, 2, 3)`` gives ``1*T_0(x) + 2*T_1(x) + 3*T_2(x)``.
1994 domain : (2,) array_like, optional
1995 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
1996 to the interval ``[window[0], window[1]]`` by shifting and scaling.
1997 The default value is [-1, 1].
1998 window : (2,) array_like, optional
1999 Window, see `domain` for its use. The default value is [-1, 1].
2001 .. versionadded:: 1.6.0
2003 """
2004 # Virtual Functions
2005 _add = staticmethod(chebadd)
2006 _sub = staticmethod(chebsub)
2007 _mul = staticmethod(chebmul)
2008 _div = staticmethod(chebdiv)
2009 _pow = staticmethod(chebpow)
2010 _val = staticmethod(chebval)
2011 _int = staticmethod(chebint)
2012 _der = staticmethod(chebder)
2013 _fit = staticmethod(chebfit)
2014 _line = staticmethod(chebline)
2015 _roots = staticmethod(chebroots)
2016 _fromroots = staticmethod(chebfromroots)
2018 @classmethod
2019 def interpolate(cls, func, deg, domain=None, args=()):
2020 """Interpolate a function at the Chebyshev points of the first kind.
2022 Returns the series that interpolates `func` at the Chebyshev points of
2023 the first kind scaled and shifted to the `domain`. The resulting series
2024 tends to a minmax approximation of `func` when the function is
2025 continuous in the domain.
2027 .. versionadded:: 1.14.0
2029 Parameters
2030 ----------
2031 func : function
2032 The function to be interpolated. It must be a function of a single
2033 variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are
2034 extra arguments passed in the `args` parameter.
2035 deg : int
2036 Degree of the interpolating polynomial.
2037 domain : {None, [beg, end]}, optional
2038 Domain over which `func` is interpolated. The default is None, in
2039 which case the domain is [-1, 1].
2040 args : tuple, optional
2041 Extra arguments to be used in the function call. Default is no
2042 extra arguments.
2044 Returns
2045 -------
2046 polynomial : Chebyshev instance
2047 Interpolating Chebyshev instance.
2049 Notes
2050 -----
2051 See `numpy.polynomial.chebfromfunction` for more details.
2053 """
2054 if domain is None:
2055 domain = cls.domain
2056 xfunc = lambda x: func(pu.mapdomain(x, cls.window, domain), *args)
2057 coef = chebinterpolate(xfunc, deg)
2058 return cls(coef, domain=domain)
2060 # Virtual properties
2061 nickname = 'cheb'
2062 domain = np.array(chebdomain)
2063 window = np.array(chebdomain)
2064 basis_name = 'T'