" /> DEMOPOLYNOM Short demo of the polynom toolbox

DEMOPOLYNOM Short demo of the polynom toolbox

Contents

format compact short
setround(0)                           % set rounding to nearest

Definition of a univariate polynomial

The simplest way to generate a univariate polynomial is (like in Matlab) by

p = polynom([1 -3 0 4])

polynom p[x] = 
    1.0000  x^3  
   -3.0000  x^2  
    4.0000

Note we use the German word Polynom, saving three letters :)

It generates a polynomial (of the INTLAB data type "polynom") with coefficients 1, -3, 0 and 4. Note that the coefficient corresponding to the highest exponent is specified first, and that the default dependent variable is "x". Another variable can be specified explicitly, e.g. by

q = polynom([1 0 -2],'y')

polynom q[y] = 
    1.0000  y^2  
   -2.0000

Access to coefficients, exponents and variables I

There is direct access to the vector of polynomial coefficients (starting with the largest exponent), the vector of exponents and the independent variable in use:

coeff = q.c
expon = q.e
vars = q.v

coeff =
     1     0    -2
expon =
     2
vars =
    'y'

The polynomial may also be specified by the individual coefficients and exponents. The polynomial p, for example, is also generated as follows:

polynom([1 -3 4],[3 2 0])

polynom ans[x] = 
    1.0000  x^3  
   -3.0000  x^2  
    4.0000

Definition of a multivariate polynomial

A multivariate polynomial is generated by specifying coefficients and corresponding exponents. An example is

P = polynom([-3 4 9],[2 3;4 0;2 2],{'a' 'b'})

polynom P[a,b] = 
    4.0000  a^4       
   -3.0000  a^2  b^3  
    9.0000  a^2  b^2

Random polynomials

A multivariate polynomial may generated randomly by

Q = randpoly(4,2)

polynom Q[x1,x2] = 
   -1.1512  x1^4  x2^4  
    0.6932  x1^4        
    0.1707  x1^3  x2^2  
   -0.4886  x1^2  x2^3  
   -0.7017        x2^2

where the first parameter specifies the degree and the second the number of variables. Note that the variables are "x1", "x2", ... by default. This may be changed by specifying other variable names explicitly:

QQ = randpoly(4,2,{'var1' 'var2'})

polynom QQ[var1,var2] = 
   -0.8400  var1^2  var2^3  
   -0.5652  var1^2          
   -0.2322  var1    var2^3  
   -0.6180  var1    var2^2  
   -0.2553  var1    var2

Access to coefficients, exponents and variables II

As before there is also direct access to the polynomial coefficients, the exponents and the independent variables for multivariate polynomials:

coeff = QQ.c
expon = QQ.e
vars = QQ.v

coeff =
   -0.8400
   -0.5652
   -0.2322
   -0.6180
   -0.2553
expon =
     2     3
     2     0
     1     3
     1     2
     1     1
vars =
  1×2 cell array
    {'var1'}    {'var2'}

Display of polynomials

Univariate polynomials may be displayed in dense or sparse mode, for example

format upolyvector
p

polynom p[x]  = 
     1    -3     0     4

format upolysparse
p

polynom p[x] = 
    1.0000  x^3  
   -3.0000  x^2  
    4.0000

Operations between polynomials

Operations between univariate polynomials are as usual

p, 3*p+1

polynom p[x] = 
    1.0000  x^3  
   -3.0000  x^2  
    4.0000       
polynom ans[x] = 
    3.0000  x^3  
   -9.0000  x^2  
   13.0000

and may produce multivariate polynomials if not depending on the same variable:

q, p+q

polynom q[y] = 
    1.0000  y^2  
   -2.0000       
polynom ans[x,y] = 
    1.0000  x^3       
   -3.0000  x^2       
    1.0000       y^2  
    2.0000

Interval polynomials

Interval polynomials are specified in the same way as before. Consider, for example (taken from Hansen/Walster: Sharp Bounds for Interval Polynomial Roots, Reliable Computing 8(2) 2002)

format infsup
r = polynom([infsup(1,2) infsup(-4,2) infsup(-3,1)])

intval polynom r[x] = 
[    1.0000,    2.0000]   x^2  
[   -4.0000,    2.0000]   x    
[   -3.0000,    1.0000]

The polynomial may be displayed using other interval formats, for example

format midrad
r

intval polynom r[x] = 
<    1.5000,   0.5000>   x^2  
<   -1.0000,   3.0000>   x    
<   -1.0000,   2.0000>

Plot of polynomials

The following plots the lower and upper bound polynomial within root bounds:

plotpoly(r)

or within specified bounds:

plotpoly(r,-2,2)

Access of coefficients I

In contrast to Matlab, coefficients of INTLAB polynomials are set and accessed as in mathematics:

q = p+1
coeff3 = q(3)
q(0) = -2
q(0:2) = 4.7

polynom q[x] = 
    1.0000  x^3  
   -3.0000  x^2  
    5.0000       
coeff3 =
     1
polynom q[x] = 
    1.0000  x^3  
   -3.0000  x^2  
   -2.0000       
