MATLAB Examples

Tutorial on Polynomials in MATLAB

This is the code documented in the tutorial on Polynomials in MATLAB. (This published MATLAB code, only presents the code implementation. The PDF document on the tutorial also includes the theory behind the implementation).

(c) Paul O'Leary and Matthew Harker 2013

Contents

clear up the workspace and close all figures before starting the computation

close all;
clear all;
setUpGraphics;

Simple Polynomial

\cellName{defCoeffs}

define the coefficients o fthe polynomial

yCfs = [1,1,0,-0.1];

\cellName{defXvals}

Define the x value for which the polynomial should be evaluated

n = 120;
x = linspace(-1,1,n)';

\cellName{polyVal}

Evaluate the polynomial which has the coefficients yCfs at the x points defined above

y = polyval( yCfs, x );

\cellName{plotPoly1}

fig1 = figure;
plot(x,y,'k');
xlabel('$$x$$');
ylabel('$$y(x)$$');
grid on;
%
%\caption{Example of the polynomial $y(x)$}

Polynomial roots

\cellName{rootsOfY}

rootsY = roots( yCfs )
%
% and add the roots to the previous plot
%
figure(fig1);
hold on;
plot( rootsY, zeros(size(rootsY)), 'ko', 'markerFaceColor', 'w');
%\caption{Example of the polynomial $y(x)$ with the roots added}
rootsY =
   -0.8670
   -0.4126
    0.2796

Multiple Roots

\cellName{multipleRoots}

There may be numerical errors when computing the roots of a polynomial which has multiple roots, i.e. more than one root at exactly the same value.

%
% Define multiple roots for a polynomial, in this example 4 roots at the
% x value of 1.
%
defR = [1, 1, 1, 1];
%
% Compute the corresponding polynomial
%
polyMRs = poly( defR );
%
% Use roots to compute the roots of this poylnomial, this should return
% exactly three 1's.
%
compR = roots( polyMRs )
%
compR =
   1.0002          
   1.0000 + 0.0002i
   1.0000 - 0.0002i
   0.9998          

Polynomial from roots

\cellName{polyFromRoots}

Consider the polynomial with roots at x = 0 and x = 1.

rs = [1, 0];
pCfs = poly( rs )
%
pCfs =
     1    -1     0

\cellName{stationaryPoints}

Compute the coefficients of the derivative of the polynomial

dyCfs = polyder( yCfs );
%
% determine the x locations of the stationary points
%
xStationary = roots( dyCfs );
%
% and now evaluate the original polynomial at the x locations
% of the stationary points to obtain the corresponding y values.
%
yStationary = polyval( yCfs, xStationary);
%
% Then add the points to the plot.
%
figure(fig1);
plot( xStationary, yStationary, 'ko', 'markerFaceColor', 'k');
%\caption{Example of the polynomial $y(x)$ with the roots and
% astationary points marked}

\cellName{multPolys}

Define the coefficients of the polynomials g(x) and f(x)

gCfs = [-1,0.1,1];
fCfs = [1, -0.2];
%
% compute the convolution of g(x) and f(x) to obtain h(x)
% and display the result
%
hCfs = conv( gCfs, fCfs )
hCfs =
   -1.0000    0.3000    0.9800   -0.2000

\cellName{convMtx}

We also compute the convolution matrix here, it will not be used until later in this program. Note the use of the transpose to ensure the correct dimensioning of the matrices

