12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556 |
- """
- Objects for dealing with polynomials.
- This module provides a number of objects (mostly functions) useful for
- dealing with polynomials, including a `Polynomial` class that
- encapsulates the usual arithmetic operations. (General information
- on how this module represents and works with polynomial objects is in
- the docstring for its "parent" sub-package, `numpy.polynomial`).
- Constants
- ---------
- - `polydomain` -- Polynomial default domain, [-1,1].
- - `polyzero` -- (Coefficients of the) "zero polynomial."
- - `polyone` -- (Coefficients of the) constant polynomial 1.
- - `polyx` -- (Coefficients of the) identity map polynomial, ``f(x) = x``.
- Arithmetic
- ----------
- - `polyadd` -- add two polynomials.
- - `polysub` -- subtract one polynomial from another.
- - `polymul` -- multiply two polynomials.
- - `polydiv` -- divide one polynomial by another.
- - `polypow` -- raise a polynomial to an positive integer power
- - `polyval` -- evaluate a polynomial at given points.
- - `polyval2d` -- evaluate a 2D polynomial at given points.
- - `polyval3d` -- evaluate a 3D polynomial at given points.
- - `polygrid2d` -- evaluate a 2D polynomial on a Cartesian product.
- - `polygrid3d` -- evaluate a 3D polynomial on a Cartesian product.
- Calculus
- --------
- - `polyder` -- differentiate a polynomial.
- - `polyint` -- integrate a polynomial.
- Misc Functions
- --------------
- - `polyfromroots` -- create a polynomial with specified roots.
- - `polyroots` -- find the roots of a polynomial.
- - `polyvander` -- Vandermonde-like matrix for powers.
- - `polyvander2d` -- Vandermonde-like matrix for 2D power series.
- - `polyvander3d` -- Vandermonde-like matrix for 3D power series.
- - `polycompanion` -- companion matrix in power series form.
- - `polyfit` -- least-squares fit returning a polynomial.
- - `polytrim` -- trim leading coefficients from a polynomial.
- - `polyline` -- polynomial representing given straight line.
- Classes
- -------
- - `Polynomial` -- polynomial class.
- See Also
- --------
- `numpy.polynomial`
- """
- from __future__ import division, absolute_import, print_function
- __all__ = [
- 'polyzero', 'polyone', 'polyx', 'polydomain', 'polyline', 'polyadd',
- 'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow', 'polyval',
- 'polyder', 'polyint', 'polyfromroots', 'polyvander', 'polyfit',
- 'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d',
- 'polygrid2d', 'polygrid3d', 'polyvander2d', 'polyvander3d']
- import warnings
- import numpy as np
- import numpy.linalg as la
- from . import polyutils as pu
- from ._polybase import ABCPolyBase
- polytrim = pu.trimcoef
- #
- # These are constant arrays are of integer type so as to be compatible
- # with the widest range of other types, such as Decimal.
- #
- # Polynomial default domain.
- polydomain = np.array([-1, 1])
- # Polynomial coefficients representing zero.
- polyzero = np.array([0])
- # Polynomial coefficients representing one.
- polyone = np.array([1])
- # Polynomial coefficients representing the identity x.
- polyx = np.array([0, 1])
- #
- # Polynomial series functions
- #
- def polyline(off, scl):
- """
- Returns an array representing a linear polynomial.
- Parameters
- ----------
- off, scl : scalars
- The "y-intercept" and "slope" of the line, respectively.
- Returns
- -------
- y : ndarray
- This module's representation of the linear polynomial ``off +
- scl*x``.
- See Also
- --------
- chebline
- Examples
- --------
- >>> from numpy.polynomial import polynomial as P
- >>> P.polyline(1,-1)
- array([ 1, -1])
- >>> P.polyval(1, P.polyline(1,-1)) # should be 0
- 0.0
- """
- if scl != 0:
- return np.array([off, scl])
- else:
- return np.array([off])
- def polyfromroots(roots):
- """
- Generate a monic polynomial with given roots.
- Return the coefficients of the polynomial
- .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
- where the `r_n` are the roots specified in `roots`. If a zero has
- multiplicity n, then it must appear in `roots` n times. For instance,
- if 2 is a root of multiplicity three and 3 is a root of multiplicity 2,
- then `roots` looks something like [2, 2, 2, 3, 3]. The roots can appear
- in any order.
- If the returned coefficients are `c`, then
- .. math:: p(x) = c_0 + c_1 * x + ... + x^n
- The coefficient of the last term is 1 for monic polynomials in this
- form.
- Parameters
- ----------
- roots : array_like
- Sequence containing the roots.
- Returns
- -------
- out : ndarray
- 1-D array of the polynomial's coefficients If all the roots are
- real, then `out` is also real, otherwise it is complex. (see
- Examples below).
- See Also
- --------
- chebfromroots, legfromroots, lagfromroots, hermfromroots
- hermefromroots
- Notes
- -----
- The coefficients are determined by multiplying together linear factors
- of the form `(x - r_i)`, i.e.
- .. math:: p(x) = (x - r_0) (x - r_1) ... (x - r_n)
- where ``n == len(roots) - 1``; note that this implies that `1` is always
- returned for :math:`a_n`.
- Examples
- --------
- >>> from numpy.polynomial import polynomial as P
- >>> P.polyfromroots((-1,0,1)) # x(x - 1)(x + 1) = x^3 - x
- array([ 0., -1., 0., 1.])
- >>> j = complex(0,1)
- >>> P.polyfromroots((-j,j)) # complex returned, though values are real
- array([ 1.+0.j, 0.+0.j, 1.+0.j])
- """
- if len(roots) == 0:
- return np.ones(1)
- else:
- [roots] = pu.as_series([roots], trim=False)
- roots.sort()
- p = [polyline(-r, 1) for r in roots]
- n = len(p)
- while n > 1:
- m, r = divmod(n, 2)
- tmp = [polymul(p[i], p[i+m]) for i in range(m)]
- if r:
- tmp[0] = polymul(tmp[0], p[-1])
- p = tmp
- n = m
- return p[0]
- def polyadd(c1, c2):
- """
- Add one polynomial to another.
- Returns the sum of two polynomials `c1` + `c2`. The arguments are
- sequences of coefficients from lowest order term to highest, i.e.,
- [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``.
- Parameters
- ----------
- c1, c2 : array_like
- 1-D arrays of polynomial coefficients ordered from low to high.
- Returns
- -------
- out : ndarray
- The coefficient array representing their sum.
- See Also
- --------
- polysub, polymul, polydiv, polypow
- Examples
- --------
- >>> from numpy.polynomial import polynomial as P
- >>> c1 = (1,2,3)
- >>> c2 = (3,2,1)
- >>> sum = P.polyadd(c1,c2); sum
- array([ 4., 4., 4.])
- >>> P.polyval(2, sum) # 4 + 4(2) + 4(2**2)
- 28.0
- """
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] += c2
- ret = c1
- else:
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
- def polysub(c1, c2):
- """
- Subtract one polynomial from another.
- Returns the difference of two polynomials `c1` - `c2`. The arguments
- are sequences of coefficients from lowest order term to highest, i.e.,
- [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``.
- Parameters
- ----------
- c1, c2 : array_like
- 1-D arrays of polynomial coefficients ordered from low to
- high.
- Returns
- -------
- out : ndarray
- Of coefficients representing their difference.
- See Also
- --------
- polyadd, polymul, polydiv, polypow
- Examples
- --------
- >>> from numpy.polynomial import polynomial as P
- >>> c1 = (1,2,3)
- >>> c2 = (3,2,1)
- >>> P.polysub(c1,c2)
- array([-2., 0., 2.])
- >>> P.polysub(c2,c1) # -P.polysub(c1,c2)
- array([ 2., 0., -2.])
- """
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if len(c1) > len(c2):
- c1[:c2.size] -= c2
- ret = c1
- else:
- c2 = -c2
- c2[:c1.size] += c1
- ret = c2
- return pu.trimseq(ret)
- def polymulx(c):
- """Multiply a polynomial by x.
- Multiply the polynomial `c` by x, where x is the independent
- variable.
- Parameters
- ----------
- c : array_like
- 1-D array of polynomial coefficients ordered from low to
- high.
- Returns
- -------
- out : ndarray
- Array representing the result of the multiplication.
- Notes
- -----
- .. versionadded:: 1.5.0
- """
- # c is a trimmed copy
- [c] = pu.as_series([c])
- # The zero series needs special treatment
- if len(c) == 1 and c[0] == 0:
- return c
- prd = np.empty(len(c) + 1, dtype=c.dtype)
- prd[0] = c[0]*0
- prd[1:] = c
- return prd
- def polymul(c1, c2):
- """
- Multiply one polynomial by another.
- Returns the product of two polynomials `c1` * `c2`. The arguments are
- sequences of coefficients, from lowest order term to highest, e.g.,
- [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2.``
- Parameters
- ----------
- c1, c2 : array_like
- 1-D arrays of coefficients representing a polynomial, relative to the
- "standard" basis, and ordered from lowest order term to highest.
- Returns
- -------
- out : ndarray
- Of the coefficients of their product.
- See Also
- --------
- polyadd, polysub, polydiv, polypow
- Examples
- --------
- >>> from numpy.polynomial import polynomial as P
- >>> c1 = (1,2,3)
- >>> c2 = (3,2,1)
- >>> P.polymul(c1,c2)
- array([ 3., 8., 14., 8., 3.])
- """
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- ret = np.convolve(c1, c2)
- return pu.trimseq(ret)
- def polydiv(c1, c2):
- """
- Divide one polynomial by another.
- Returns the quotient-with-remainder of two polynomials `c1` / `c2`.
- The arguments are sequences of coefficients, from lowest order term
- to highest, e.g., [1,2,3] represents ``1 + 2*x + 3*x**2``.
- Parameters
- ----------
- c1, c2 : array_like
- 1-D arrays of polynomial coefficients ordered from low to high.
- Returns
- -------
- [quo, rem] : ndarrays
- Of coefficient series representing the quotient and remainder.
- See Also
- --------
- polyadd, polysub, polymul, polypow
- Examples
- --------
- >>> from numpy.polynomial import polynomial as P
- >>> c1 = (1,2,3)
- >>> c2 = (3,2,1)
- >>> P.polydiv(c1,c2)
- (array([ 3.]), array([-8., -4.]))
- >>> P.polydiv(c2,c1)
- (array([ 0.33333333]), array([ 2.66666667, 1.33333333]))
- """
- # c1, c2 are trimmed copies
- [c1, c2] = pu.as_series([c1, c2])
- if c2[-1] == 0:
- raise ZeroDivisionError()
- len1 = len(c1)
- len2 = len(c2)
- if len2 == 1:
- return c1/c2[-1], c1[:1]*0
- elif len1 < len2:
- return c1[:1]*0, c1
- else:
- dlen = len1 - len2
- scl = c2[-1]
- c2 = c2[:-1]/scl
- i = dlen
- j = len1 - 1
- while i >= 0:
- c1[i:j] -= c2*c1[j]
- i -= 1
- j -= 1
- return c1[j+1:]/scl, pu.trimseq(c1[:j+1])
- def polypow(c, pow, maxpower=None):
- """Raise a polynomial to a power.
- Returns the polynomial `c` raised to the power `pow`. The argument
- `c` is a sequence of coefficients ordered from low to high. i.e.,
- [1,2,3] is the series ``1 + 2*x + 3*x**2.``
- Parameters
- ----------
- c : array_like
- 1-D array of array of series coefficients ordered from low to
- high degree.
- pow : integer
- Power to which the series will be raised
- maxpower : integer, optional
- Maximum power allowed. This is mainly to limit growth of the series
- to unmanageable size. Default is 16
- Returns
- -------
- coef : ndarray
- Power series of power.
- See Also
- --------
- polyadd, polysub, polymul, polydiv
- Examples
- --------
- """
- # c is a trimmed copy
- [c] = pu.as_series([c])
- power = int(pow)
- if power != pow or power < 0:
- raise ValueError("Power must be a non-negative integer.")
- elif maxpower is not None and power > maxpower:
- raise ValueError("Power is too large")
- elif power == 0:
- return np.array([1], dtype=c.dtype)
- elif power == 1:
- return c
- else:
- # This can be made more efficient by using powers of two
- # in the usual way.
- prd = c
- for i in range(2, power + 1):
- prd = np.convolve(prd, c)
- return prd
- def polyder(c, m=1, scl=1, axis=0):
- """
- Differentiate a polynomial.
- Returns the polynomial coefficients `c` differentiated `m` times along
- `axis`. At each iteration the result is multiplied by `scl` (the
- scaling factor is for use in a linear change of variable). The
- argument `c` is an array of coefficients from low to high degree along
- each axis, e.g., [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``
- while [[1,2],[1,2]] represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is
- ``x`` and axis=1 is ``y``.
- Parameters
- ----------
- c : array_like
- Array of polynomial coefficients. If c is multidimensional the
- different axis correspond to different variables with the degree
- in each axis given by the corresponding index.
- m : int, optional
- Number of derivatives taken, must be non-negative. (Default: 1)
- scl : scalar, optional
- Each differentiation is multiplied by `scl`. The end result is
- multiplication by ``scl**m``. This is for use in a linear change
- of variable. (Default: 1)
- axis : int, optional
- Axis over which the derivative is taken. (Default: 0).
- .. versionadded:: 1.7.0
- Returns
- -------
- der : ndarray
- Polynomial coefficients of the derivative.
- See Also
- --------
- polyint
- Examples
- --------
- >>> from numpy.polynomial import polynomial as P
- >>> c = (1,2,3,4) # 1 + 2x + 3x**2 + 4x**3
- >>> P.polyder(c) # (d/dx)(c) = 2 + 6x + 12x**2
- array([ 2., 6., 12.])
- >>> P.polyder(c,3) # (d**3/dx**3)(c) = 24
- array([ 24.])
- >>> P.polyder(c,scl=-1) # (d/d(-x))(c) = -2 - 6x - 12x**2
- array([ -2., -6., -12.])
- >>> P.polyder(c,2,-1) # (d**2/d(-x)**2)(c) = 6 + 24x
- array([ 6., 24.])
- """
- c = np.array(c, ndmin=1, copy=1)
- if c.dtype.char in '?bBhHiIlLqQpP':
- # astype fails with NA
- c = c + 0.0
- cdt = c.dtype
- cnt, iaxis = [int(t) for t in [m, axis]]
- if cnt != m:
- raise ValueError("The order of derivation must be integer")
- if cnt < 0:
- raise ValueError("The order of derivation must be non-negative")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
- if not -c.ndim <= iaxis < c.ndim:
- raise ValueError("The axis is out of range")
- if iaxis < 0:
- iaxis += c.ndim
- if cnt == 0:
- return c
- c = np.rollaxis(c, iaxis)
- n = len(c)
- if cnt >= n:
- c = c[:1]*0
- else:
- for i in range(cnt):
- n = n - 1
- c *= scl
- der = np.empty((n,) + c.shape[1:], dtype=cdt)
- for j in range(n, 0, -1):
- der[j - 1] = j*c[j]
- c = der
- c = np.rollaxis(c, 0, iaxis + 1)
- return c
- def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
- """
- Integrate a polynomial.
- Returns the polynomial coefficients `c` integrated `m` times from
- `lbnd` along `axis`. At each iteration the resulting series is
- **multiplied** by `scl` and an integration constant, `k`, is added.
- The scaling factor is for use in a linear change of variable. ("Buyer
- beware": note that, depending on what one is doing, one may want `scl`
- to be the reciprocal of what one might expect; for more information,
- see the Notes section below.) The argument `c` is an array of
- coefficients, from low to high degree along each axis, e.g., [1,2,3]
- represents the polynomial ``1 + 2*x + 3*x**2`` while [[1,2],[1,2]]
- represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is ``x`` and axis=1 is
- ``y``.
- Parameters
- ----------
- c : array_like
- 1-D array of polynomial coefficients, ordered from low to high.
- m : int, optional
- Order of integration, must be positive. (Default: 1)
- k : {[], list, scalar}, optional
- Integration constant(s). The value of the first integral at zero
- is the first value in the list, the value of the second integral
- at zero is the second value, etc. If ``k == []`` (the default),
- all constants are set to zero. If ``m == 1``, a single scalar can
- be given instead of a list.
- lbnd : scalar, optional
- The lower bound of the integral. (Default: 0)
- scl : scalar, optional
- Following each integration the result is *multiplied* by `scl`
- before the integration constant is added. (Default: 1)
- axis : int, optional
- Axis over which the integral is taken. (Default: 0).
- .. versionadded:: 1.7.0
- Returns
- -------
- S : ndarray
- Coefficient array of the integral.
- Raises
- ------
- ValueError
- If ``m < 1``, ``len(k) > m``.
- See Also
- --------
- polyder
- Notes
- -----
- Note that the result of each integration is *multiplied* by `scl`. Why
- is this important to note? Say one is making a linear change of
- variable :math:`u = ax + b` in an integral relative to `x`. Then
- .. math::`dx = du/a`, so one will need to set `scl` equal to
- :math:`1/a` - perhaps not what one would have first thought.
- Examples
- --------
- >>> from numpy.polynomial import polynomial as P
- >>> c = (1,2,3)
- >>> P.polyint(c) # should return array([0, 1, 1, 1])
- array([ 0., 1., 1., 1.])
- >>> P.polyint(c,3) # should return array([0, 0, 0, 1/6, 1/12, 1/20])
- array([ 0. , 0. , 0. , 0.16666667, 0.08333333,
- 0.05 ])
- >>> P.polyint(c,k=3) # should return array([3, 1, 1, 1])
- array([ 3., 1., 1., 1.])
- >>> P.polyint(c,lbnd=-2) # should return array([6, 1, 1, 1])
- array([ 6., 1., 1., 1.])
- >>> P.polyint(c,scl=-2) # should return array([0, -2, -2, -2])
- array([ 0., -2., -2., -2.])
- """
- c = np.array(c, ndmin=1, copy=1)
- if c.dtype.char in '?bBhHiIlLqQpP':
- # astype doesn't preserve mask attribute.
- c = c + 0.0
- cdt = c.dtype
- if not np.iterable(k):
- k = [k]
- cnt, iaxis = [int(t) for t in [m, axis]]
- if cnt != m:
- raise ValueError("The order of integration must be integer")
- if cnt < 0:
- raise ValueError("The order of integration must be non-negative")
- if len(k) > cnt:
- raise ValueError("Too many integration constants")
- if iaxis != axis:
- raise ValueError("The axis must be integer")
- if not -c.ndim <= iaxis < c.ndim:
- raise ValueError("The axis is out of range")
- if iaxis < 0:
- iaxis += c.ndim
- if cnt == 0:
- return c
- k = list(k) + [0]*(cnt - len(k))
- c = np.rollaxis(c, iaxis)
- for i in range(cnt):
- n = len(c)
- c *= scl
- if n == 1 and np.all(c[0] == 0):
- c[0] += k[i]
- else:
- tmp = np.empty((n + 1,) + c.shape[1:], dtype=cdt)
- tmp[0] = c[0]*0
- tmp[1] = c[0]
- for j in range(1, n):
- tmp[j + 1] = c[j]/(j + 1)
- tmp[0] += k[i] - polyval(lbnd, tmp)
- c = tmp
- c = np.rollaxis(c, 0, iaxis + 1)
- return c
- def polyval(x, c, tensor=True):
- """
- Evaluate a polynomial at points x.
- If `c` is of length `n + 1`, this function returns the value
- .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n
- The parameter `x` is converted to an array only if it is a tuple or a
- list, otherwise it is treated as a scalar. In either case, either `x`
- or its elements must support multiplication and addition both with
- themselves and with the elements of `c`.
- If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
- `c` is multidimensional, then the shape of the result depends on the
- value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
- x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
- scalars have shape (,).
- Trailing zeros in the coefficients will be used in the evaluation, so
- they should be avoided if efficiency is a concern.
- Parameters
- ----------
- x : array_like, compatible object
- If `x` is a list or tuple, it is converted to an ndarray, otherwise
- it is left unchanged and treated as a scalar. In either case, `x`
- or its elements must support addition and multiplication with
- with themselves and with the elements of `c`.
- c : array_like
- Array of coefficients ordered so that the coefficients for terms of
- degree n are contained in c[n]. If `c` is multidimensional the
- remaining indices enumerate multiple polynomials. In the two
- dimensional case the coefficients may be thought of as stored in
- the columns of `c`.
- tensor : boolean, optional
- If True, the shape of the coefficient array is extended with ones
- on the right, one for each dimension of `x`. Scalars have dimension 0
- for this action. The result is that every column of coefficients in
- `c` is evaluated for every element of `x`. If False, `x` is broadcast
- over the columns of `c` for the evaluation. This keyword is useful
- when `c` is multidimensional. The default value is True.
- .. versionadded:: 1.7.0
- Returns
- -------
- values : ndarray, compatible object
- The shape of the returned array is described above.
- See Also
- --------
- polyval2d, polygrid2d, polyval3d, polygrid3d
- Notes
- -----
- The evaluation uses Horner's method.
- Examples
- --------
- >>> from numpy.polynomial.polynomial import polyval
- >>> polyval(1, [1,2,3])
- 6.0
- >>> a = np.arange(4).reshape(2,2)
- >>> a
- array([[0, 1],
- [2, 3]])
- >>> polyval(a, [1,2,3])
- array([[ 1., 6.],
- [ 17., 34.]])
- >>> coef = np.arange(4).reshape(2,2) # multidimensional coefficients
- >>> coef
- array([[0, 1],
- [2, 3]])
- >>> polyval([1,2], coef, tensor=True)
- array([[ 2., 4.],
- [ 4., 7.]])
- >>> polyval([1,2], coef, tensor=False)
- array([ 2., 7.])
- """
- c = np.array(c, ndmin=1, copy=0)
- if c.dtype.char in '?bBhHiIlLqQpP':
- # astype fails with NA
- c = c + 0.0
- if isinstance(x, (tuple, list)):
- x = np.asarray(x)
- if isinstance(x, np.ndarray) and tensor:
- c = c.reshape(c.shape + (1,)*x.ndim)
- c0 = c[-1] + x*0
- for i in range(2, len(c) + 1):
- c0 = c[-i] + c0*x
- return c0
- def polyval2d(x, y, c):
- """
- Evaluate a 2-D polynomial at points (x, y).
- This function returns the value
- .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * x^i * y^j
- The parameters `x` and `y` are converted to arrays only if they are
- tuples or a lists, otherwise they are treated as a scalars and they
- must have the same shape after conversion. In either case, either `x`
- and `y` or their elements must support multiplication and addition both
- with themselves and with the elements of `c`.
- If `c` has fewer than two dimensions, ones are implicitly appended to
- its shape to make it 2-D. The shape of the result will be c.shape[2:] +
- x.shape.
- Parameters
- ----------
- x, y : array_like, compatible objects
- The two dimensional series is evaluated at the points `(x, y)`,
- where `x` and `y` must have the same shape. If `x` or `y` is a list
- or tuple, it is first converted to an ndarray, otherwise it is left
- unchanged and, if it isn't an ndarray, it is treated as a scalar.
- c : array_like
- Array of coefficients ordered so that the coefficient of the term
- of multi-degree i,j is contained in `c[i,j]`. If `c` has
- dimension greater than two the remaining indices enumerate multiple
- sets of coefficients.
- Returns
- -------
- values : ndarray, compatible object
- The values of the two dimensional polynomial at points formed with
- pairs of corresponding values from `x` and `y`.
- See Also
- --------
- polyval, polygrid2d, polyval3d, polygrid3d
- Notes
- -----
- .. versionadded:: 1.7.0
- """
- try:
- x, y = np.array((x, y), copy=0)
- except:
- raise ValueError('x, y are incompatible')
- c = polyval(x, c)
- c = polyval(y, c, tensor=False)
- return c
- def polygrid2d(x, y, c):
- """
- Evaluate a 2-D polynomial on the Cartesian product of x and y.
- This function returns the values:
- .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * a^i * b^j
- where the points `(a, b)` consist of all pairs formed by taking
- `a` from `x` and `b` from `y`. The resulting points form a grid with
- `x` in the first dimension and `y` in the second.
- The parameters `x` and `y` are converted to arrays only if they are
- tuples or a lists, otherwise they are treated as a scalars. In either
- case, either `x` and `y` or their elements must support multiplication
- and addition both with themselves and with the elements of `c`.
- If `c` has fewer than two dimensions, ones are implicitly appended to
- its shape to make it 2-D. The shape of the result will be c.shape[2:] +
- x.shape + y.shape.
- Parameters
- ----------
- x, y : array_like, compatible objects
- The two dimensional series is evaluated at the points in the
- Cartesian product of `x` and `y`. If `x` or `y` is a list or
- tuple, it is first converted to an ndarray, otherwise it is left
- unchanged and, if it isn't an ndarray, it is treated as a scalar.
- c : array_like
- Array of coefficients ordered so that the coefficients for terms of
- degree i,j are contained in ``c[i,j]``. If `c` has dimension
- greater than two the remaining indices enumerate multiple sets of
- coefficients.
- Returns
- -------
- values : ndarray, compatible object
- The values of the two dimensional polynomial at points in the Cartesian
- product of `x` and `y`.
- See Also
- --------
- polyval, polyval2d, polyval3d, polygrid3d
- Notes
- -----
- .. versionadded:: 1.7.0
- """
- c = polyval(x, c)
- c = polyval(y, c)
- return c
- def polyval3d(x, y, z, c):
- """
- Evaluate a 3-D polynomial at points (x, y, z).
- This function returns the values:
- .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * x^i * y^j * z^k
- The parameters `x`, `y`, and `z` are converted to arrays only if
- they are tuples or a lists, otherwise they are treated as a scalars and
- they must have the same shape after conversion. In either case, either
- `x`, `y`, and `z` or their elements must support multiplication and
- addition both with themselves and with the elements of `c`.
- If `c` has fewer than 3 dimensions, ones are implicitly appended to its
- shape to make it 3-D. The shape of the result will be c.shape[3:] +
- x.shape.
- Parameters
- ----------
- x, y, z : array_like, compatible object
- The three dimensional series is evaluated at the points
- `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
- any of `x`, `y`, or `z` is a list or tuple, it is first converted
- to an ndarray, otherwise it is left unchanged and if it isn't an
- ndarray it is treated as a scalar.
- c : array_like
- Array of coefficients ordered so that the coefficient of the term of
- multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
- greater than 3 the remaining indices enumerate multiple sets of
- coefficients.
- Returns
- -------
- values : ndarray, compatible object
- The values of the multidimensional polynomial on points formed with
- triples of corresponding values from `x`, `y`, and `z`.
- See Also
- --------
- polyval, polyval2d, polygrid2d, polygrid3d
- Notes
- -----
- .. versionadded:: 1.7.0
- """
- try:
- x, y, z = np.array((x, y, z), copy=0)
- except:
- raise ValueError('x, y, z are incompatible')
- c = polyval(x, c)
- c = polyval(y, c, tensor=False)
- c = polyval(z, c, tensor=False)
- return c
- def polygrid3d(x, y, z, c):
- """
- Evaluate a 3-D polynomial on the Cartesian product of x, y and z.
- This function returns the values:
- .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * a^i * b^j * c^k
- where the points `(a, b, c)` consist of all triples formed by taking
- `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
- a grid with `x` in the first dimension, `y` in the second, and `z` in
- the third.
- The parameters `x`, `y`, and `z` are converted to arrays only if they
- are tuples or a lists, otherwise they are treated as a scalars. In
- either case, either `x`, `y`, and `z` or their elements must support
- multiplication and addition both with themselves and with the elements
- of `c`.
- If `c` has fewer than three dimensions, ones are implicitly appended to
- its shape to make it 3-D. The shape of the result will be c.shape[3:] +
- x.shape + y.shape + z.shape.
- Parameters
- ----------
- x, y, z : array_like, compatible objects
- The three dimensional series is evaluated at the points in the
- Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
- list or tuple, it is first converted to an ndarray, otherwise it is
- left unchanged and, if it isn't an ndarray, it is treated as a
- scalar.
- c : array_like
- Array of coefficients ordered so that the coefficients for terms of
- degree i,j are contained in ``c[i,j]``. If `c` has dimension
- greater than two the remaining indices enumerate multiple sets of
- coefficients.
- Returns
- -------
- values : ndarray, compatible object
- The values of the two dimensional polynomial at points in the Cartesian
- product of `x` and `y`.
- See Also
- --------
- polyval, polyval2d, polygrid2d, polyval3d
- Notes
- -----
- .. versionadded:: 1.7.0
- """
- c = polyval(x, c)
- c = polyval(y, c)
- c = polyval(z, c)
- return c
- def polyvander(x, deg):
- """Vandermonde matrix of given degree.
- Returns the Vandermonde matrix of degree `deg` and sample points
- `x`. The Vandermonde matrix is defined by
- .. math:: V[..., i] = x^i,
- where `0 <= i <= deg`. The leading indices of `V` index the elements of
- `x` and the last index is the power of `x`.
- If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
- matrix ``V = polyvander(x, n)``, then ``np.dot(V, c)`` and
- ``polyval(x, c)`` are the same up to roundoff. This equivalence is
- useful both for least squares fitting and for the evaluation of a large
- number of polynomials of the same degree and sample points.
- Parameters
- ----------
- x : array_like
- Array of points. The dtype is converted to float64 or complex128
- depending on whether any of the elements are complex. If `x` is
- scalar it is converted to a 1-D array.
- deg : int
- Degree of the resulting matrix.
- Returns
- -------
- vander : ndarray.
- The Vandermonde matrix. The shape of the returned matrix is
- ``x.shape + (deg + 1,)``, where the last index is the power of `x`.
- The dtype will be the same as the converted `x`.
- See Also
- --------
- polyvander2d, polyvander3d
- """
- ideg = int(deg)
- if ideg != deg:
- raise ValueError("deg must be integer")
- if ideg < 0:
- raise ValueError("deg must be non-negative")
- x = np.array(x, copy=0, ndmin=1) + 0.0
- dims = (ideg + 1,) + x.shape
- dtyp = x.dtype
- v = np.empty(dims, dtype=dtyp)
- v[0] = x*0 + 1
- if ideg > 0:
- v[1] = x
- for i in range(2, ideg + 1):
- v[i] = v[i-1]*x
- return np.rollaxis(v, 0, v.ndim)
- def polyvander2d(x, y, deg):
- """Pseudo-Vandermonde matrix of given degrees.
- Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
- points `(x, y)`. The pseudo-Vandermonde matrix is defined by
- .. math:: V[..., deg[1]*i + j] = x^i * y^j,
- where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
- `V` index the points `(x, y)` and the last index encodes the powers of
- `x` and `y`.
- If ``V = polyvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
- correspond to the elements of a 2-D coefficient array `c` of shape
- (xdeg + 1, ydeg + 1) in the order
- .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
- and ``np.dot(V, c.flat)`` and ``polyval2d(x, y, c)`` will be the same
- up to roundoff. This equivalence is useful both for least squares
- fitting and for the evaluation of a large number of 2-D polynomials
- of the same degrees and sample points.
- Parameters
- ----------
- x, y : array_like
- Arrays of point coordinates, all of the same shape. The dtypes
- will be converted to either float64 or complex128 depending on
- whether any of the elements are complex. Scalars are converted to
- 1-D arrays.
- deg : list of ints
- List of maximum degrees of the form [x_deg, y_deg].
- Returns
- -------
- vander2d : ndarray
- The shape of the returned matrix is ``x.shape + (order,)``, where
- :math:`order = (deg[0]+1)*(deg([1]+1)`. The dtype will be the same
- as the converted `x` and `y`.
- See Also
- --------
- polyvander, polyvander3d. polyval2d, polyval3d
- """
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy = ideg
- x, y = np.array((x, y), copy=0) + 0.0
- vx = polyvander(x, degx)
- vy = polyvander(y, degy)
- v = vx[..., None]*vy[..., None,:]
- # einsum bug
- #v = np.einsum("...i,...j->...ij", vx, vy)
- return v.reshape(v.shape[:-2] + (-1,))
- def polyvander3d(x, y, z, deg):
- """Pseudo-Vandermonde matrix of given degrees.
- Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
- points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
- then The pseudo-Vandermonde matrix is defined by
- .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = x^i * y^j * z^k,
- where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
- indices of `V` index the points `(x, y, z)` and the last index encodes
- the powers of `x`, `y`, and `z`.
- If ``V = polyvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
- of `V` correspond to the elements of a 3-D coefficient array `c` of
- shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
- .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
- and ``np.dot(V, c.flat)`` and ``polyval3d(x, y, z, c)`` will be the
- same up to roundoff. This equivalence is useful both for least squares
- fitting and for the evaluation of a large number of 3-D polynomials
- of the same degrees and sample points.
- Parameters
- ----------
- x, y, z : array_like
- Arrays of point coordinates, all of the same shape. The dtypes will
- be converted to either float64 or complex128 depending on whether
- any of the elements are complex. Scalars are converted to 1-D
- arrays.
- deg : list of ints
- List of maximum degrees of the form [x_deg, y_deg, z_deg].
- Returns
- -------
- vander3d : ndarray
- The shape of the returned matrix is ``x.shape + (order,)``, where
- :math:`order = (deg[0]+1)*(deg([1]+1)*(deg[2]+1)`. The dtype will
- be the same as the converted `x`, `y`, and `z`.
- See Also
- --------
- polyvander, polyvander3d. polyval2d, polyval3d
- Notes
- -----
- .. versionadded:: 1.7.0
- """
- ideg = [int(d) for d in deg]
- is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
- if is_valid != [1, 1, 1]:
- raise ValueError("degrees must be non-negative integers")
- degx, degy, degz = ideg
- x, y, z = np.array((x, y, z), copy=0) + 0.0
- vx = polyvander(x, degx)
- vy = polyvander(y, degy)
- vz = polyvander(z, degz)
- v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
- # einsum bug
- #v = np.einsum("...i, ...j, ...k->...ijk", vx, vy, vz)
- return v.reshape(v.shape[:-3] + (-1,))
- def polyfit(x, y, deg, rcond=None, full=False, w=None):
- """
- Least-squares fit of a polynomial to data.
- Return the coefficients of a polynomial of degree `deg` that is the
- least squares fit to the data values `y` given at points `x`. If `y` is
- 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
- fits are done, one for each column of `y`, and the resulting
- coefficients are stored in the corresponding columns of a 2-D return.
- The fitted polynomial(s) are in the form
- .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n,
- where `n` is `deg`.
- Parameters
- ----------
- x : array_like, shape (`M`,)
- x-coordinates of the `M` sample (data) points ``(x[i], y[i])``.
- y : array_like, shape (`M`,) or (`M`, `K`)
- y-coordinates of the sample points. Several sets of sample points
- sharing the same x-coordinates can be (independently) fit with one
- call to `polyfit` by passing in for `y` a 2-D array that contains
- one data set per column.
- deg : int or 1-D array_like
- Degree(s) of the fitting polynomials. If `deg` is a single integer
- all terms up to and including the `deg`'th term are included in the
- fit. For Numpy versions >= 1.11 a list of integers specifying the
- degrees of the terms to include may be used instead.
- rcond : float, optional
- Relative condition number of the fit. Singular values smaller
- than `rcond`, relative to the largest singular value, will be
- ignored. The default value is ``len(x)*eps``, where `eps` is the
- relative precision of the platform's float type, about 2e-16 in
- most cases.
- full : bool, optional
- Switch determining the nature of the return value. When ``False``
- (the default) just the coefficients are returned; when ``True``,
- diagnostic information from the singular value decomposition (used
- to solve the fit's matrix equation) is also returned.
- w : array_like, shape (`M`,), optional
- Weights. If not None, the contribution of each point
- ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the
- weights are chosen so that the errors of the products ``w[i]*y[i]``
- all have the same variance. The default value is None.
- .. versionadded:: 1.5.0
- Returns
- -------
- coef : ndarray, shape (`deg` + 1,) or (`deg` + 1, `K`)
- Polynomial coefficients ordered from low to high. If `y` was 2-D,
- the coefficients in column `k` of `coef` represent the polynomial
- fit to the data in `y`'s `k`-th column.
- [residuals, rank, singular_values, rcond] : list
- These values are only returned if `full` = True
- resid -- sum of squared residuals of the least squares fit
- rank -- the numerical rank of the scaled Vandermonde matrix
- sv -- singular values of the scaled Vandermonde matrix
- rcond -- value of `rcond`.
- For more details, see `linalg.lstsq`.
- Raises
- ------
- RankWarning
- Raised if the matrix in the least-squares fit is rank deficient.
- The warning is only raised if `full` == False. The warnings can
- be turned off by:
- >>> import warnings
- >>> warnings.simplefilter('ignore', RankWarning)
- See Also
- --------
- chebfit, legfit, lagfit, hermfit, hermefit
- polyval : Evaluates a polynomial.
- polyvander : Vandermonde matrix for powers.
- linalg.lstsq : Computes a least-squares fit from the matrix.
- scipy.interpolate.UnivariateSpline : Computes spline fits.
- Notes
- -----
- The solution is the coefficients of the polynomial `p` that minimizes
- the sum of the weighted squared errors
- .. math :: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
- where the :math:`w_j` are the weights. This problem is solved by
- setting up the (typically) over-determined matrix equation:
- .. math :: V(x) * c = w * y,
- where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
- coefficients to be solved for, `w` are the weights, and `y` are the
- observed values. This equation is then solved using the singular value
- decomposition of `V`.
- If some of the singular values of `V` are so small that they are
- neglected (and `full` == ``False``), a `RankWarning` will be raised.
- This means that the coefficient values may be poorly determined.
- Fitting to a lower order polynomial will usually get rid of the warning
- (but may not be what you want, of course; if you have independent
- reason(s) for choosing the degree which isn't working, you may have to:
- a) reconsider those reasons, and/or b) reconsider the quality of your
- data). The `rcond` parameter can also be set to a value smaller than
- its default, but the resulting fit may be spurious and have large
- contributions from roundoff error.
- Polynomial fits using double precision tend to "fail" at about
- (polynomial) degree 20. Fits using Chebyshev or Legendre series are
- generally better conditioned, but much can still depend on the
- distribution of the sample points and the smoothness of the data. If
- the quality of the fit is inadequate, splines may be a good
- alternative.
- Examples
- --------
- >>> from numpy.polynomial import polynomial as P
- >>> x = np.linspace(-1,1,51) # x "data": [-1, -0.96, ..., 0.96, 1]
- >>> y = x**3 - x + np.random.randn(len(x)) # x^3 - x + N(0,1) "noise"
- >>> c, stats = P.polyfit(x,y,3,full=True)
- >>> c # c[0], c[2] should be approx. 0, c[1] approx. -1, c[3] approx. 1
- array([ 0.01909725, -1.30598256, -0.00577963, 1.02644286])
- >>> stats # note the large SSR, explaining the rather poor results
- [array([ 38.06116253]), 4, array([ 1.38446749, 1.32119158, 0.50443316,
- 0.28853036]), 1.1324274851176597e-014]
- Same thing without the added noise
- >>> y = x**3 - x
- >>> c, stats = P.polyfit(x,y,3,full=True)
- >>> c # c[0], c[2] should be "very close to 0", c[1] ~= -1, c[3] ~= 1
- array([ -1.73362882e-17, -1.00000000e+00, -2.67471909e-16,
- 1.00000000e+00])
- >>> stats # note the minuscule SSR
- [array([ 7.46346754e-31]), 4, array([ 1.38446749, 1.32119158,
- 0.50443316, 0.28853036]), 1.1324274851176597e-014]
- """
- x = np.asarray(x) + 0.0
- y = np.asarray(y) + 0.0
- deg = np.asarray(deg)
- # check arguments.
- if deg.ndim > 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
- raise TypeError("deg must be an int or non-empty 1-D array of int")
- if deg.min() < 0:
- raise ValueError("expected deg >= 0")
- if x.ndim != 1:
- raise TypeError("expected 1D vector for x")
- if x.size == 0:
- raise TypeError("expected non-empty vector for x")
- if y.ndim < 1 or y.ndim > 2:
- raise TypeError("expected 1D or 2D array for y")
- if len(x) != len(y):
- raise TypeError("expected x and y to have same length")
- if deg.ndim == 0:
- lmax = deg
- order = lmax + 1
- van = polyvander(x, lmax)
- else:
- deg = np.sort(deg)
- lmax = deg[-1]
- order = len(deg)
- van = polyvander(x, lmax)[:, deg]
- # set up the least squares matrices in transposed form
- lhs = van.T
- rhs = y.T
- if w is not None:
- w = np.asarray(w) + 0.0
- if w.ndim != 1:
- raise TypeError("expected 1D vector for w")
- if len(x) != len(w):
- raise TypeError("expected x and w to have same length")
- # apply weights. Don't use inplace operations as they
- # can cause problems with NA.
- lhs = lhs * w
- rhs = rhs * w
- # set rcond
- if rcond is None:
- rcond = len(x)*np.finfo(x.dtype).eps
- # Determine the norms of the design matrix columns.
- if issubclass(lhs.dtype.type, np.complexfloating):
- scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
- else:
- scl = np.sqrt(np.square(lhs).sum(1))
- scl[scl == 0] = 1
- # Solve the least squares problem.
- c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
- c = (c.T/scl).T
- # Expand c to include non-fitted coefficients which are set to zero
- if deg.ndim == 1:
- if c.ndim == 2:
- cc = np.zeros((lmax + 1, c.shape[1]), dtype=c.dtype)
- else:
- cc = np.zeros(lmax + 1, dtype=c.dtype)
- cc[deg] = c
- c = cc
- # warn on rank reduction
- if rank != order and not full:
- msg = "The fit may be poorly conditioned"
- warnings.warn(msg, pu.RankWarning)
- if full:
- return c, [resids, rank, s, rcond]
- else:
- return c
- def polycompanion(c):
- """
- Return the companion matrix of c.
- The companion matrix for power series cannot be made symmetric by
- scaling the basis, so this function differs from those for the
- orthogonal polynomials.
- Parameters
- ----------
- c : array_like
- 1-D array of polynomial coefficients ordered from low to high
- degree.
- Returns
- -------
- mat : ndarray
- Companion matrix of dimensions (deg, deg).
- Notes
- -----
- .. versionadded:: 1.7.0
- """
- # c is a trimmed copy
- [c] = pu.as_series([c])
- if len(c) < 2:
- raise ValueError('Series must have maximum degree of at least 1.')
- if len(c) == 2:
- return np.array([[-c[0]/c[1]]])
- n = len(c) - 1
- mat = np.zeros((n, n), dtype=c.dtype)
- bot = mat.reshape(-1)[n::n+1]
- bot[...] = 1
- mat[:, -1] -= c[:-1]/c[-1]
- return mat
- def polyroots(c):
- """
- Compute the roots of a polynomial.
- Return the roots (a.k.a. "zeros") of the polynomial
- .. math:: p(x) = \\sum_i c[i] * x^i.
- Parameters
- ----------
- c : 1-D array_like
- 1-D array of polynomial coefficients.
- Returns
- -------
- out : ndarray
- Array of the roots of the polynomial. If all the roots are real,
- then `out` is also real, otherwise it is complex.
- See Also
- --------
- chebroots
- Notes
- -----
- The root estimates are obtained as the eigenvalues of the companion
- matrix, Roots far from the origin of the complex plane may have large
- errors due to the numerical instability of the power series for such
- values. Roots with multiplicity greater than 1 will also show larger
- errors as the value of the series near such points is relatively
- insensitive to errors in the roots. Isolated roots near the origin can
- be improved by a few iterations of Newton's method.
- Examples
- --------
- >>> import numpy.polynomial.polynomial as poly
- >>> poly.polyroots(poly.polyfromroots((-1,0,1)))
- array([-1., 0., 1.])
- >>> poly.polyroots(poly.polyfromroots((-1,0,1))).dtype
- dtype('float64')
- >>> j = complex(0,1)
- >>> poly.polyroots(poly.polyfromroots((-j,0,j)))
- array([ 0.00000000e+00+0.j, 0.00000000e+00+1.j, 2.77555756e-17-1.j])
- """
- # c is a trimmed copy
- [c] = pu.as_series([c])
- if len(c) < 2:
- return np.array([], dtype=c.dtype)
- if len(c) == 2:
- return np.array([-c[0]/c[1]])
- m = polycompanion(c)
- r = la.eigvals(m)
- r.sort()
- return r
- #
- # polynomial class
- #
- class Polynomial(ABCPolyBase):
- """A power series class.
- The Polynomial class provides the standard Python numerical methods
- '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
- attributes and methods listed in the `ABCPolyBase` documentation.
- Parameters
- ----------
- coef : array_like
- Polynomial coefficients in order of increasing degree, i.e.,
- ``(1, 2, 3)`` give ``1 + 2*x + 3*x**2``.
- domain : (2,) array_like, optional
- Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
- to the interval ``[window[0], window[1]]`` by shifting and scaling.
- The default value is [-1, 1].
- window : (2,) array_like, optional
- Window, see `domain` for its use. The default value is [-1, 1].
- .. versionadded:: 1.6.0
- """
- # Virtual Functions
- _add = staticmethod(polyadd)
- _sub = staticmethod(polysub)
- _mul = staticmethod(polymul)
- _div = staticmethod(polydiv)
- _pow = staticmethod(polypow)
- _val = staticmethod(polyval)
- _int = staticmethod(polyint)
- _der = staticmethod(polyder)
- _fit = staticmethod(polyfit)
- _line = staticmethod(polyline)
- _roots = staticmethod(polyroots)
- _fromroots = staticmethod(polyfromroots)
- # Virtual properties
- nickname = 'poly'
- domain = np.array(polydomain)
- window = np.array(polydomain)
|