polynom q[x] = 
    1.0000  x^3  
    4.7000  x^2  
    4.7000  x    
    4.7000

Access of coefficients II

Access of coefficients for multivariate polynomials works the same way by specifying the position for the individual variables:

P = polynom([-3 4 9],[2 3;4 0;2 2],{'a' 'b'})
coeff23 = P(2,3)
P(1,4) = -9

polynom P[a,b] = 
    4.0000  a^4       
   -3.0000  a^2  b^3  
    9.0000  a^2  b^2  
coeff23 =
    -3
polynom P[a,b] = 
    4.0000  a^4       
   -3.0000  a^2  b^3  
    9.0000  a^2  b^2  
   -9.0000  a    b^4

Subpolynomials

Subpolynomials may be accessed by specifying certain unknowns as []. This corresponds to a distributive representation of the polynomial:

P
Q = P(2,[])

polynom P[a,b] = 
    4.0000  a^4       
   -3.0000  a^2  b^3  
    9.0000  a^2  b^2  
   -9.0000  a    b^4  
polynom Q[b] = 
   -3.0000  b^3  
    9.0000  b^2

Polynomial evaluation

There are two (equivalent) possibilities of polynomial evaluation, by polyval or by {}:

p = polynom([1 -3 0 4])
polyval(p,2)
p{2}

polynom p[x] = 
    1.0000  x^3  
   -3.0000  x^2  
    4.0000       
ans =
     0
ans =
     0

Interval polynomial evaluation

Of course, verified bounds are obtained in the well-known ways:

polyval(intval(p),2)
p{intval(2)}

intval ans = 
<    0.0000,   0.0000> 
intval ans = 
<    0.0000,   0.0000>

Polynomial evaluation for multivariate polynomials works the same way:

P = polynom([-3 4 9],[2 3;4 0;2 2],{'a' 'b'})
polyval(P,2,3)
P{2,intval(3)}

polynom P[a,b] = 
    4.0000  a^4       
   -3.0000  a^2  b^3  
    9.0000  a^2  b^2  
ans =
    64
intval ans = 
<   64.0000,   0.0000>

Evaluation of subpolynomials

In addition, evaluation of sub- (or coefficient-) polynomials is possible by specifying certain unknowns as []. Unknowns specified by [] are still treated as independent variables. In this case the argument list must be one cell array:

polyval(P,{2,[]})
P{{[],intval(3)}}

polynom ans[b] = 
  -12.0000  b^3  
   36.0000  b^2  
   64.0000       
intval polynom ans[a] = 
<    4.0000,   0.0000>   a^4

Derivatives of polynomials

First and higher polynomial derivatives are calculated by

p
p'
pderiv(p,2)

polynom p[x] = 
    1.0000  x^3  
   -3.0000  x^2  
    4.0000       
polynom ans[x] = 
    3.0000  x^2  
   -6.0000  x    
polynom ans[x] = 
    6.0000  x    
   -6.0000

or, for multivariate polynomials, by specifying the variable:

P = polynom([-3 4 9],[2 3;4 0;2 2],{'a' 'b'})
pderiv(P,'a')
pderiv(P,'b',2)

polynom P[a,b] = 
    4.0000  a^4       
   -3.0000  a^2  b^3  
    9.0000  a^2  b^2  
polynom ans[a,b] = 
   16.0000  a^3       
   -6.0000  a    b^3  
   18.0000  a    b^2  
polynom ans[a,b] = 
  -18.0000  a^2  b    
   18.0000  a^2

Bernstein polynomials

A simple application is the computation of Bernstein coefficients. Consider

P = polynom([2 -3 0 3 1 -2])

polynom P[x] = 
    2.0000  x^5  
   -3.0000  x^4  
    3.0000  x^2  
    1.0000  x    
   -2.0000

Suppose, we wish to expand the polynomial in the interval [-1,1]:

B = bernsteincoeff(P,infsup(-1,1))

polynom B[x] = 
    1.0000  x^5  
   -1.0000  x^4  
   -1.0000  x^3  
   -5.4000  x^2  
    1.8000  x    
   -5.0000

Note that Bernstein coefficients are stored in vector (that is a change with Version 12 of INTLAB). This new implementation follows

[TG] J. Titi, J. Garloff, Matrix methods for the tensorial Bernstein form,
       Applied Mathematics and Computation 346, p.254-271, 2019.

It is much faster than my previous implementation! That is in particular beneficial in the taylormodel toolbox.

Polynomial evaluation and Bernstein polynomials

The convex hull of the Bernstein points B contains the convex hull of the polynomial:

plotpoly(P,-1,1);
hold on
plotbernstein(B,-1,1)
hold off

This picture is not untypical. The Bernstein points overestimate the true range, but sometimes not too much. To obtain a true inclusion of the range, we perform the computation with verified bounds:

format infsup
B = bernsteincoeff(intval(P),infsup(-1,1))
X = P{infsup(-1,1)},
Y = infsup(min(B.c.inf),max(B.c.sup))

