Polynomials and MATLAB

 

 

 

Basic terminology and ideas:

 

·       f(x) is a polynomial function of the variable x.

·       the polynomial is of degree or order n

·       the a1….an are the polynomial coefficients

·       f(x) has a numerical value for any numerical value of x

·       the roots of f(x) are the values of x that make f(x)=0

·       a plot of f(x) vs. x shows these roots where the plot of f(x) crosses the x-axis

·       these roots depend on the coefficients of f(x)

 

 

In MATLAB:

 

·       Polynomials are described by using a row vector of the coefficients of the polynomial beginning with the highest power of x and inserting zeros for “missing” terms:

 

 

f=[9 –5 3 7];

g=[6  –1 2];

 

 

·       We can add and subtract polynomial functions in MATLAB. To do this we must “zero” pad the polynomials so that their row vector representations are the same length:

 

 

h=f+[0 g];

 

 

·       We can multiply and divide polynomials by using the conv(.) and deconv(.) functions in MATLAB:

 

 

y=conv(f,g)

r=deconv(y,f)

 

·       We can evaluate a polynomial at any value of x:

 

 

p=[1 3 -4];

x=[-5:.1:5];

px=polyval(p,x);

plot(x,px);grid

 

·       Alternatively, we can evaluate the polynomial directly as:

 

 px2=p(1)*x.^2+p(2)*x+p(3)      

 

·       We can calculate the roots of p(x) and they should correspond to the zero axis crossings on the plot of p(x) as shown in Figure 1:

 

rpx=roots(p)

 

 

 

 

Figure 1: A plot of px(x) vs. x . the roots of p(x) are evident in the plot.

 

 

·       From a polynomial’s roots, we can determine the polynomial using ‘poly’:

 

poly(rpx)

 

Polynomial Curve Fitting

 

·       Often experimental data is collected and you may suspect that the underlying relationship between the variables you are observing is described by a polynomial.  In this case, you can use polyfit( ) to find the best polynomial for your data.

 

·       Let’s say you have a spring to be used in an automotive suspension.  You need to understand how much compression you get for various forces applied.  You do this by putting a variety of masses on the spring and measure the displacement of the spring.

 

 

Measurement

1

2

3

4

5

Mass (kg)

10

100

200

300

400

Force=mass x 9.8 (Newtons)

98

980

1960

2940

3920

Displacement (meters)

0.0092

0.0953

0.1802

0.3012

0.3909

 

·       Let’s enter these data into MATLAB, plot it and find a good polynomial to describe the “underlying” relationship between force and spring displacement.

 

M=[10,100,200,300,400];

F=M*9.8;

D=[0.0092,0.0953,0.1802,0.3012,0.3909];

plot(D,F, 'r*');

xlabel('Displacement (m) ');

ylabel('Force (N) ');

title('Experimental Spring Data');

 

 

>> help polyfit

 

 POLYFIT Fit polynomial to data.

    POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of

    degree N that fits the data, P(X(I))~=Y(I), in a least-squares sense.

 

p=polyfit(D,F,1)  % try degree 1 polynomial, data looks like a line!

 

d=linspace(.0092,.3909,1000);

f=polyval(p,d);

 

plot(D,F,'r*',d,f, 'b-');

legend('Experimental Data', 'Curve Fitting',2)    % The 2 puts legend in upper left

xlabel('Displacement (m) ');

ylabel('Force (N) ');

title('Experimental Spring Data With Curve Fit');

 

·       It does appear that a straight line (degree 1 polynomial) does a good job describing our experimental data.  The slope of the line is something called the “spring constant.”

 

·        Based on experiments like this, it is understood that an ideal spring follows this relationship: f=kd, where f is force in Newtons, k is the spring constant and d is the displacement of the spring in meters.

 

·       Thus, our spring constant is estimated to be k=p(1).  Now we can predict how the spring will respond to most any force.  Keep in mind that when the force gets sufficiently large, the spring will be fully compressed and will no longer obey this idealized linear relationship.

 

Interactive Curve Fitting

 

·       Try this…

 

     clear; close all; clc

M=[10,100,200,300,400];

F=M*9.8;

D=[0.0092,0.0953,0.1802,0.3012,0.3909];

plot(D,F, 'r*');

xlabel('Displacement (m) ');

ylabel('Force (N) ');

title('Experimental Spring Data');

 

·       Now in the figure window, select “Tools” then “Basic Fitting”.  Click on the right arrow at the bottom and then try selecting types of curves to fit to your data in the upper left hand window.