Shortcuts on this page
r m x 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 (https://web.archive.org/web/20080221202153/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 """Convert Chebyshev series to z-series.
136 Convert 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 """Convert z-series to a Chebyshev series.
161 Convert 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.
480 Parameters
481 ----------
482 off, scl : scalars
483 The specified line is given by ``off + scl*x``.
485 Returns
486 -------
487 y : ndarray
488 This module's representation of the Chebyshev series for
489 ``off + scl*x``.
491 See Also
492 --------
493 numpy.polynomial.polynomial.polyline
494 numpy.polynomial.legendre.legline
495 numpy.polynomial.laguerre.lagline
496 numpy.polynomial.hermite.hermline
497 numpy.polynomial.hermite_e.hermeline
499 Examples
500 --------
501 >>> import numpy.polynomial.chebyshev as C
502 >>> C.chebline(3,2)
503 array([3, 2])
504 >>> C.chebval(-3, C.chebline(3,2)) # should be -3
505 -3.0
507 """
508 if scl != 0:
509 return np.array([off, scl])
510 else:
511 return np.array([off])
514def chebfromroots(roots):
515 """
516 Generate a Chebyshev series with given roots.
518 The function returns the coefficients of the polynomial
520 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
522 in Chebyshev form, where the `r_n` are the roots specified in `roots`.
523 If a zero has multiplicity n, then it must appear in `roots` n times.
524 For instance, if 2 is a root of multiplicity three and 3 is a root of
525 multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The
526 roots can appear in any order.
528 If the returned coefficients are `c`, then
530 .. math:: p(x) = c_0 + c_1 * T_1(x) + ... + c_n * T_n(x)
532 The coefficient of the last term is not generally 1 for monic
533 polynomials in Chebyshev form.
535 Parameters
536 ----------
537 roots : array_like
538 Sequence containing the roots.
540 Returns
541 -------
542 out : ndarray
543 1-D array of coefficients. If all roots are real then `out` is a
544 real array, if some of the roots are complex, then `out` is complex
545 even if all the coefficients in the result are real (see Examples
546 below).
548 See Also
549 --------
550 numpy.polynomial.polynomial.polyfromroots
551 numpy.polynomial.legendre.legfromroots
552 numpy.polynomial.laguerre.lagfromroots
553 numpy.polynomial.hermite.hermfromroots
554 numpy.polynomial.hermite_e.hermefromroots
556 Examples
557 --------
558 >>> import numpy.polynomial.chebyshev as C
559 >>> C.chebfromroots((-1,0,1)) # x^3 - x relative to the standard basis
560 array([ 0. , -0.25, 0. , 0.25])
561 >>> j = complex(0,1)
562 >>> C.chebfromroots((-j,j)) # x^2 + 1 relative to the standard basis
563 array([1.5+0.j, 0. +0.j, 0.5+0.j])
565 """
566 return pu._fromroots(chebline, chebmul, roots)
569def chebadd(c1, c2):
570 """
571 Add one Chebyshev series to another.
573 Returns the sum of two Chebyshev series `c1` + `c2`. The arguments
574 are sequences of coefficients ordered from lowest order term to
575 highest, i.e., [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
577 Parameters
578 ----------
579 c1, c2 : array_like
580 1-D arrays of Chebyshev series coefficients ordered from low to
581 high.
583 Returns
584 -------
585 out : ndarray
586 Array representing the Chebyshev series of their sum.
588 See Also
589 --------
590 chebsub, chebmulx, chebmul, chebdiv, chebpow
592 Notes
593 -----
594 Unlike multiplication, division, etc., the sum of two Chebyshev series
595 is a Chebyshev series (without having to "reproject" the result onto
596 the basis set) so addition, just like that of "standard" polynomials,
597 is simply "component-wise."
599 Examples
600 --------
601 >>> from numpy.polynomial import chebyshev as C
602 >>> c1 = (1,2,3)
603 >>> c2 = (3,2,1)
604 >>> C.chebadd(c1,c2)
605 array([4., 4., 4.])
607 """
608 return pu._add(c1, c2)
611def chebsub(c1, c2):
612 """
613 Subtract one Chebyshev series from another.
615 Returns the difference of two Chebyshev series `c1` - `c2`. The
616 sequences of coefficients are from lowest order term to highest, i.e.,
617 [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
619 Parameters
620 ----------
621 c1, c2 : array_like
622 1-D arrays of Chebyshev series coefficients ordered from low to
623 high.
625 Returns
626 -------
627 out : ndarray
628 Of Chebyshev series coefficients representing their difference.
630 See Also
631 --------
632 chebadd, chebmulx, chebmul, chebdiv, chebpow
634 Notes
635 -----
636 Unlike multiplication, division, etc., the difference of two Chebyshev
637 series is a Chebyshev series (without having to "reproject" the result
638 onto the basis set) so subtraction, just like that of "standard"
639 polynomials, is simply "component-wise."
641 Examples
642 --------
643 >>> from numpy.polynomial import chebyshev as C
644 >>> c1 = (1,2,3)
645 >>> c2 = (3,2,1)
646 >>> C.chebsub(c1,c2)
647 array([-2., 0., 2.])
648 >>> C.chebsub(c2,c1) # -C.chebsub(c1,c2)
649 array([ 2., 0., -2.])
651 """
652 return pu._sub(c1, c2)
655def chebmulx(c):
656 """Multiply a Chebyshev series by x.
658 Multiply the polynomial `c` by x, where x is the independent
659 variable.
662 Parameters
663 ----------
664 c : array_like
665 1-D array of Chebyshev series coefficients ordered from low to
666 high.
668 Returns
669 -------
670 out : ndarray
671 Array representing the result of the multiplication.
673 Notes
674 -----
676 .. versionadded:: 1.5.0
678 Examples
679 --------
680 >>> from numpy.polynomial import chebyshev as C
681 >>> C.chebmulx([1,2,3])
682 array([1. , 2.5, 1. , 1.5])
684 """
685 # c is a trimmed copy
686 [c] = pu.as_series([c])
687 # The zero series needs special treatment
688 if len(c) == 1 and c[0] == 0:
689 return c
691 prd = np.empty(len(c) + 1, dtype=c.dtype)
692 prd[0] = c[0]*0
693 prd[1] = c[0]
694 if len(c) > 1:
695 tmp = c[1:]/2
696 prd[2:] = tmp
697 prd[0:-2] += tmp
698 return prd
701def chebmul(c1, c2):
702 """
703 Multiply one Chebyshev series by another.
705 Returns the product of two Chebyshev series `c1` * `c2`. The arguments
706 are sequences of coefficients, from lowest order "term" to highest,
707 e.g., [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
709 Parameters
710 ----------
711 c1, c2 : array_like
712 1-D arrays of Chebyshev series coefficients ordered from low to
713 high.
715 Returns
716 -------
717 out : ndarray
718 Of Chebyshev series coefficients representing their product.
720 See Also
721 --------
722 chebadd, chebsub, chebmulx, chebdiv, chebpow
724 Notes
725 -----
726 In general, the (polynomial) product of two C-series results in terms
727 that are not in the Chebyshev polynomial basis set. Thus, to express
728 the product as a C-series, it is typically necessary to "reproject"
729 the product onto said basis set, which typically produces
730 "unintuitive live" (but correct) results; see Examples section below.
732 Examples
733 --------
734 >>> from numpy.polynomial import chebyshev as C
735 >>> c1 = (1,2,3)
736 >>> c2 = (3,2,1)
737 >>> C.chebmul(c1,c2) # multiplication requires "reprojection"
738 array([ 6.5, 12. , 12. , 4. , 1.5])
740 """
741 # c1, c2 are trimmed copies
742 [c1, c2] = pu.as_series([c1, c2])
743 z1 = _cseries_to_zseries(c1)
744 z2 = _cseries_to_zseries(c2)
745 prd = _zseries_mul(z1, z2)
746 ret = _zseries_to_cseries(prd)
747 return pu.trimseq(ret)
750def chebdiv(c1, c2):
751 """
752 Divide one Chebyshev series by another.
754 Returns the quotient-with-remainder of two Chebyshev series
755 `c1` / `c2`. The arguments are sequences of coefficients from lowest
756 order "term" to highest, e.g., [1,2,3] represents the series
757 ``T_0 + 2*T_1 + 3*T_2``.
759 Parameters
760 ----------
761 c1, c2 : array_like
762 1-D arrays of Chebyshev series coefficients ordered from low to
763 high.
765 Returns
766 -------
767 [quo, rem] : ndarrays
768 Of Chebyshev series coefficients representing the quotient and
769 remainder.
771 See Also
772 --------
773 chebadd, chebsub, chebmulx, chebmul, chebpow
775 Notes
776 -----
777 In general, the (polynomial) division of one C-series by another
778 results in quotient and remainder terms that are not in the Chebyshev
779 polynomial basis set. Thus, to express these results as C-series, it
780 is typically necessary to "reproject" the results onto said basis
781 set, which typically produces "unintuitive" (but correct) results;
782 see Examples section below.
784 Examples
785 --------
786 >>> from numpy.polynomial import chebyshev as C
787 >>> c1 = (1,2,3)
788 >>> c2 = (3,2,1)
789 >>> C.chebdiv(c1,c2) # quotient "intuitive," remainder not
790 (array([3.]), array([-8., -4.]))
791 >>> c2 = (0,1,2,3)
792 >>> C.chebdiv(c2,c1) # neither "intuitive"
793 (array([0., 2.]), array([-2., -4.]))
795 """
796 # c1, c2 are trimmed copies
797 [c1, c2] = pu.as_series([c1, c2])
798 if c2[-1] == 0:
799 raise ZeroDivisionError()
801 # note: this is more efficient than `pu._div(chebmul, c1, c2)`
802 lc1 = len(c1)
803 lc2 = len(c2)
804 if lc1 < lc2:
805 return c1[:1]*0, c1
806 elif lc2 == 1:
807 return c1/c2[-1], c1[:1]*0
808 else:
809 z1 = _cseries_to_zseries(c1)
810 z2 = _cseries_to_zseries(c2)
811 quo, rem = _zseries_div(z1, z2)
812 quo = pu.trimseq(_zseries_to_cseries(quo))
813 rem = pu.trimseq(_zseries_to_cseries(rem))
814 return quo, rem
817def chebpow(c, pow, maxpower=16):
818 """Raise a Chebyshev series to a power.
820 Returns the Chebyshev series `c` raised to the power `pow`. The
821 argument `c` is a sequence of coefficients ordered from low to high.
822 i.e., [1,2,3] is the series ``T_0 + 2*T_1 + 3*T_2.``
824 Parameters
825 ----------
826 c : array_like
827 1-D array of Chebyshev series coefficients ordered from low to
828 high.
829 pow : integer
830 Power to which the series will be raised
831 maxpower : integer, optional
832 Maximum power allowed. This is mainly to limit growth of the series
833 to unmanageable size. Default is 16
835 Returns
836 -------
837 coef : ndarray
838 Chebyshev series of power.
840 See Also
841 --------
842 chebadd, chebsub, chebmulx, chebmul, chebdiv
844 Examples
845 --------
846 >>> from numpy.polynomial import chebyshev as C
847 >>> C.chebpow([1, 2, 3, 4], 2)
848 array([15.5, 22. , 16. , ..., 12.5, 12. , 8. ])
850 """
851 # note: this is more efficient than `pu._pow(chebmul, c1, c2)`, as it
852 # avoids converting between z and c series repeatedly
854 # c is a trimmed copy
855 [c] = pu.as_series([c])
856 power = int(pow)
857 if power != pow or power < 0:
858 raise ValueError("Power must be a non-negative integer.")
859 elif maxpower is not None and power > maxpower:
860 raise ValueError("Power is too large")
861 elif power == 0:
862 return np.array([1], dtype=c.dtype)
863 elif power == 1:
864 return c
865 else:
866 # This can be made more efficient by using powers of two
867 # in the usual way.
868 zs = _cseries_to_zseries(c)
869 prd = zs
870 for i in range(2, power + 1):
871 prd = np.convolve(prd, zs)
872 return _zseries_to_cseries(prd)
875def chebder(c, m=1, scl=1, axis=0):
876 """
877 Differentiate a Chebyshev series.
879 Returns the Chebyshev series coefficients `c` differentiated `m` times
880 along `axis`. At each iteration the result is multiplied by `scl` (the
881 scaling factor is for use in a linear change of variable). The argument
882 `c` is an array of coefficients from low to high degree along each
883 axis, e.g., [1,2,3] represents the series ``1*T_0 + 2*T_1 + 3*T_2``
884 while [[1,2],[1,2]] represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) +
885 2*T_0(x)*T_1(y) + 2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is
886 ``y``.
888 Parameters
889 ----------
890 c : array_like
891 Array of Chebyshev series coefficients. If c is multidimensional
892 the different axis correspond to different variables with the
893 degree in each axis given by the corresponding index.
894 m : int, optional
895 Number of derivatives taken, must be non-negative. (Default: 1)
896 scl : scalar, optional
897 Each differentiation is multiplied by `scl`. The end result is
898 multiplication by ``scl**m``. This is for use in a linear change of
899 variable. (Default: 1)
900 axis : int, optional
901 Axis over which the derivative is taken. (Default: 0).
903 .. versionadded:: 1.7.0
905 Returns
906 -------
907 der : ndarray
908 Chebyshev series of the derivative.
910 See Also
911 --------
912 chebint
914 Notes
915 -----
916 In general, the result of differentiating a C-series needs to be
917 "reprojected" onto the C-series basis set. Thus, typically, the
918 result of this function is "unintuitive," albeit correct; see Examples
919 section below.
921 Examples
922 --------
923 >>> from numpy.polynomial import chebyshev as C
924 >>> c = (1,2,3,4)
925 >>> C.chebder(c)
926 array([14., 12., 24.])
927 >>> C.chebder(c,3)
928 array([96.])
929 >>> C.chebder(c,scl=-1)
930 array([-14., -12., -24.])
931 >>> C.chebder(c,2,-1)
932 array([12., 96.])
934 """
935 c = np.array(c, ndmin=1, copy=True)
936 if c.dtype.char in '?bBhHiIlLqQpP':
937 c = c.astype(np.double)
938 cnt = pu._deprecate_as_int(m, "the order of derivation")
939 iaxis = pu._deprecate_as_int(axis, "the axis")
940 if cnt < 0:
941 raise ValueError("The order of derivation must be non-negative")
942 iaxis = normalize_axis_index(iaxis, c.ndim)
944 if cnt == 0:
945 return c
947 c = np.moveaxis(c, iaxis, 0)
948 n = len(c)
949 if cnt >= n:
950 c = c[:1]*0
951 else:
952 for i in range(cnt):
953 n = n - 1
954 c *= scl
955 der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
956 for j in range(n, 2, -1):
957 der[j - 1] = (2*j)*c[j]
958 c[j - 2] += (j*c[j])/(j - 2)
959 if n > 1:
960 der[1] = 4*c[2]
961 der[0] = c[1]
962 c = der
963 c = np.moveaxis(c, 0, iaxis)
964 return c
967def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
968 """
969 Integrate a Chebyshev series.
971 Returns the Chebyshev series coefficients `c` integrated `m` times from
972 `lbnd` along `axis`. At each iteration the resulting series is
973 **multiplied** by `scl` and an integration constant, `k`, is added.
974 The scaling factor is for use in a linear change of variable. ("Buyer
975 beware": note that, depending on what one is doing, one may want `scl`
976 to be the reciprocal of what one might expect; for more information,
977 see the Notes section below.) The argument `c` is an array of
978 coefficients from low to high degree along each axis, e.g., [1,2,3]
979 represents the series ``T_0 + 2*T_1 + 3*T_2`` while [[1,2],[1,2]]
980 represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) + 2*T_0(x)*T_1(y) +
981 2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
983 Parameters
984 ----------
985 c : array_like
986 Array of Chebyshev series coefficients. If c is multidimensional
987 the different axis correspond to different variables with the
988 degree in each axis given by the corresponding index.
989 m : int, optional
990 Order of integration, must be positive. (Default: 1)
991 k : {[], list, scalar}, optional
992 Integration constant(s). The value of the first integral at zero
993 is the first value in the list, the value of the second integral
994 at zero is the second value, etc. If ``k == []`` (the default),
995 all constants are set to zero. If ``m == 1``, a single scalar can
996 be given instead of a list.
997 lbnd : scalar, optional
998 The lower bound of the integral. (Default: 0)
999 scl : scalar, optional
1000 Following each integration the result is *multiplied* by `scl`
1001 before the integration constant is added. (Default: 1)
1002 axis : int, optional
1003 Axis over which the integral is taken. (Default: 0).
1005 .. versionadded:: 1.7.0
1007 Returns
1008 -------
1009 S : ndarray
1010 C-series coefficients of the integral.
1012 Raises
1013 ------
1014 ValueError
1015 If ``m < 1``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
1016 ``np.ndim(scl) != 0``.
1018 See Also
1019 --------
1020 chebder
1022 Notes
1023 -----
1024 Note that the result of each integration is *multiplied* by `scl`.
1025 Why is this important to note? Say one is making a linear change of
1026 variable :math:`u = ax + b` in an integral relative to `x`. Then
1027 :math:`dx = du/a`, so one will need to set `scl` equal to
1028 :math:`1/a`- perhaps not what one would have first thought.
1030 Also note that, in general, the result of integrating a C-series needs
1031 to be "reprojected" onto the C-series basis set. Thus, typically,
1032 the result of this function is "unintuitive," albeit correct; see
1033 Examples section below.
1035 Examples
1036 --------
1037 >>> from numpy.polynomial import chebyshev as C
1038 >>> c = (1,2,3)
1039 >>> C.chebint(c)
1040 array([ 0.5, -0.5, 0.5, 0.5])
1041 >>> C.chebint(c,3)
1042 array([ 0.03125 , -0.1875 , 0.04166667, -0.05208333, 0.01041667, # may vary
1043 0.00625 ])
1044 >>> C.chebint(c, k=3)
1045 array([ 3.5, -0.5, 0.5, 0.5])
1046 >>> C.chebint(c,lbnd=-2)
1047 array([ 8.5, -0.5, 0.5, 0.5])
1048 >>> C.chebint(c,scl=-2)
1049 array([-1., 1., -1., -1.])
1051 """
1052 c = np.array(c, ndmin=1, copy=True)
1053 if c.dtype.char in '?bBhHiIlLqQpP':
1054 c = c.astype(np.double)
1055 if not np.iterable(k):
1056 k = [k]
1057 cnt = pu._deprecate_as_int(m, "the order of integration")
1058 iaxis = pu._deprecate_as_int(axis, "the axis")
1059 if cnt < 0:
1060 raise ValueError("The order of integration must be non-negative")
1061 if len(k) > cnt:
1062 raise ValueError("Too many integration constants")
1063 if np.ndim(lbnd) != 0:
1064 raise ValueError("lbnd must be a scalar.")
1065 if np.ndim(scl) != 0:
1066 raise ValueError("scl must be a scalar.")
1067 iaxis = normalize_axis_index(iaxis, c.ndim)
1069 if cnt == 0:
1070 return c
1072 c = np.moveaxis(c, iaxis, 0)
1073 k = list(k) + [0]*(cnt - len(k))
1074 for i in range(cnt):
1075 n = len(c)
1076 c *= scl
1077 if n == 1 and np.all(c[0] == 0):
1078 c[0] += k[i]
1079 else:
1080 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
1081 tmp[0] = c[0]*0
1082 tmp[1] = c[0]
1083 if n > 1:
1084 tmp[2] = c[1]/4
1085 for j in range(2, n):
1086 tmp[j + 1] = c[j]/(2*(j + 1))
1087 tmp[j - 1] -= c[j]/(2*(j - 1))
1088 tmp[0] += k[i] - chebval(lbnd, tmp)
1089 c = tmp
1090 c = np.moveaxis(c, 0, iaxis)
1091 return c
1094def chebval(x, c, tensor=True):
1095 """
1096 Evaluate a Chebyshev series at points x.
1098 If `c` is of length `n + 1`, this function returns the value:
1100 .. math:: p(x) = c_0 * T_0(x) + c_1 * T_1(x) + ... + c_n * T_n(x)
1102 The parameter `x` is converted to an array only if it is a tuple or a
1103 list, otherwise it is treated as a scalar. In either case, either `x`
1104 or its elements must support multiplication and addition both with
1105 themselves and with the elements of `c`.
1107 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
1108 `c` is multidimensional, then the shape of the result depends on the
1109 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
1110 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
1111 scalars have shape (,).
1113 Trailing zeros in the coefficients will be used in the evaluation, so
1114 they should be avoided if efficiency is a concern.
1116 Parameters
1117 ----------
1118 x : array_like, compatible object
1119 If `x` is a list or tuple, it is converted to an ndarray, otherwise
1120 it is left unchanged and treated as a scalar. In either case, `x`
1121 or its elements must support addition and multiplication with
1122 themselves and with the elements of `c`.
1123 c : array_like
1124 Array of coefficients ordered so that the coefficients for terms of
1125 degree n are contained in c[n]. If `c` is multidimensional the
1126 remaining indices enumerate multiple polynomials. In the two
1127 dimensional case the coefficients may be thought of as stored in
1128 the columns of `c`.
1129 tensor : boolean, optional
1130 If True, the shape of the coefficient array is extended with ones
1131 on the right, one for each dimension of `x`. Scalars have dimension 0
1132 for this action. The result is that every column of coefficients in
1133 `c` is evaluated for every element of `x`. If False, `x` is broadcast
1134 over the columns of `c` for the evaluation. This keyword is useful
1135 when `c` is multidimensional. The default value is True.
1137 .. versionadded:: 1.7.0
1139 Returns
1140 -------
1141 values : ndarray, algebra_like
1142 The shape of the return value is described above.
1144 See Also
1145 --------
1146 chebval2d, chebgrid2d, chebval3d, chebgrid3d
1148 Notes
1149 -----
1150 The evaluation uses Clenshaw recursion, aka synthetic division.
1152 """
1153 c = np.array(c, ndmin=1, copy=True)
1154 if c.dtype.char in '?bBhHiIlLqQpP':
1155 c = c.astype(np.double)
1156 if isinstance(x, (tuple, list)):
1157 x = np.asarray(x)
1158 if isinstance(x, np.ndarray) and tensor:
1159 c = c.reshape(c.shape + (1,)*x.ndim)
1161 if len(c) == 1:
1162 c0 = c[0]
1163 c1 = 0
1164 elif len(c) == 2:
1165 c0 = c[0]
1166 c1 = c[1]
1167 else:
1168 x2 = 2*x
1169 c0 = c[-2]
1170 c1 = c[-1]
1171 for i in range(3, len(c) + 1):
1172 tmp = c0
1173 c0 = c[-i] - c1
1174 c1 = tmp + c1*x2
1175 return c0 + c1*x
1178def chebval2d(x, y, c):
1179 """
1180 Evaluate a 2-D Chebyshev series at points (x, y).
1182 This function returns the values:
1184 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * T_i(x) * T_j(y)
1186 The parameters `x` and `y` are converted to arrays only if they are
1187 tuples or a lists, otherwise they are treated as a scalars and they
1188 must have the same shape after conversion. In either case, either `x`
1189 and `y` or their elements must support multiplication and addition both
1190 with themselves and with the elements of `c`.
1192 If `c` is a 1-D array a one is implicitly appended to its shape to make
1193 it 2-D. The shape of the result will be c.shape[2:] + x.shape.
1195 Parameters
1196 ----------
1197 x, y : array_like, compatible objects
1198 The two dimensional series is evaluated at the points `(x, y)`,
1199 where `x` and `y` must have the same shape. If `x` or `y` is a list
1200 or tuple, it is first converted to an ndarray, otherwise it is left
1201 unchanged and if it isn't an ndarray it is treated as a scalar.
1202 c : array_like
1203 Array of coefficients ordered so that the coefficient of the term
1204 of multi-degree i,j is contained in ``c[i,j]``. If `c` has
1205 dimension greater than 2 the remaining indices enumerate multiple
1206 sets of coefficients.
1208 Returns
1209 -------
1210 values : ndarray, compatible object
1211 The values of the two dimensional Chebyshev series at points formed
1212 from pairs of corresponding values from `x` and `y`.
1214 See Also
1215 --------
1216 chebval, chebgrid2d, chebval3d, chebgrid3d
1218 Notes
1219 -----
1221 .. versionadded:: 1.7.0
1223 """
1224 return pu._valnd(chebval, c, x, y)
1227def chebgrid2d(x, y, c):
1228 """
1229 Evaluate a 2-D Chebyshev series on the Cartesian product of x and y.
1231 This function returns the values:
1233 .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * T_i(a) * T_j(b),
1235 where the points `(a, b)` consist of all pairs formed by taking
1236 `a` from `x` and `b` from `y`. The resulting points form a grid with
1237 `x` in the first dimension and `y` in the second.
1239 The parameters `x` and `y` are converted to arrays only if they are
1240 tuples or a lists, otherwise they are treated as a scalars. In either
1241 case, either `x` and `y` or their elements must support multiplication
1242 and addition both with themselves and with the elements of `c`.
1244 If `c` has fewer than two dimensions, ones are implicitly appended to
1245 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
1246 x.shape + y.shape.
1248 Parameters
1249 ----------
1250 x, y : array_like, compatible objects
1251 The two dimensional series is evaluated at the points in the
1252 Cartesian product of `x` and `y`. If `x` or `y` is a list or
1253 tuple, it is first converted to an ndarray, otherwise it is left
1254 unchanged and, if it isn't an ndarray, it is treated as a scalar.
1255 c : array_like
1256 Array of coefficients ordered so that the coefficient of the term of
1257 multi-degree i,j is contained in `c[i,j]`. If `c` has dimension
1258 greater than two the remaining indices enumerate multiple sets of
1259 coefficients.
1261 Returns
1262 -------
1263 values : ndarray, compatible object
1264 The values of the two dimensional Chebyshev series at points in the
1265 Cartesian product of `x` and `y`.
1267 See Also
1268 --------
1269 chebval, chebval2d, chebval3d, chebgrid3d
1271 Notes
1272 -----
1274 .. versionadded:: 1.7.0
1276 """
1277 return pu._gridnd(chebval, c, x, y)
1280def chebval3d(x, y, z, c):
1281 """
1282 Evaluate a 3-D Chebyshev series at points (x, y, z).
1284 This function returns the values:
1286 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * T_i(x) * T_j(y) * T_k(z)
1288 The parameters `x`, `y`, and `z` are converted to arrays only if
1289 they are tuples or a lists, otherwise they are treated as a scalars and
1290 they must have the same shape after conversion. In either case, either
1291 `x`, `y`, and `z` or their elements must support multiplication and
1292 addition both with themselves and with the elements of `c`.
1294 If `c` has fewer than 3 dimensions, ones are implicitly appended to its
1295 shape to make it 3-D. The shape of the result will be c.shape[3:] +
1296 x.shape.
1298 Parameters
1299 ----------
1300 x, y, z : array_like, compatible object
1301 The three dimensional series is evaluated at the points
1302 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
1303 any of `x`, `y`, or `z` is a list or tuple, it is first converted
1304 to an ndarray, otherwise it is left unchanged and if it isn't an
1305 ndarray it is treated as a scalar.
1306 c : array_like
1307 Array of coefficients ordered so that the coefficient of the term of
1308 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
1309 greater than 3 the remaining indices enumerate multiple sets of
1310 coefficients.
1312 Returns
1313 -------
1314 values : ndarray, compatible object
1315 The values of the multidimensional polynomial on points formed with
1316 triples of corresponding values from `x`, `y`, and `z`.
1318 See Also
1319 --------
1320 chebval, chebval2d, chebgrid2d, chebgrid3d
1322 Notes
1323 -----
1325 .. versionadded:: 1.7.0
1327 """
1328 return pu._valnd(chebval, c, x, y, z)
1331def chebgrid3d(x, y, z, c):
1332 """
1333 Evaluate a 3-D Chebyshev series on the Cartesian product of x, y, and z.
1335 This function returns the values:
1337 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * T_i(a) * T_j(b) * T_k(c)
1339 where the points `(a, b, c)` consist of all triples formed by taking
1340 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
1341 a grid with `x` in the first dimension, `y` in the second, and `z` in
1342 the third.
1344 The parameters `x`, `y`, and `z` are converted to arrays only if they
1345 are tuples or a lists, otherwise they are treated as a scalars. In
1346 either case, either `x`, `y`, and `z` or their elements must support
1347 multiplication and addition both with themselves and with the elements
1348 of `c`.
1350 If `c` has fewer than three dimensions, ones are implicitly appended to
1351 its shape to make it 3-D. The shape of the result will be c.shape[3:] +
1352 x.shape + y.shape + z.shape.
1354 Parameters
1355 ----------
1356 x, y, z : array_like, compatible objects
1357 The three dimensional series is evaluated at the points in the
1358 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
1359 list or tuple, it is first converted to an ndarray, otherwise it is
1360 left unchanged and, if it isn't an ndarray, it is treated as a
1361 scalar.
1362 c : array_like
1363 Array of coefficients ordered so that the coefficients for terms of
1364 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1365 greater than two the remaining indices enumerate multiple sets of
1366 coefficients.
1368 Returns
1369 -------
1370 values : ndarray, compatible object
1371 The values of the two dimensional polynomial at points in the Cartesian
1372 product of `x` and `y`.
1374 See Also
1375 --------
1376 chebval, chebval2d, chebgrid2d, chebval3d
1378 Notes
1379 -----
1381 .. versionadded:: 1.7.0
1383 """
1384 return pu._gridnd(chebval, c, x, y, z)
1387def chebvander(x, deg):
1388 """Pseudo-Vandermonde matrix of given degree.
1390 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
1391 `x`. The pseudo-Vandermonde matrix is defined by
1393 .. math:: V[..., i] = T_i(x),
1395 where `0 <= i <= deg`. The leading indices of `V` index the elements of
1396 `x` and the last index is the degree of the Chebyshev polynomial.
1398 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
1399 matrix ``V = chebvander(x, n)``, then ``np.dot(V, c)`` and
1400 ``chebval(x, c)`` are the same up to roundoff. This equivalence is
1401 useful both for least squares fitting and for the evaluation of a large
1402 number of Chebyshev series of the same degree and sample points.
1404 Parameters
1405 ----------
1406 x : array_like
1407 Array of points. The dtype is converted to float64 or complex128
1408 depending on whether any of the elements are complex. If `x` is
1409 scalar it is converted to a 1-D array.
1410 deg : int
1411 Degree of the resulting matrix.
1413 Returns
1414 -------
1415 vander : ndarray
1416 The pseudo Vandermonde matrix. The shape of the returned matrix is
1417 ``x.shape + (deg + 1,)``, where The last index is the degree of the
1418 corresponding Chebyshev polynomial. The dtype will be the same as
1419 the converted `x`.
1421 """
1422 ideg = pu._deprecate_as_int(deg, "deg")
1423 if ideg < 0:
1424 raise ValueError("deg must be non-negative")
1426 x = np.array(x, copy=False, ndmin=1) + 0.0
1427 dims = (ideg + 1,) + x.shape
1428 dtyp = x.dtype
1429 v = np.empty(dims, dtype=dtyp)
1430 # Use forward recursion to generate the entries.
1431 v[0] = x*0 + 1
1432 if ideg > 0:
1433 x2 = 2*x
1434 v[1] = x
1435 for i in range(2, ideg + 1):
1436 v[i] = v[i-1]*x2 - v[i-2]
1437 return np.moveaxis(v, 0, -1)
1440def chebvander2d(x, y, deg):
1441 """Pseudo-Vandermonde matrix of given degrees.
1443 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1444 points `(x, y)`. The pseudo-Vandermonde matrix is defined by
1446 .. math:: V[..., (deg[1] + 1)*i + j] = T_i(x) * T_j(y),
1448 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
1449 `V` index the points `(x, y)` and the last index encodes the degrees of
1450 the Chebyshev polynomials.
1452 If ``V = chebvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
1453 correspond to the elements of a 2-D coefficient array `c` of shape
1454 (xdeg + 1, ydeg + 1) in the order
1456 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
1458 and ``np.dot(V, c.flat)`` and ``chebval2d(x, y, c)`` will be the same
1459 up to roundoff. This equivalence is useful both for least squares
1460 fitting and for the evaluation of a large number of 2-D Chebyshev
1461 series of the same degrees and sample points.
1463 Parameters
1464 ----------
1465 x, y : array_like
1466 Arrays of point coordinates, all of the same shape. The dtypes
1467 will be converted to either float64 or complex128 depending on
1468 whether any of the elements are complex. Scalars are converted to
1469 1-D arrays.
1470 deg : list of ints
1471 List of maximum degrees of the form [x_deg, y_deg].
1473 Returns
1474 -------
1475 vander2d : ndarray
1476 The shape of the returned matrix is ``x.shape + (order,)``, where
1477 :math:`order = (deg[0]+1)*(deg[1]+1)`. The dtype will be the same
1478 as the converted `x` and `y`.
1480 See Also
1481 --------
1482 chebvander, chebvander3d, chebval2d, chebval3d
1484 Notes
1485 -----
1487 .. versionadded:: 1.7.0
1489 """
1490 return pu._vander_nd_flat((chebvander, chebvander), (x, y), deg)
1493def chebvander3d(x, y, z, deg):
1494 """Pseudo-Vandermonde matrix of given degrees.
1496 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1497 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
1498 then The pseudo-Vandermonde matrix is defined by
1500 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = T_i(x)*T_j(y)*T_k(z),
1502 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
1503 indices of `V` index the points `(x, y, z)` and the last index encodes
1504 the degrees of the Chebyshev polynomials.
1506 If ``V = chebvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
1507 of `V` correspond to the elements of a 3-D coefficient array `c` of
1508 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
1510 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
1512 and ``np.dot(V, c.flat)`` and ``chebval3d(x, y, z, c)`` will be the
1513 same up to roundoff. This equivalence is useful both for least squares
1514 fitting and for the evaluation of a large number of 3-D Chebyshev
1515 series of the same degrees and sample points.
1517 Parameters
1518 ----------
1519 x, y, z : array_like
1520 Arrays of point coordinates, all of the same shape. The dtypes will
1521 be converted to either float64 or complex128 depending on whether
1522 any of the elements are complex. Scalars are converted to 1-D
1523 arrays.
1524 deg : list of ints
1525 List of maximum degrees of the form [x_deg, y_deg, z_deg].
1527 Returns
1528 -------
1529 vander3d : ndarray
1530 The shape of the returned matrix is ``x.shape + (order,)``, where
1531 :math:`order = (deg[0]+1)*(deg[1]+1)*(deg[2]+1)`. The dtype will
1532 be the same as the converted `x`, `y`, and `z`.
1534 See Also
1535 --------
1536 chebvander, chebvander3d, chebval2d, chebval3d
1538 Notes
1539 -----
1541 .. versionadded:: 1.7.0
1543 """
1544 return pu._vander_nd_flat((chebvander, chebvander, chebvander), (x, y, z), deg)
1547def chebfit(x, y, deg, rcond=None, full=False, w=None):
1548 """
1549 Least squares fit of Chebyshev series to data.
1551 Return the coefficients of a Chebyshev series of degree `deg` that is the
1552 least squares fit to the data values `y` given at points `x`. If `y` is
1553 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
1554 fits are done, one for each column of `y`, and the resulting
1555 coefficients are stored in the corresponding columns of a 2-D return.
1556 The fitted polynomial(s) are in the form
1558 .. math:: p(x) = c_0 + c_1 * T_1(x) + ... + c_n * T_n(x),
1560 where `n` is `deg`.
1562 Parameters
1563 ----------
1564 x : array_like, shape (M,)
1565 x-coordinates of the M sample points ``(x[i], y[i])``.
1566 y : array_like, shape (M,) or (M, K)
1567 y-coordinates of the sample points. Several data sets of sample
1568 points sharing the same x-coordinates can be fitted at once by
1569 passing in a 2D-array that contains one dataset per column.
1570 deg : int or 1-D array_like
1571 Degree(s) of the fitting polynomials. If `deg` is a single integer,
1572 all terms up to and including the `deg`'th term are included in the
1573 fit. For NumPy versions >= 1.11.0 a list of integers specifying the
1574 degrees of the terms to include may be used instead.
1575 rcond : float, optional
1576 Relative condition number of the fit. Singular values smaller than
1577 this relative to the largest singular value will be ignored. The
1578 default value is len(x)*eps, where eps is the relative precision of
1579 the float type, about 2e-16 in most cases.
1580 full : bool, optional
1581 Switch determining nature of return value. When it is False (the
1582 default) just the coefficients are returned, when True diagnostic
1583 information from the singular value decomposition is also returned.
1584 w : array_like, shape (`M`,), optional
1585 Weights. If not None, the weight ``w[i]`` applies to the unsquared
1586 residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
1587 chosen so that the errors of the products ``w[i]*y[i]`` all have the
1588 same variance. When using inverse-variance weighting, use
1589 ``w[i] = 1/sigma(y[i])``. The default value is None.
1591 .. versionadded:: 1.5.0
1593 Returns
1594 -------
1595 coef : ndarray, shape (M,) or (M, K)
1596 Chebyshev coefficients ordered from low to high. If `y` was 2-D,
1597 the coefficients for the data in column k of `y` are in column
1598 `k`.
1600 [residuals, rank, singular_values, rcond] : list
1601 These values are only returned if ``full == True``
1603 - residuals -- sum of squared residuals of the least squares fit
1604 - rank -- the numerical rank of the scaled Vandermonde matrix
1605 - singular_values -- singular values of the scaled Vandermonde matrix
1606 - rcond -- value of `rcond`.
1608 For more details, see `numpy.linalg.lstsq`.
1610 Warns
1611 -----
1612 RankWarning
1613 The rank of the coefficient matrix in the least-squares fit is
1614 deficient. The warning is only raised if ``full == False``. The
1615 warnings can be turned off by
1617 >>> import warnings
1618 >>> warnings.simplefilter('ignore', np.RankWarning)
1620 See Also
1621 --------
1622 numpy.polynomial.polynomial.polyfit
1623 numpy.polynomial.legendre.legfit
1624 numpy.polynomial.laguerre.lagfit
1625 numpy.polynomial.hermite.hermfit
1626 numpy.polynomial.hermite_e.hermefit
1627 chebval : Evaluates a Chebyshev series.
1628 chebvander : Vandermonde matrix of Chebyshev series.
1629 chebweight : Chebyshev weight function.
1630 numpy.linalg.lstsq : Computes a least-squares fit from the matrix.
1631 scipy.interpolate.UnivariateSpline : Computes spline fits.
1633 Notes
1634 -----
1635 The solution is the coefficients of the Chebyshev series `p` that
1636 minimizes the sum of the weighted squared errors
1638 .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
1640 where :math:`w_j` are the weights. This problem is solved by setting up
1641 as the (typically) overdetermined matrix equation
1643 .. math:: V(x) * c = w * y,
1645 where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
1646 coefficients to be solved for, `w` are the weights, and `y` are the
1647 observed values. This equation is then solved using the singular value
1648 decomposition of `V`.
1650 If some of the singular values of `V` are so small that they are
1651 neglected, then a `RankWarning` will be issued. This means that the
1652 coefficient values may be poorly determined. Using a lower order fit
1653 will usually get rid of the warning. The `rcond` parameter can also be
1654 set to a value smaller than its default, but the resulting fit may be
1655 spurious and have large contributions from roundoff error.
1657 Fits using Chebyshev series are usually better conditioned than fits
1658 using power series, but much can depend on the distribution of the
1659 sample points and the smoothness of the data. If the quality of the fit
1660 is inadequate splines may be a good alternative.
1662 References
1663 ----------
1664 .. [1] Wikipedia, "Curve fitting",
1665 https://en.wikipedia.org/wiki/Curve_fitting
1667 Examples
1668 --------
1670 """
1671 return pu._fit(chebvander, x, y, deg, rcond, full, w)
1674def chebcompanion(c):
1675 """Return the scaled companion matrix of c.
1677 The basis polynomials are scaled so that the companion matrix is
1678 symmetric when `c` is a Chebyshev basis polynomial. This provides
1679 better eigenvalue estimates than the unscaled case and for basis
1680 polynomials the eigenvalues are guaranteed to be real if
1681 `numpy.linalg.eigvalsh` is used to obtain them.
1683 Parameters
1684 ----------
1685 c : array_like
1686 1-D array of Chebyshev series coefficients ordered from low to high
1687 degree.
1689 Returns
1690 -------
1691 mat : ndarray
1692 Scaled companion matrix of dimensions (deg, deg).
1694 Notes
1695 -----
1697 .. versionadded:: 1.7.0
1699 """
1700 # c is a trimmed copy
1701 [c] = pu.as_series([c])
1702 if len(c) < 2:
1703 raise ValueError('Series must have maximum degree of at least 1.')
1704 if len(c) == 2:
1705 return np.array([[-c[0]/c[1]]])
1707 n = len(c) - 1
1708 mat = np.zeros((n, n), dtype=c.dtype)
1709 scl = np.array([1.] + [np.sqrt(.5)]*(n-1))
1710 top = mat.reshape(-1)[1::n+1]
1711 bot = mat.reshape(-1)[n::n+1]
1712 top[0] = np.sqrt(.5)
1713 top[1:] = 1/2
1714 bot[...] = top
1715 mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*.5
1716 return mat
1719def chebroots(c):
1720 """
1721 Compute the roots of a Chebyshev series.
1723 Return the roots (a.k.a. "zeros") of the polynomial
1725 .. math:: p(x) = \\sum_i c[i] * T_i(x).
1727 Parameters
1728 ----------
1729 c : 1-D array_like
1730 1-D array of coefficients.
1732 Returns
1733 -------
1734 out : ndarray
1735 Array of the roots of the series. If all the roots are real,
1736 then `out` is also real, otherwise it is complex.
1738 See Also
1739 --------
1740 numpy.polynomial.polynomial.polyroots
1741 numpy.polynomial.legendre.legroots
1742 numpy.polynomial.laguerre.lagroots
1743 numpy.polynomial.hermite.hermroots
1744 numpy.polynomial.hermite_e.hermeroots
1746 Notes
1747 -----
1748 The root estimates are obtained as the eigenvalues of the companion
1749 matrix, Roots far from the origin of the complex plane may have large
1750 errors due to the numerical instability of the series for such
1751 values. Roots with multiplicity greater than 1 will also show larger
1752 errors as the value of the series near such points is relatively
1753 insensitive to errors in the roots. Isolated roots near the origin can
1754 be improved by a few iterations of Newton's method.
1756 The Chebyshev series basis polynomials aren't powers of `x` so the
1757 results of this function may seem unintuitive.
1759 Examples
1760 --------
1761 >>> import numpy.polynomial.chebyshev as cheb
1762 >>> cheb.chebroots((-1, 1,-1, 1)) # T3 - T2 + T1 - T0 has real roots
1763 array([ -5.00000000e-01, 2.60860684e-17, 1.00000000e+00]) # may vary
1765 """
1766 # c is a trimmed copy
1767 [c] = pu.as_series([c])
1768 if len(c) < 2:
1769 return np.array([], dtype=c.dtype)
1770 if len(c) == 2:
1771 return np.array([-c[0]/c[1]])
1773 # rotated companion matrix reduces error
1774 m = chebcompanion(c)[::-1,::-1]
1775 r = la.eigvals(m)
1776 r.sort()
1777 return r
1780def chebinterpolate(func, deg, args=()):
1781 """Interpolate a function at the Chebyshev points of the first kind.
1783 Returns the Chebyshev series that interpolates `func` at the Chebyshev
1784 points of the first kind in the interval [-1, 1]. The interpolating
1785 series tends to a minmax approximation to `func` with increasing `deg`
1786 if the function is continuous in the interval.
1788 .. versionadded:: 1.14.0
1790 Parameters
1791 ----------
1792 func : function
1793 The function to be approximated. It must be a function of a single
1794 variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are
1795 extra arguments passed in the `args` parameter.
1796 deg : int
1797 Degree of the interpolating polynomial
1798 args : tuple, optional
1799 Extra arguments to be used in the function call. Default is no extra
1800 arguments.
1802 Returns
1803 -------
1804 coef : ndarray, shape (deg + 1,)
1805 Chebyshev coefficients of the interpolating series ordered from low to
1806 high.
1808 Examples
1809 --------
1810 >>> import numpy.polynomial.chebyshev as C
1811 >>> C.chebfromfunction(lambda x: np.tanh(x) + 0.5, 8)
1812 array([ 5.00000000e-01, 8.11675684e-01, -9.86864911e-17,
1813 -5.42457905e-02, -2.71387850e-16, 4.51658839e-03,
1814 2.46716228e-17, -3.79694221e-04, -3.26899002e-16])
1816 Notes
1817 -----
1819 The Chebyshev polynomials used in the interpolation are orthogonal when
1820 sampled at the Chebyshev points of the first kind. If it is desired to
1821 constrain some of the coefficients they can simply be set to the desired
1822 value after the interpolation, no new interpolation or fit is needed. This
1823 is especially useful if it is known apriori that some of coefficients are
1824 zero. For instance, if the function is even then the coefficients of the
1825 terms of odd degree in the result can be set to zero.
1827 """
1828 deg = np.asarray(deg)
1830 # check arguments.
1831 if deg.ndim > 0 or deg.dtype.kind not in 'iu' or deg.size == 0:
1832 raise TypeError("deg must be an int")
1833 if deg < 0:
1834 raise ValueError("expected deg >= 0")
1836 order = deg + 1
1837 xcheb = chebpts1(order)
1838 yfunc = func(xcheb, *args)
1839 m = chebvander(xcheb, deg)
1840 c = np.dot(m.T, yfunc)
1841 c[0] /= order
1842 c[1:] /= 0.5*order
1844 return c
1847def chebgauss(deg):
1848 """
1849 Gauss-Chebyshev quadrature.
1851 Computes the sample points and weights for Gauss-Chebyshev quadrature.
1852 These sample points and weights will correctly integrate polynomials of
1853 degree :math:`2*deg - 1` or less over the interval :math:`[-1, 1]` with
1854 the weight function :math:`f(x) = 1/\\sqrt{1 - x^2}`.
1856 Parameters
1857 ----------
1858 deg : int
1859 Number of sample points and weights. It must be >= 1.
1861 Returns
1862 -------
1863 x : ndarray
1864 1-D ndarray containing the sample points.
1865 y : ndarray
1866 1-D ndarray containing the weights.
1868 Notes
1869 -----
1871 .. versionadded:: 1.7.0
1873 The results have only been tested up to degree 100, higher degrees may
1874 be problematic. For Gauss-Chebyshev there are closed form solutions for
1875 the sample points and weights. If n = `deg`, then
1877 .. math:: x_i = \\cos(\\pi (2 i - 1) / (2 n))
1879 .. math:: w_i = \\pi / n
1881 """
1882 ideg = pu._deprecate_as_int(deg, "deg")
1883 if ideg <= 0:
1884 raise ValueError("deg must be a positive integer")
1886 x = np.cos(np.pi * np.arange(1, 2*ideg, 2) / (2.0*ideg))
1887 w = np.ones(ideg)*(np.pi/ideg)
1889 return x, w
1892def chebweight(x):
1893 """
1894 The weight function of the Chebyshev polynomials.
1896 The weight function is :math:`1/\\sqrt{1 - x^2}` and the interval of
1897 integration is :math:`[-1, 1]`. The Chebyshev polynomials are
1898 orthogonal, but not normalized, with respect to this weight function.
1900 Parameters
1901 ----------
1902 x : array_like
1903 Values at which the weight function will be computed.
1905 Returns
1906 -------
1907 w : ndarray
1908 The weight function at `x`.
1910 Notes
1911 -----
1913 .. versionadded:: 1.7.0
1915 """
1916 w = 1./(np.sqrt(1. + x) * np.sqrt(1. - x))
1917 return w
1920def chebpts1(npts):
1921 """
1922 Chebyshev points of the first kind.
1924 The Chebyshev points of the first kind are the points ``cos(x)``,
1925 where ``x = [pi*(k + .5)/npts for k in range(npts)]``.
1927 Parameters
1928 ----------
1929 npts : int
1930 Number of sample points desired.
1932 Returns
1933 -------
1934 pts : ndarray
1935 The Chebyshev points of the first kind.
1937 See Also
1938 --------
1939 chebpts2
1941 Notes
1942 -----
1944 .. versionadded:: 1.5.0
1946 """
1947 _npts = int(npts)
1948 if _npts != npts:
1949 raise ValueError("npts must be integer")
1950 if _npts < 1:
1951 raise ValueError("npts must be >= 1")
1953 x = 0.5 * np.pi / _npts * np.arange(-_npts+1, _npts+1, 2)
1954 return np.sin(x)
1957def chebpts2(npts):
1958 """
1959 Chebyshev points of the second kind.
1961 The Chebyshev points of the second kind are the points ``cos(x)``,
1962 where ``x = [pi*k/(npts - 1) for k in range(npts)]`` sorted in ascending
1963 order.
1965 Parameters
1966 ----------
1967 npts : int
1968 Number of sample points desired.
1970 Returns
1971 -------
1972 pts : ndarray
1973 The Chebyshev points of the second kind.
1975 Notes
1976 -----
1978 .. versionadded:: 1.5.0
1980 """
1981 _npts = int(npts)
1982 if _npts != npts:
1983 raise ValueError("npts must be integer")
1984 if _npts < 2:
1985 raise ValueError("npts must be >= 2")
1987 x = np.linspace(-np.pi, 0, _npts)
1988 return np.cos(x)
1991#
1992# Chebyshev series class
1993#
1995class Chebyshev(ABCPolyBase):
1996 """A Chebyshev series class.
1998 The Chebyshev class provides the standard Python numerical methods
1999 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
2000 methods listed below.
2002 Parameters
2003 ----------
2004 coef : array_like
2005 Chebyshev coefficients in order of increasing degree, i.e.,
2006 ``(1, 2, 3)`` gives ``1*T_0(x) + 2*T_1(x) + 3*T_2(x)``.
2007 domain : (2,) array_like, optional
2008 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
2009 to the interval ``[window[0], window[1]]`` by shifting and scaling.
2010 The default value is [-1, 1].
2011 window : (2,) array_like, optional
2012 Window, see `domain` for its use. The default value is [-1, 1].
2014 .. versionadded:: 1.6.0
2016 """
2017 # Virtual Functions
2018 _add = staticmethod(chebadd)
2019 _sub = staticmethod(chebsub)
2020 _mul = staticmethod(chebmul)
2021 _div = staticmethod(chebdiv)
2022 _pow = staticmethod(chebpow)
2023 _val = staticmethod(chebval)
2024 _int = staticmethod(chebint)
2025 _der = staticmethod(chebder)
2026 _fit = staticmethod(chebfit)
2027 _line = staticmethod(chebline)
2028 _roots = staticmethod(chebroots)
2029 _fromroots = staticmethod(chebfromroots)
2031 @classmethod
2032 def interpolate(cls, func, deg, domain=None, args=()):
2033 """Interpolate a function at the Chebyshev points of the first kind.
2035 Returns the series that interpolates `func` at the Chebyshev points of
2036 the first kind scaled and shifted to the `domain`. The resulting series
2037 tends to a minmax approximation of `func` when the function is
2038 continuous in the domain.
2040 .. versionadded:: 1.14.0
2042 Parameters
2043 ----------
2044 func : function
2045 The function to be interpolated. It must be a function of a single
2046 variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are
2047 extra arguments passed in the `args` parameter.
2048 deg : int
2049 Degree of the interpolating polynomial.
2050 domain : {None, [beg, end]}, optional
2051 Domain over which `func` is interpolated. The default is None, in
2052 which case the domain is [-1, 1].
2053 args : tuple, optional
2054 Extra arguments to be used in the function call. Default is no
2055 extra arguments.
2057 Returns
2058 -------
2059 polynomial : Chebyshev instance
2060 Interpolating Chebyshev instance.
2062 Notes
2063 -----
2064 See `numpy.polynomial.chebfromfunction` for more details.
2066 """
2067 if domain is None:
2068 domain = cls.domain
2069 xfunc = lambda x: func(pu.mapdomain(x, cls.window, domain), *args)
2070 coef = chebinterpolate(xfunc, deg)
2071 return cls(coef, domain=domain)
2073 # Virtual properties
2074 domain = np.array(chebdomain)
2075 window = np.array(chebdomain)
2076 basis_name = 'T'