intval polynom B[x] = 
[    0.9999,    1.0001]   x^5  
[   -1.0001,   -0.9999]   x^4  
[   -1.0001,   -0.9999]   x^3  
[   -5.4001,   -5.3999]   x^2  
[    1.7999,    1.8001]   x    
[   -5.0000,   -5.0000]        
intval X = 
[  -11.0000,    7.0000] 
intval Y = 
[   -5.4001,    1.8001]

From the picture we read the true range [-5,1] which is slightly overestimated by Y computed by the Bernstein approach. The same principle can be applied to multivariate polynomials.

Inclusion of roots of polynomials

Roots of a univariate polynomial can approximated and included. Consider a polynomial with roots 1,2,...,7:

format long
p = polynom(poly(1:7))
roots(p)                  % approximations of the roots

polynom p[x] = 
  1.0e+004 *
   0.00010000000000  x^7  
  -0.00280000000000  x^6  
   0.03220000000000  x^5  
  -0.19600000000000  x^4  
   0.67690000000000  x^3  
  -1.31320000000000  x^2  
   1.30680000000000  x    
  -0.50400000000000       
ans =
   6.999999999999008
   6.000000000002801
   4.999999999997638
   4.000000000000223
   3.000000000000553
   1.999999999999781
   1.000000000000020

Based on some approximation (in this case near 4.1), verified bounds for a root are obtained by

verifypoly(p,4.1)

intval ans = 
[   3.99999999999739,   4.00000000000329]

Note that the accuracy of the bounds is of the order of the (usually unknown) sensitivity of the root.

Inclusion of clustered or multiple roots of polynomials

The routine "verifypoly" calculates verified bounds for multiple roots as well. Consider the polynomial with three 4-fold roots at x=1, x=2 and x=3:

format short midrad
p = polynom(poly([1 1 1 1 2 2 2 2 3 3 3 3])); roots(p)

ans =
   3.0008 + 0.0011i
   3.0008 - 0.0011i
   2.9992 + 0.0004i
   2.9992 - 0.0004i
   2.0066 + 0.0000i
   2.0000 + 0.0066i
   2.0000 - 0.0066i
   1.9934 + 0.0000i
   1.0016 + 0.0000i
   1.0000 + 0.0016i
   1.0000 - 0.0016i
   0.9984 + 0.0000i

Based on some approximation (in this case 2.001), verified bounds for a multiple root are obtained by

verifypoly(p,2.001)

intval ans = 
<    2.0000,   0.0055>

Note that the accuracy of the bounds is of the order of the sensitivity of the root, i.e. of the order eps^(1/4) = 1.2e-4. The multiplicity of the root is determined as well (for details, see "verifypoly").

[X,k] = verifypoly(p,2.999)

intval X = 
<    3.0000,   0.0050> 
k =
     4

Quality of the computed bounds

One may argue that the previous inclusion is rather broad. But look at the graph of the polynomial near x=3:

p = polynom(poly([1 1 1 1 2 2 2 2 3 3 3 3]));
plotpoly(p,2.99,3.01)
title('Floating point evaluation of p near x=3')

Quality of the computed bounds for coefficients with tolerances

Due to rounding and cancellation errors the accuracy of the inclusion is about optimal. Things are even more drastic when the polynomial coefficients are afflicted with tolerances. Put smallest possible intervals around the coefficients of p and look at the graph. Between about 2.993 and 3.007 the lower and upper bound enclose zero.

P = polynom(midrad(p.c,1e-16));
plotpoly(P,2.99,3.01)
title('Interval polynomial P near x=3')

Sylvester matrix

The Sylvester matrix of two polynomials p,q is singular if and only if p and q have a root in common. Therefore, the Sylvester matrix of p and p' may determine whether p has multiple roots. Note, however, that this test is not numerically stable.

format short
p = polynom(poly([2-3i 2-3i randn(1,3)]))       % polynomial with double root 2-3i
roots(p)
S = sylvester(p);         % Sylvester matrix of p and p'
format short e
svd(S)

polynom p[x] = 
   1.0000 +   0.0000i  x^5  
  -5.4585 +   6.0000i  x^4  
  -0.6680 -  20.7508i  x^3  
  13.1249 +   8.4901i  x^2  
   8.2104 +  16.9721i  x    
   0.8759 +   2.1022i       
ans =
   2.0000 - 3.0000i
   2.0000 - 3.0000i
   2.1832 + 0.0000i
  -0.5883 + 0.0000i
  -0.1364 + 0.0000i
ans =
   1.3202e+02
   1.0042e+02
   7.2270e+01
   4.7457e+01
   2.6990e+01
   1.1919e+01
   1.6622e+00
   8.1720e-01
   1.3388e-15

Predefined polynomials

There are a number of predefined polynomials such as Chebyshev, Gegenbauer, Hermite polynomials etc.

Enjoy INTLAB

INTLAB was designed and written by S.M. Rump, head of the Institute for Reliable Computing, Hamburg University of Technology. Suggestions are always welcome to rump (at) tuhh.de


Published with MATLAB® R2020a