Cg = convmtx( gCfs', length( fCfs ) )
Cf = convmtx( fCfs', length( gCfs ) )
%
Cg =
   -1.0000         0
    0.1000   -1.0000
    1.0000    0.1000
         0    1.0000
Cf =
    1.0000         0         0
   -0.2000    1.0000         0
         0   -0.2000    1.0000
         0         0   -0.2000

\cellName{polyPlots}

Evaluate the polynomials and plot the resulting curves

x = linspace(0,1,n);
%
gx = polyval( gCfs, x );
fx = polyval( fCfs, x );
hx = polyval( hCfs, x );
%
fig2 = figure;
plot( x, gx, 'r');
grid on;
hold on;
plot( x, fx, 'b');
plot( x, hx, 'k');
xlabel('$$x$$');
ylabel('$$ f(x), g(x), h(x) $$');
legend('$$g(x)$$','$$f(x)$$','$$h(x)$$');
% \caption{Plot o fthe values of $g(x)$, $f(x)$ and $h(x)$.
% The coefficients of the polynomial $h(x)$ were generated using convolution.}
%

Weight and Tower Example

\cellName{towerEqns}

Define the known constants

mass = 10;
height = 30;
g = 9.8;
%
% Acceleration
%
aCfs = -g;
%
% Velocity
%
vCfs = polyint( aCfs );
%
% position
%
pCfs = polyint( vCfs, height );

\cellName{evalEqnsMotion}

td = max(roots( pCfs ));
%
% generate the time vector
%
nTimePoints = 200;
time = linspace(0,td,nTimePoints )';
%
% Now evaluate the polynomials for velocity and time
%
vt = polyval( vCfs, time );
pt = polyval( pCfs, time );
%
% and compute the kenetic energy. Please note the ".*"
% notation, i.e., each entry in the vector is squared.
%
ke = (mass * vt.^2 ) / 2;

\cellName{plotEqnsMotion}

Plot the results

fig3 = figure;
subplot(3,1,1);
plot( time, vt, 'r');
grid on;
ylabel('$$v(t)$$');
%
subplot(3,1,2);
plot( time, pt, 'r');
grid on;
ylabel('$$p(t)$$');
%
subplot(3,1,3);
plot( time, ke, 'r');
grid on;
ylabel('$$k_e(t)$$');
%
xlabel('Time [s]');
%
%\caption{The velovity $v(t)$, position $p(t)$ and
% kinetic energy $k_e(t)$ as a function of time for the
% weight dropped from the tower.}
%

Polynomial Approximation: An Example

\cellName{loadFitData}

load fitDataSet;
%
fig4 = figure;
plot(x,y,'k.');
hold on;
grid on;
xlabel('$$x$$');
ylabel('$$y(x)$$');
%
% \caption{Test data for a simple fit.}

\cellName{simpleFit}

Fit the polynomial

degree = 2;
yCfs = polyfit( x, y, degree );
%
% Evaluate the polynomial
%
yx = polyval( yCfs, x );
%
figure(fig4);
plot( x, yx, 'k');
text(0.1, 225, ['$$y(x) = ',poly2latex(yCfs,'x'),'$$'], 'BackgroundColor', 'w' );
%
% \caption{Simple polynomial fit with the result of the fit.}
signs =
    -1    -1     1
ans =
    -1

Covariance of the Polynomial Coefficients

\cellName{polyCov}

Use the function polyfitcov to obtain the coefficients and their variance

[yCfs, S, Ly] = polyfitcov( x, y, degree )
%
%
yCfs =
 -272.0261  -25.3260  502.8972
S = 
        R: [3x3 double]
       df: 97
    normr: 121.0598
Ly =
  261.3718 -261.3718   43.1220
 -261.3718  279.1433  -52.0077
   43.1220  -52.0077   13.0681

Evaluation a Polynomial with Confidence Estimation

\cellName{polyFitDelta}

% Calling polyfit and polyval with the necessary paramaters to generate
% the estimation bounds
%
[yCfs, S] = polyfit( x, y, degree );
[yx, delat] = polyval( yCfs, x, S );
%
fig5 = figure;
plot(x,y,'k.');
hold on;
plot( x, yx, 'k');
plot( x, yx + delat, 'r');
plot( x, yx - delat, 'r');
grid on;
xlabel('$$x$$');
ylabel('$$y(x)$$');
%
% \caption{Polynomial fit with estimation of the $50 \%$ estimation bound.}

Shift and Scale x

\cellName{polyShift1}

p1 = polyfit( x, y, degree );
[p2, S, mu] = polyfit( x, y, degree );
%
poly1 = poly2latex( p1, 'x' )
poly2 = poly2latex( p2, 'x' )
signs =
    -1    -1     1
ans =
    -1
poly1 =
 - 272.0261 \, x^{2} - 25.326 \, x + 502.8972
signs =
    -1    -1     1
ans =
    -1
poly2 =
 - 23.3604 \, x^{2} - 87.1376 \, x + 422.2277

\cellName{polyShift2}

p3 = polytrans( p2, mu );
poly3 = poly2latex( p3, 'x' )
signs =
    -1    -1     1
ans =
    -1
poly3 =
 - 272.0261 \, x^{2} - 25.326 \, x + 502.8972

Example Polymult with Covariance

\cellName{fitMult}

%
% Fit the polynomial and determine the covariance of the coefficients
%
degree = 2;
[yCfs, S, Ly] = polyfitcov( x, y, degree );
%
% Define the equation relating price/vol to profit/vol, note these
% coefficients are free from error
%
kf = 0.2;
gCfs = [1, -kf];
%
% Multiply the two polynomials
%
[pCfs, Lp] = polymult( yCfs, gCfs, Ly, []);
%
% Find the stationary points
%
xStat = max(roots( polyder( pCfs )));
yStat = polyval( pCfs, xStat );
%
% Now evaluate the polynomial for specific values of x and determine the
% covariance matrix for the final result
%
px = polyval( pCfs, x );
V = polyvander( x, length(pCfs) - 1 );
%
Lpx = V * Lp * V';
%
% Use the diagonal of Lps as a 1 sigma estimate
%
delta = sqrt( diag( Lpx ));

\cellName{fitMultDisp}

fig6 = figure;
plot( x, px, 'k');
hold on;
plot( x, px - delta*3, 'r');
plot( x, px + delta*3, 'r');
plot( xStat, yStat, 'ko', 'MarkerFaceColor', 'k');
grid on;
xlabel('Price per liter $$(x)$$');
ylabel('Total profit $$p(x)$$');
axis([0.2,1,0,200]);
latexPoly = poly2latex( pCfs, 'x');
title( ['$$p(x) = ',latexPoly,'$$'] );
text(0.55, 100, ['$$x_{max} = ',num2str(xStat),...
    '\, y_{max} = ',num2str(yStat),'$$'], ...
    'BackGroundColor','w');
%
% \caption{Polynomial relating prive per liter to profit per liter
% with the upper and lower error bound estimates.}
signs =
    -1     1     1    -1
ans =
    -1

Fitting Polynomials with Homogeneous Value Constraints

\cellName{hConstraints}

Load the example data which contains x, y data and the position of the root i.e. xr

load constPoly;
%
% The coefficients of the constraining polynomial are computed as
%
cr = poly( xr );
%
% Determine the number of free coefficients
%
d = 3;
m = d - length( xr ) + 1;
%
% Setup the convolution matrix
%
Cr = convmtx( cr', m );
%
% generate a Vandermonde matrix and the constrained basis functions
%
V = polyvander( x, d );
Vr = V * Cr;
%
% Use Qr decomposition to performe the fitt
%
[Q,R] = qr( Vr, 0 );
cm = R\(Q'*y);
%
% Compute the final polynomial
%
cf = Cr * cm;
%

\cellName{hCPlot}

xs = linspace(0,1,100)';
yf = polyval( cf, xs );
%
fig1 = figure;
plot( xs, yf, 'k');
xlabel('Price per liter $$(x)$$');
ylabel('Total profit $$p(x)$$');
hold on;
plot(x,y,'ko','MarkerFaceColor','w');
plot(xr,0,'ko','MarkerFaceColor','k');
grid on;
range = axis;
plot( range(1:2), [0,0], 'k');
%
% \caption{Fitting the polynomial with homogeneous constraint $ y(0.2) = $ using matrix
% convolution to implement the constrained fit.}

Alternative Solution to Constrained Fitting

\cellName{hConstraints2}

d = 3;
%
% Generate a linear function in x shich fulfills the constraint
%
xc = x - xr;
%
% The coefficients for this polynomial are
%
cxc = [1, -xr];
%
% The corresponding convolution matrix is
%
Cxc = convmtx( cxc', d );
%
% Since there is 1 constraing generate a vandermonde matrix with one degre
% less than originally required
%
V = polyvander( x, d - 1);
%
% Compute the constrained set of basis functions
%
Vc = diag( xc) * V;
%
% Compute a fit
%
[Q,R] = qr( Vc, 0 );
cc = R\(Q'*y);
yf = Vc * cc;
%
% compute the complete polynomial cofficients
%
ccf = Cxc * cc;

\cellName{hCPlot2}

xs = linspace(0,1,100)';
yf = polyval( ccf, xs );
%
fig1 = figure;
plot( xs, yf, 'k');
xlabel('Price per liter $$(x)$$');
ylabel('Total profit $$p(x)$$');
hold on;
plot( x, y, 'ko', 'MarkerFaceColor', 'w');
grid on;
range = axis;
plot( range(1:2), [0,0], 'k');
plot( xr, 0, 'ko', 'MarkerFaceColor', 'k');
%
% \caption{Fitting the polynomial with homogeneous constraint $ y(0.2) = $ using
% premultiplication by a polynomial with the desired roots.}