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

```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.} ```