codekingpro/portable-devtools
114k
1"""
2====================================================
3Chebyshev Series (:mod:`numpy.polynomial.chebyshev`)
4====================================================
5
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`).
11
12Classes
13-------
14
15.. autosummary::
16 :toctree: generated/
17
18 Chebyshev
19
20
21Constants
22---------
23
24.. autosummary::
25 :toctree: generated/
26
27 chebdomain
28 chebzero
29 chebone
30 chebx
31
32Arithmetic
33----------
34
35.. autosummary::
36 :toctree: generated/
37
38 chebadd
39 chebsub
40 chebmulx
41 chebmul
42 chebdiv
43 chebpow
44 chebval
45 chebval2d
46 chebval3d
47 chebgrid2d
48 chebgrid3d
49
50Calculus
51--------
52
53.. autosummary::
54 :toctree: generated/
55
56 chebder
57 chebint
58
59Misc Functions
60--------------
61
62.. autosummary::
63 :toctree: generated/
64
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
81
82See also
83--------
84`numpy.polynomial`
85
86Notes
87-----
88The implementations of multiplication, division, integration, and
89differentiation use the algebraic identities [1]_:
90
91.. math::
92 T_n(x) = \\frac{z^n + z^{-n}}{2} \\\\
93 z\\frac{dx}{dz} = \\frac{z - z^{-1}}{2}.
94
95where
96
97.. math:: x = \\frac{z + z^{-1}}{2}.
98
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."
102
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)
108
109""" # noqa: E501
110import numpy as np
111
112from . import polyutils as pu
113from ._polybase import ABCPolyBase
114
115__all__ = [
116 'chebzero', 'chebone', 'chebx', 'chebdomain', 'chebline', 'chebadd',
117 'chebsub', 'chebmulx', 'chebmul', 'chebdiv', 'chebpow', 'chebval',
118 'chebder', 'chebint', 'cheb2poly', 'poly2cheb', 'chebfromroots',
119 'chebvander', 'chebfit', 'chebtrim', 'chebroots', 'chebpts1',
120 'chebpts2', 'Chebyshev', 'chebval2d', 'chebval3d', 'chebgrid2d',
121 'chebgrid3d', 'chebvander2d', 'chebvander3d', 'chebcompanion',
122 'chebgauss', 'chebweight', 'chebinterpolate']
123
124chebtrim = pu.trimcoef
125
126#
127# A collection of functions for manipulating z-series. These are private
128# functions and do minimal error checking.
129#
130
131def _cseries_to_zseries(c):
132 """Convert Chebyshev series to z-series.
133
134 Convert a Chebyshev series to the equivalent z-series. The result is
135 never an empty array. The dtype of the return is the same as that of
136 the input. No checks are run on the arguments as this routine is for
137 internal use.
138
139 Parameters
140 ----------
141 c : 1-D ndarray
142 Chebyshev coefficients, ordered from low to high
143
144 Returns
145 -------
146 zs : 1-D ndarray
147 Odd length symmetric z-series, ordered from low to high.
148
149 """
150 n = c.size
151 zs = np.zeros(2 * n - 1, dtype=c.dtype)
152 zs[n - 1:] = c / 2
153 return zs + zs[::-1]
154
155
156def _zseries_to_cseries(zs):
157 """Convert z-series to a Chebyshev series.
158
159 Convert a z series to the equivalent Chebyshev series. The result is
160 never an empty array. The dtype of the return is the same as that of
161 the input. No checks are run on the arguments as this routine is for
162 internal use.
163
164 Parameters
165 ----------
166 zs : 1-D ndarray
167 Odd length symmetric z-series, ordered from low to high.
168
169 Returns
170 -------
171 c : 1-D ndarray
172 Chebyshev coefficients, ordered from low to high.
173
174 """
175 n = (zs.size + 1) // 2
176 c = zs[n - 1:].copy()
177 c[1:n] *= 2
178 return c
179
180
181def _zseries_mul(z1, z2):
182 """Multiply two z-series.
183
184 Multiply two z-series to produce a z-series.
185
186 Parameters
187 ----------
188 z1, z2 : 1-D ndarray
189 The arrays must be 1-D but this is not checked.
190
191 Returns
192 -------
193 product : 1-D ndarray
194 The product z-series.
195
196 Notes
197 -----
198 This is simply convolution. If symmetric/anti-symmetric z-series are
199 denoted by S/A then the following rules apply:
200
201 S*S, A*A -> S
202 S*A, A*S -> A
203
204 """
205 return np.convolve(z1, z2)
206
207
208def _zseries_div(z1, z2):
209 """Divide the first z-series by the second.
210
211 Divide `z1` by `z2` and return the quotient and remainder as z-series.
212 Warning: this implementation only applies when both z1 and z2 have the
213 same symmetry, which is sufficient for present purposes.
214
215 Parameters
216 ----------
217 z1, z2 : 1-D ndarray
218 The arrays must be 1-D and have the same symmetry, but this is not
219 checked.
220
221 Returns
222 -------
223
224 (quotient, remainder) : 1-D ndarrays
225 Quotient and remainder as z-series.
226
227 Notes
228 -----
229 This is not the same as polynomial division on account of the desired form
230 of the remainder. If symmetric/anti-symmetric z-series are denoted by S/A
231 then the following rules apply:
232
233 S/S -> S,S
234 A/A -> S,A
235
236 The restriction to types of the same symmetry could be fixed but seems like
237 unneeded generality. There is no natural form for the remainder in the case
238 where there is no symmetry.
239
240 """
241 z1 = z1.copy()
242 z2 = z2.copy()
243 lc1 = len(z1)
244 lc2 = len(z2)
245 if lc2 == 1:
246 z1 /= z2
247 return z1, z1[:1] * 0
248 elif lc1 < lc2:
249 return z1[:1] * 0, z1
250 else:
251 dlen = lc1 - lc2
252 scl = z2[0]
253 z2 /= scl
254 quo = np.empty(dlen + 1, dtype=z1.dtype)
255 i = 0
256 j = dlen
257 while i < j:
258 r = z1[i]
259 quo[i] = z1[i]
260 quo[dlen - i] = r
261 tmp = r * z2
262 z1[i:i + lc2] -= tmp
263 z1[j:j + lc2] -= tmp
264 i += 1
265 j -= 1
266 r = z1[i]
267 quo[i] = r
268 tmp = r * z2
269 z1[i:i + lc2] -= tmp
270 quo /= scl
271 rem = z1[i + 1:i - 1 + lc2].copy()
272 return quo, rem
273
274
275def _zseries_der(zs):
276 """Differentiate a z-series.
277
278 The derivative is with respect to x, not z. This is achieved using the
279 chain rule and the value of dx/dz given in the module notes.
280
281 Parameters
282 ----------
283 zs : z-series
284 The z-series to differentiate.
285
286 Returns
287 -------
288 derivative : z-series
289 The derivative
290
291 Notes
292 -----
293 The zseries for x (ns) has been multiplied by two in order to avoid
294 using floats that are incompatible with Decimal and likely other
295 specialized scalar types. This scaling has been compensated by
296 multiplying the value of zs by two also so that the two cancels in the
297 division.
298
299 """
300 n = len(zs) // 2
301 ns = np.array([-1, 0, 1], dtype=zs.dtype)
302 zs *= np.arange(-n, n + 1) * 2
303 d, r = _zseries_div(zs, ns)
304 return d
305
306
307def _zseries_int(zs):
308 """Integrate a z-series.
309
310 The integral is with respect to x, not z. This is achieved by a change
311 of variable using dx/dz given in the module notes.
312
313 Parameters
314 ----------
315 zs : z-series
316 The z-series to integrate
317
318 Returns
319 -------
320 integral : z-series
321 The indefinite integral
322
323 Notes
324 -----
325 The zseries for x (ns) has been multiplied by two in order to avoid
326 using floats that are incompatible with Decimal and likely other
327 specialized scalar types. This scaling has been compensated by
328 dividing the resulting zs by two.
329
330 """
331 n = 1 + len(zs) // 2
332 ns = np.array([-1, 0, 1], dtype=zs.dtype)
333 zs = _zseries_mul(zs, ns)
334 div = np.arange(-n, n + 1) * 2
335 zs[:n] /= div[:n]
336 zs[n + 1:] /= div[n + 1:]
337 zs[n] = 0
338 return zs
339
340#
341# Chebyshev series functions
342#
343
344
345def poly2cheb(pol):
346 """
347 Convert a polynomial to a Chebyshev series.
348
349 Convert an array representing the coefficients of a polynomial (relative
350 to the "standard" basis) ordered from lowest degree to highest, to an
351 array of the coefficients of the equivalent Chebyshev series, ordered
352 from lowest to highest degree.
353
354 Parameters
355 ----------
356 pol : array_like
357 1-D array containing the polynomial coefficients
358
359 Returns
360 -------
361 c : ndarray
362 1-D array containing the coefficients of the equivalent Chebyshev
363 series.
364
365 See Also
366 --------
367 cheb2poly
368
369 Notes
370 -----
371 The easy way to do conversions between polynomial basis sets
372 is to use the convert method of a class instance.
373
374 Examples
375 --------
376 >>> from numpy import polynomial as P
377 >>> p = P.Polynomial(range(4))
378 >>> p
379 Polynomial([0., 1., 2., 3.], domain=[-1., 1.], window=[-1., 1.], symbol='x')
380 >>> c = p.convert(kind=P.Chebyshev)
381 >>> c
382 Chebyshev([1. , 3.25, 1. , 0.75], domain=[-1., 1.], window=[-1., ...
383 >>> P.chebyshev.poly2cheb(range(4))
384 array([1. , 3.25, 1. , 0.75])
385
386 """
387 [pol] = pu.as_series([pol])
388 deg = len(pol) - 1
389 res = 0
390 for i in range(deg, -1, -1):
391 res = chebadd(chebmulx(res), pol[i])
392 return res
393
394
395def cheb2poly(c):
396 """
397 Convert a Chebyshev series to a polynomial.
398
399 Convert an array representing the coefficients of a Chebyshev series,
400 ordered from lowest degree to highest, to an array of the coefficients
401 of the equivalent polynomial (relative to the "standard" basis) ordered
402 from lowest to highest degree.
403
404 Parameters
405 ----------
406 c : array_like
407 1-D array containing the Chebyshev series coefficients, ordered
408 from lowest order term to highest.
409
410 Returns
411 -------
412 pol : ndarray
413 1-D array containing the coefficients of the equivalent polynomial
414 (relative to the "standard" basis) ordered from lowest order term
415 to highest.
416
417 See Also
418 --------
419 poly2cheb
420
421 Notes
422 -----
423 The easy way to do conversions between polynomial basis sets
424 is to use the convert method of a class instance.
425
426 Examples
427 --------
428 >>> from numpy import polynomial as P
429 >>> c = P.Chebyshev(range(4))
430 >>> c
431 Chebyshev([0., 1., 2., 3.], domain=[-1., 1.], window=[-1., 1.], symbol='x')
432 >>> p = c.convert(kind=P.Polynomial)
433 >>> p
434 Polynomial([-2., -8., 4., 12.], domain=[-1., 1.], window=[-1., 1.], ...
435 >>> P.chebyshev.cheb2poly(range(4))
436 array([-2., -8., 4., 12.])
437
438 """
439 from .polynomial import polyadd, polymulx, polysub
440
441 [c] = pu.as_series([c])
442 n = len(c)
443 if n < 3:
444 return c
445 else:
446 c0 = c[-2]
447 c1 = c[-1]
448 # i is the current degree of c1
449 for i in range(n - 1, 1, -1):
450 tmp = c0
451 c0 = polysub(c[i - 2], c1)
452 c1 = polyadd(tmp, polymulx(c1) * 2)
453 return polyadd(c0, polymulx(c1))
454
455
456#
457# These are constant arrays are of integer type so as to be compatible
458# with the widest range of other types, such as Decimal.
459#
460
461# Chebyshev default domain.
462chebdomain = np.array([-1., 1.])
463
464# Chebyshev coefficients representing zero.
465chebzero = np.array([0])
466
467# Chebyshev coefficients representing one.
468chebone = np.array([1])
469
470# Chebyshev coefficients representing the identity x.
471chebx = np.array([0, 1])
472
473
474def chebline(off, scl):
475 """
476 Chebyshev series whose graph is a straight line.
477
478 Parameters
479 ----------
480 off, scl : scalars
481 The specified line is given by ``off + scl*x``.
482
483 Returns
484 -------
485 y : ndarray
486 This module's representation of the Chebyshev series for
487 ``off + scl*x``.
488
489 See Also
490 --------
491 numpy.polynomial.polynomial.polyline
492 numpy.polynomial.legendre.legline
493 numpy.polynomial.laguerre.lagline
494 numpy.polynomial.hermite.hermline
495 numpy.polynomial.hermite_e.hermeline
496
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
504
505 """
506 if scl != 0:
507 return np.array([off, scl])
508 else:
509 return np.array([off])
510
511
512def chebfromroots(roots):
513 """
514 Generate a Chebyshev series with given roots.
515
516 The function returns the coefficients of the polynomial
517
518 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
519
520 in Chebyshev form, where the :math:`r_n` are the roots specified in
521 `roots`. If a zero has multiplicity n, then it must appear in `roots`
522 n times. For instance, if 2 is a root of multiplicity three and 3 is a
523 root of multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3].
524 The roots can appear in any order.
525
526 If the returned coefficients are `c`, then
527
528 .. math:: p(x) = c_0 + c_1 * T_1(x) + ... + c_n * T_n(x)
529
530 The coefficient of the last term is not generally 1 for monic
531 polynomials in Chebyshev form.
532
533 Parameters
534 ----------
535 roots : array_like
536 Sequence containing the roots.
537
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).
545
546 See Also
547 --------
548 numpy.polynomial.polynomial.polyfromroots
549 numpy.polynomial.legendre.legfromroots
550 numpy.polynomial.laguerre.lagfromroots
551 numpy.polynomial.hermite.hermfromroots
552 numpy.polynomial.hermite_e.hermefromroots
553
554 Examples
555 --------
556 >>> import numpy.polynomial.chebyshev as C
557 >>> C.chebfromroots((-1,0,1)) # x^3 - x relative to the standard basis
558 array([ 0. , -0.25, 0. , 0.25])
559 >>> j = complex(0,1)
560 >>> C.chebfromroots((-j,j)) # x^2 + 1 relative to the standard basis
561 array([1.5+0.j, 0. +0.j, 0.5+0.j])
562
563 """
564 return pu._fromroots(chebline, chebmul, roots)
565
566
567def chebadd(c1, c2):
568 """
569 Add one Chebyshev series to another.
570
571 Returns the sum of two Chebyshev series `c1` + `c2`. The arguments
572 are sequences of coefficients ordered from lowest order term to
573 highest, i.e., [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
574
575 Parameters
576 ----------
577 c1, c2 : array_like
578 1-D arrays of Chebyshev series coefficients ordered from low to
579 high.
580
581 Returns
582 -------
583 out : ndarray
584 Array representing the Chebyshev series of their sum.
585
586 See Also
587 --------
588 chebsub, chebmulx, chebmul, chebdiv, chebpow
589
590 Notes
591 -----
592 Unlike multiplication, division, etc., the sum of two Chebyshev series
593 is a Chebyshev series (without having to "reproject" the result onto
594 the basis set) so addition, just like that of "standard" polynomials,
595 is simply "component-wise."
596
597 Examples
598 --------
599 >>> from numpy.polynomial import chebyshev as C
600 >>> c1 = (1,2,3)
601 >>> c2 = (3,2,1)
602 >>> C.chebadd(c1,c2)
603 array([4., 4., 4.])
604
605 """
606 return pu._add(c1, c2)
607
608
609def chebsub(c1, c2):
610 """
611 Subtract one Chebyshev series from another.
612
613 Returns the difference of two Chebyshev series `c1` - `c2`. The
614 sequences of coefficients are from lowest order term to highest, i.e.,
615 [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
616
617 Parameters
618 ----------
619 c1, c2 : array_like
620 1-D arrays of Chebyshev series coefficients ordered from low to
621 high.
622
623 Returns
624 -------
625 out : ndarray
626 Of Chebyshev series coefficients representing their difference.
627
628 See Also
629 --------
630 chebadd, chebmulx, chebmul, chebdiv, chebpow
631
632 Notes
633 -----
634 Unlike multiplication, division, etc., the difference of two Chebyshev
635 series is a Chebyshev series (without having to "reproject" the result
636 onto the basis set) so subtraction, just like that of "standard"
637 polynomials, is simply "component-wise."
638
639 Examples
640 --------
641 >>> from numpy.polynomial import chebyshev as C
642 >>> c1 = (1,2,3)
643 >>> c2 = (3,2,1)
644 >>> C.chebsub(c1,c2)
645 array([-2., 0., 2.])
646 >>> C.chebsub(c2,c1) # -C.chebsub(c1,c2)
647 array([ 2., 0., -2.])
648
649 """
650 return pu._sub(c1, c2)
651
652
653def chebmulx(c):
654 """Multiply a Chebyshev series by x.
655
656 Multiply the polynomial `c` by x, where x is the independent
657 variable.
658
659
660 Parameters
661 ----------
662 c : array_like
663 1-D array of Chebyshev series coefficients ordered from low to
664 high.
665
666 Returns
667 -------
668 out : ndarray
669 Array representing the result of the multiplication.
670
671 See Also
672 --------
673 chebadd, chebsub, chebmul, chebdiv, chebpow
674
675 Examples
676 --------
677 >>> from numpy.polynomial import chebyshev as C
678 >>> C.chebmulx([1,2,3])
679 array([1. , 2.5, 1. , 1.5])
680
681 """
682 # c is a trimmed copy
683 [c] = pu.as_series([c])
684 # The zero series needs special treatment
685 if len(c) == 1 and c[0] == 0:
686 return c
687
688 prd = np.empty(len(c) + 1, dtype=c.dtype)
689 prd[0] = c[0] * 0
690 prd[1] = c[0]
691 if len(c) > 1:
692 tmp = c[1:] / 2
693 prd[2:] = tmp
694 prd[0:-2] += tmp
695 return prd
696
697
698def chebmul(c1, c2):
699 """
700 Multiply one Chebyshev series by another.
701
702 Returns the product of two Chebyshev series `c1` * `c2`. The arguments
703 are sequences of coefficients, from lowest order "term" to highest,
704 e.g., [1,2,3] represents the series ``T_0 + 2*T_1 + 3*T_2``.
705
706 Parameters
707 ----------
708 c1, c2 : array_like
709 1-D arrays of Chebyshev series coefficients ordered from low to
710 high.
711
712 Returns
713 -------
714 out : ndarray
715 Of Chebyshev series coefficients representing their product.
716
717 See Also
718 --------
719 chebadd, chebsub, chebmulx, chebdiv, chebpow
720
721 Notes
722 -----
723 In general, the (polynomial) product of two C-series results in terms
724 that are not in the Chebyshev polynomial basis set. Thus, to express
725 the product as a C-series, it is typically necessary to "reproject"
726 the product onto said basis set, which typically produces
727 "unintuitive live" (but correct) results; see Examples section below.
728
729 Examples
730 --------
731 >>> from numpy.polynomial import chebyshev as C
732 >>> c1 = (1,2,3)
733 >>> c2 = (3,2,1)
734 >>> C.chebmul(c1,c2) # multiplication requires "reprojection"
735 array([ 6.5, 12. , 12. , 4. , 1.5])
736
737 """
738 # c1, c2 are trimmed copies
739 [c1, c2] = pu.as_series([c1, c2])
740 z1 = _cseries_to_zseries(c1)
741 z2 = _cseries_to_zseries(c2)
742 prd = _zseries_mul(z1, z2)
743 ret = _zseries_to_cseries(prd)
744 return pu.trimseq(ret)
745
746
747def chebdiv(c1, c2):
748 """
749 Divide one Chebyshev series by another.
750
751 Returns the quotient-with-remainder of two Chebyshev series
752 `c1` / `c2`. The arguments are sequences of coefficients from lowest
753 order "term" to highest, e.g., [1,2,3] represents the series
754 ``T_0 + 2*T_1 + 3*T_2``.
755
756 Parameters
757 ----------
758 c1, c2 : array_like
759 1-D arrays of Chebyshev series coefficients ordered from low to
760 high.
761
762 Returns
763 -------
764 [quo, rem] : ndarrays
765 Of Chebyshev series coefficients representing the quotient and
766 remainder.
767
768 See Also
769 --------
770 chebadd, chebsub, chebmulx, chebmul, chebpow
771
772 Notes
773 -----
774 In general, the (polynomial) division of one C-series by another
775 results in quotient and remainder terms that are not in the Chebyshev
776 polynomial basis set. Thus, to express these results as C-series, it
777 is typically necessary to "reproject" the results onto said basis
778 set, which typically produces "unintuitive" (but correct) results;
779 see Examples section below.
780
781 Examples
782 --------
783 >>> from numpy.polynomial import chebyshev as C
784 >>> c1 = (1,2,3)
785 >>> c2 = (3,2,1)
786 >>> C.chebdiv(c1,c2) # quotient "intuitive," remainder not
787 (array([3.]), array([-8., -4.]))
788 >>> c2 = (0,1,2,3)
789 >>> C.chebdiv(c2,c1) # neither "intuitive"
790 (array([0., 2.]), array([-2., -4.]))
791
792 """
793 # c1, c2 are trimmed copies
794 [c1, c2] = pu.as_series([c1, c2])
795 if c2[-1] == 0:
796 raise ZeroDivisionError # FIXME: add message with details to exception
797
798 # note: this is more efficient than `pu._div(chebmul, c1, c2)`
799 lc1 = len(c1)
800 lc2 = len(c2)
801 if lc1 < lc2:
802 return c1[:1] * 0, c1
803 elif lc2 == 1:
804 return c1 / c2[-1], c1[:1] * 0
805 else:
806 z1 = _cseries_to_zseries(c1)
807 z2 = _cseries_to_zseries(c2)
808 quo, rem = _zseries_div(z1, z2)
809 quo = pu.trimseq(_zseries_to_cseries(quo))
810 rem = pu.trimseq(_zseries_to_cseries(rem))
811 return quo, rem
812
813
814def chebpow(c, pow, maxpower=16):
815 """Raise a Chebyshev series to a power.
816
817 Returns the Chebyshev series `c` raised to the power `pow`. The
818 argument `c` is a sequence of coefficients ordered from low to high.
819 i.e., [1,2,3] is the series ``T_0 + 2*T_1 + 3*T_2.``
820
821 Parameters
822 ----------
823 c : array_like
824 1-D array of Chebyshev series coefficients ordered from low to
825 high.
826 pow : integer
827 Power to which the series will be raised
828 maxpower : integer, optional
829 Maximum power allowed. This is mainly to limit growth of the series
830 to unmanageable size. Default is 16
831
832 Returns
833 -------
834 coef : ndarray
835 Chebyshev series of power.
836
837 See Also
838 --------
839 chebadd, chebsub, chebmulx, chebmul, chebdiv
840
841 Examples
842 --------
843 >>> from numpy.polynomial import chebyshev as C
844 >>> C.chebpow([1, 2, 3, 4], 2)
845 array([15.5, 22. , 16. , ..., 12.5, 12. , 8. ])
846
847 """
848 # note: this is more efficient than `pu._pow(chebmul, c1, c2)`, as it
849 # avoids converting between z and c series repeatedly
850
851 # c is a trimmed copy
852 [c] = pu.as_series([c])
853 power = int(pow)
854 if power != pow or power < 0:
855 raise ValueError("Power must be a non-negative integer.")
856 elif maxpower is not None and power > maxpower:
857 raise ValueError("Power is too large")
858 elif power == 0:
859 return np.array([1], dtype=c.dtype)
860 elif power == 1:
861 return c
862 else:
863 # This can be made more efficient by using powers of two
864 # in the usual way.
865 zs = _cseries_to_zseries(c)
866 prd = zs
867 for i in range(2, power + 1):
868 prd = np.convolve(prd, zs)
869 return _zseries_to_cseries(prd)
870
871
872def chebder(c, m=1, scl=1, axis=0):
873 """
874 Differentiate a Chebyshev series.
875
876 Returns the Chebyshev series coefficients `c` differentiated `m` times
877 along `axis`. At each iteration the result is multiplied by `scl` (the
878 scaling factor is for use in a linear change of variable). The argument
879 `c` is an array of coefficients from low to high degree along each
880 axis, e.g., [1,2,3] represents the series ``1*T_0 + 2*T_1 + 3*T_2``
881 while [[1,2],[1,2]] represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) +
882 2*T_0(x)*T_1(y) + 2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is
883 ``y``.
884
885 Parameters
886 ----------
887 c : array_like
888 Array of Chebyshev series coefficients. If c is multidimensional
889 the different axis correspond to different variables with the
890 degree in each axis given by the corresponding index.
891 m : int, optional
892 Number of derivatives taken, must be non-negative. (Default: 1)
893 scl : scalar, optional
894 Each differentiation is multiplied by `scl`. The end result is
895 multiplication by ``scl**m``. This is for use in a linear change of
896 variable. (Default: 1)
897 axis : int, optional
898 Axis over which the derivative is taken. (Default: 0).
899
900 Returns
901 -------
902 der : ndarray
903 Chebyshev series of the derivative.
904
905 See Also
906 --------
907 chebint
908
909 Notes
910 -----
911 In general, the result of differentiating a C-series needs to be
912 "reprojected" onto the C-series basis set. Thus, typically, the
913 result of this function is "unintuitive," albeit correct; see Examples
914 section below.
915
916 Examples
917 --------
918 >>> from numpy.polynomial import chebyshev as C
919 >>> c = (1,2,3,4)
920 >>> C.chebder(c)
921 array([14., 12., 24.])
922 >>> C.chebder(c,3)
923 array([96.])
924 >>> C.chebder(c,scl=-1)
925 array([-14., -12., -24.])
926 >>> C.chebder(c,2,-1)
927 array([12., 96.])
928
929 """
930 c = np.array(c, ndmin=1, copy=True)
931 if c.dtype.char in '?bBhHiIlLqQpP':
932 c = c.astype(np.double)
933 cnt = pu._as_int(m, "the order of derivation")
934 iaxis = pu._as_int(axis, "the axis")
935 if cnt < 0:
936 raise ValueError("The order of derivation must be non-negative")
937 iaxis = np.lib.array_utils.normalize_axis_index(iaxis, c.ndim)
938
939 if cnt == 0:
940 return c
941
942 c = np.moveaxis(c, iaxis, 0)
943 n = len(c)
944 if cnt >= n:
945 c = c[:1] * 0
946 else:
947 for i in range(cnt):
948 n = n - 1
949 c *= scl
950 der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
951 for j in range(n, 2, -1):
952 der[j - 1] = (2 * j) * c[j]
953 c[j - 2] += (j * c[j]) / (j - 2)
954 if n > 1:
955 der[1] = 4 * c[2]
956 der[0] = c[1]
957 c = der
958 c = np.moveaxis(c, 0, iaxis)
959 return c
960
961
962def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
963 """
964 Integrate a Chebyshev series.
965
966 Returns the Chebyshev series coefficients `c` integrated `m` times from
967 `lbnd` along `axis`. At each iteration the resulting series is
968 **multiplied** by `scl` and an integration constant, `k`, is added.
969 The scaling factor is for use in a linear change of variable. ("Buyer
970 beware": note that, depending on what one is doing, one may want `scl`
971 to be the reciprocal of what one might expect; for more information,
972 see the Notes section below.) The argument `c` is an array of
973 coefficients from low to high degree along each axis, e.g., [1,2,3]
974 represents the series ``T_0 + 2*T_1 + 3*T_2`` while [[1,2],[1,2]]
975 represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) + 2*T_0(x)*T_1(y) +
976 2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
977
978 Parameters
979 ----------
980 c : array_like
981 Array of Chebyshev series coefficients. If c is multidimensional
982 the different axis correspond to different variables with the
983 degree in each axis given by the corresponding index.
984 m : int, optional
985 Order of integration, must be positive. (Default: 1)
986 k : {[], list, scalar}, optional
987 Integration constant(s). The value of the first integral at zero
988 is the first value in the list, the value of the second integral
989 at zero is the second value, etc. If ``k == []`` (the default),
990 all constants are set to zero. If ``m == 1``, a single scalar can
991 be given instead of a list.
992 lbnd : scalar, optional
993 The lower bound of the integral. (Default: 0)
994 scl : scalar, optional
995 Following each integration the result is *multiplied* by `scl`
996 before the integration constant is added. (Default: 1)
997 axis : int, optional
998 Axis over which the integral is taken. (Default: 0).
999
1000 Returns
1001 -------
1002 S : ndarray
1003 C-series coefficients of the integral.
1004
1005 Raises
1006 ------
1007 ValueError
1008 If ``m < 1``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
1009 ``np.ndim(scl) != 0``.
1010
1011 See Also
1012 --------
1013 chebder
1014
1015 Notes
1016 -----
1017 Note that the result of each integration is *multiplied* by `scl`.
1018 Why is this important to note? Say one is making a linear change of
1019 variable :math:`u = ax + b` in an integral relative to `x`. Then
1020 :math:`dx = du/a`, so one will need to set `scl` equal to
1021 :math:`1/a`- perhaps not what one would have first thought.
1022
1023 Also note that, in general, the result of integrating a C-series needs
1024 to be "reprojected" onto the C-series basis set. Thus, typically,
1025 the result of this function is "unintuitive," albeit correct; see
1026 Examples section below.
1027
1028 Examples
1029 --------
1030 >>> from numpy.polynomial import chebyshev as C
1031 >>> c = (1,2,3)
1032 >>> C.chebint(c)
1033 array([ 0.5, -0.5, 0.5, 0.5])
1034 >>> C.chebint(c,3)
1035 array([ 0.03125 , -0.1875 , 0.04166667, -0.05208333, 0.01041667, # may vary
1036 0.00625 ])
1037 >>> C.chebint(c, k=3)
1038 array([ 3.5, -0.5, 0.5, 0.5])
1039 >>> C.chebint(c,lbnd=-2)
1040 array([ 8.5, -0.5, 0.5, 0.5])
1041 >>> C.chebint(c,scl=-2)
1042 array([-1., 1., -1., -1.])
1043
1044 """
1045 c = np.array(c, ndmin=1, copy=True)
1046 if c.dtype.char in '?bBhHiIlLqQpP':
1047 c = c.astype(np.double)
1048 if not np.iterable(k):
1049 k = [k]
1050 cnt = pu._as_int(m, "the order of integration")
1051 iaxis = pu._as_int(axis, "the axis")
1052 if cnt < 0:
1053 raise ValueError("The order of integration must be non-negative")
1054 if len(k) > cnt:
1055 raise ValueError("Too many integration constants")
1056 if np.ndim(lbnd) != 0:
1057 raise ValueError("lbnd must be a scalar.")
1058 if np.ndim(scl) != 0:
1059 raise ValueError("scl must be a scalar.")
1060 iaxis = np.lib.array_utils.normalize_axis_index(iaxis, c.ndim)
1061
1062 if cnt == 0:
1063 return c
1064
1065 c = np.moveaxis(c, iaxis, 0)
1066 k = list(k) + [0] * (cnt - len(k))
1067 for i in range(cnt):
1068 n = len(c)
1069 c *= scl
1070 if n == 1 and np.all(c[0] == 0):
1071 c[0] += k[i]
1072 else:
1073 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
1074 tmp[0] = c[0] * 0
1075 tmp[1] = c[0]
1076 if n > 1:
1077 tmp[2] = c[1] / 4
1078 for j in range(2, n):
1079 tmp[j + 1] = c[j] / (2 * (j + 1))
1080 tmp[j - 1] -= c[j] / (2 * (j - 1))
1081 tmp[0] += k[i] - chebval(lbnd, tmp)
1082 c = tmp
1083 c = np.moveaxis(c, 0, iaxis)
1084 return c
1085
1086
1087def chebval(x, c, tensor=True):
1088 """
1089 Evaluate a Chebyshev series at points x.
1090
1091 If `c` is of length `n + 1`, this function returns the value:
1092
1093 .. math:: p(x) = c_0 * T_0(x) + c_1 * T_1(x) + ... + c_n * T_n(x)
1094
1095 The parameter `x` is converted to an array only if it is a tuple or a
1096 list, otherwise it is treated as a scalar. In either case, either `x`
1097 or its elements must support multiplication and addition both with
1098 themselves and with the elements of `c`.
1099
1100 If `c` is a 1-D array, then ``p(x)`` will have the same shape as `x`. If
1101 `c` is multidimensional, then the shape of the result depends on the
1102 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
1103 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
1104 scalars have shape (,).
1105
1106 Trailing zeros in the coefficients will be used in the evaluation, so
1107 they should be avoided if efficiency is a concern.
1108
1109 Parameters
1110 ----------
1111 x : array_like, compatible object
1112 If `x` is a list or tuple, it is converted to an ndarray, otherwise
1113 it is left unchanged and treated as a scalar. In either case, `x`
1114 or its elements must support addition and multiplication with
1115 themselves and with the elements of `c`.
1116 c : array_like
1117 Array of coefficients ordered so that the coefficients for terms of
1118 degree n are contained in c[n]. If `c` is multidimensional the
1119 remaining indices enumerate multiple polynomials. In the two
1120 dimensional case the coefficients may be thought of as stored in
1121 the columns of `c`.
1122 tensor : boolean, optional
1123 If True, the shape of the coefficient array is extended with ones
1124 on the right, one for each dimension of `x`. Scalars have dimension 0
1125 for this action. The result is that every column of coefficients in
1126 `c` is evaluated for every element of `x`. If False, `x` is broadcast
1127 over the columns of `c` for the evaluation. This keyword is useful
1128 when `c` is multidimensional. The default value is True.
1129
1130 Returns
1131 -------
1132 values : ndarray, algebra_like
1133 The shape of the return value is described above.
1134
1135 See Also
1136 --------
1137 chebval2d, chebgrid2d, chebval3d, chebgrid3d
1138
1139 Notes
1140 -----
1141 The evaluation uses Clenshaw recursion, aka synthetic division.
1142
1143 """
1144 c = np.array(c, ndmin=1, copy=True)
1145 if c.dtype.char in '?bBhHiIlLqQpP':
1146 c = c.astype(np.double)
1147 if isinstance(x, (tuple, list)):
1148 x = np.asarray(x)
1149 if isinstance(x, np.ndarray) and tensor:
1150 c = c.reshape(c.shape + (1,) * x.ndim)
1151
1152 if len(c) == 1:
1153 c0 = c[0]
1154 c1 = 0
1155 elif len(c) == 2:
1156 c0 = c[0]
1157 c1 = c[1]
1158 else:
1159 x2 = 2 * x
1160 c0 = c[-2]
1161 c1 = c[-1]
1162 for i in range(3, len(c) + 1):
1163 tmp = c0
1164 c0 = c[-i] - c1
1165 c1 = tmp + c1 * x2
1166 return c0 + c1 * x
1167
1168
1169def chebval2d(x, y, c):
1170 """
1171 Evaluate a 2-D Chebyshev series at points (x, y).
1172
1173 This function returns the values:
1174
1175 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * T_i(x) * T_j(y)
1176
1177 The parameters `x` and `y` are converted to arrays only if they are
1178 tuples or a lists, otherwise they are treated as a scalars and they
1179 must have the same shape after conversion. In either case, either `x`
1180 and `y` or their elements must support multiplication and addition both
1181 with themselves and with the elements of `c`.
1182
1183 If `c` is a 1-D array a one is implicitly appended to its shape to make
1184 it 2-D. The shape of the result will be c.shape[2:] + x.shape.
1185
1186 Parameters
1187 ----------
1188 x, y : array_like, compatible objects
1189 The two dimensional series is evaluated at the points ``(x, y)``,
1190 where `x` and `y` must have the same shape. If `x` or `y` is a list
1191 or tuple, it is first converted to an ndarray, otherwise it is left
1192 unchanged and if it isn't an ndarray it is treated as a scalar.
1193 c : array_like
1194 Array of coefficients ordered so that the coefficient of the term
1195 of multi-degree i,j is contained in ``c[i,j]``. If `c` has
1196 dimension greater than 2 the remaining indices enumerate multiple
1197 sets of coefficients.
1198
1199 Returns
1200 -------