Download Polynomials in MATLAB: Roots, Multiplication, Addition, Division, and Differentiation and more Study notes Engineering in PDF only on Docsity!
Chapter 18 Polynomials
18.1 Roots
MATLAB performs numerous operations involving polynomials. The most
common situation dealing with polynomials is finding its roots. A polynomial is
uniquely defined by its coefficients. MATLAB represents a polynomial by a vector
containing the coefficients of the powers in descending order. For example, the
polynomial
1
n n
p x a xn an x a x a
−
is saved in MATLAB as a vector p where p = [ an , an − 1 ,......, a 1 , a 0 ]. If a coefficient is
zero, it must appear in the vector. After the coefficient vector is defined, the roots of the
polynomial are returned in a column vector from a call to 'roots(p)'.
The roots of a polynomial are found in the example below. The polynomial is
evaluated at the roots (where it should be zero) and then plotted to verify it passes
through the real roots.
Example 18.1.
p=[1,-10,42,-98,65] % Define coefficient vector for polynomial % p(x) = x^4 - 10x^3 + 42x^2 - 98x + 65 proots=roots(p); % Find roots of p(x) x=proots % Set x equal to row vector of roots of p(x)= y=p(1)x.^4 + p(2)x.^3 + p(3)x.^2 + p(4)x + p(5) % Evaluate p(x) % at each of its 4 roots xi=linspace(0,10,250); % Create xi vector for evaluating yi=p(xi) yi=p(1)xi.^4 + p(2)xi.^3 + p(3)xi.^2 + p(4)xi + p(5); % yi=p(xi) plot(xi,yi) % Plot y=p(x) vs x hold on axis([0 6 -50 50]) % Set scale for plot plot([0 10],[0 0],'k') % Plot x-axis xlabel('x'), ylabel('y') title('y=p(x) vs. x')
p = 1 -10 42 -98 65
x = 2.0000 + 3.0000i 2.0000 - 3.0000i
y = 1.0e-012 * -0.0568 + 0.0568i -0.0568 - 0.0568i -0. 0
0 1 2 3 4 5 6
0
1 0
2 0
3 0
4 0
5 0
x
y
y = p ( x ) v s. x
The reverse process, i.e. obtaining the polynomial coefficients from the roots is
easily accomplished using the function 'poly'.
Example 18.1.
pcoeffs=poly(proots) % Find coeff's of polynomial with roots 'proots'
pcoeffs = 1.0000 -10.0000 42.0000 -98.0000 65.
18.3 Addition
Addition of polynomials is a bit more complicated than multiplication when the
polynomials involved are not the same order. That is, MATLAB cannot add polynomials
of different order directly because the coefficient vectors representing them are not the
same size. The solution is to pad the lower order polynomial coefficient vector with
leading zeros until the two coefficient vectors are the same size. The function
'polyadd' in Example 18.3.1 below adds two polynomials. Several test cases are
included after the function.
Example 18.3.
function h=polyadd(f,g) % This function adds two polynomials. If the polynomials are not the % same order, it pads the coefficient vector of the lower order % polynomial with zeros so the two vectors can be added.
if length(f)==length(g) % Check if polynomials are the same order h=f+g elseif length(f)>length(g) % Check if f is higher order than g diff=length(f)-length(g); gg=[zeros(1:diff) g] % Pad g with 'diff' zeros h=f+gg % Compute sum else % g is higher order than f diff=length(g)-lenght(f) ff=[zeros(1:diff) f] % Pad f with 'diff' zeros h=ff+g % Compute sum end % if length(f)==length(g)
disp(' '), disp('Case 1') g=[1 2 -3 4] % Define g(x) = x^3 + 2x^2 - 3x + 4 h=[-1 -2 3 -4] % Define h(x) = -x^3 - 2x^2 + 3x - 4 z=polyadd(g,h) % Add polynomials g and h disp(' '), disp('Case 2') f=[1 3 -2 0 5] % Define f(x) = x^4 +3x^3 -2x^2 + 5x g=[4 -1 6] % Define g(x) = 4x^2 - x + y=polyadd(f,g) % Add polynomials f and g disp(' '), disp('Case 3') a=[1] % a(x)= b=[1 0 0 0 -1] % b(x) = x^4 - 1 y=polyadd(a,b)
Case 1
g = 1 2 -3 4 h = -1 -2 3 - z = 0 0 0 0
Case 2
f = 1 3 -2 0 5 g = 4 -1 6 y = 1 3 2 -1 11
Case 3
a = 1 b = 1 0 0 0 - y = 1 0 0 0 0
Note the use of the local variable 'diff' in the function workspace. Despite the
fact a built-in MATLAB function named 'diff' exists, there is no problem assigning it a
value and refering to it in the next statement. Recall that MATLAB first looks in the
current workspace for named variables before checking if its a built-in function.
[Cq,Cr]=deconv(A,B) % Find A(x)/B(x)
- Example 18.4.
- A=[1 2 3] % Define A(x) = x^2 + 2x +
- B=[1 0 1 0 5] % Define B(x) = x^4 + x^2 +
- A =
- B =
- Cq =
- Cr =
18.5 Differentiation
Differentiation of polynomial functions is a straightforward procedure, one that
MATLAB implements with the function 'polyder'.
Example 18.5.
f(1)=3; % Coefficient of highest order term f(2)=5; % Coefficient of next highest order term f(3)=0; f(4)=-1; f(5)=2; % Coefficient of zero order term, i.e. the constant f % f(x) = 3x^4 + 5x^3 - x + 2 dfdx=polyder(f) % Find df/dx coefficient vector
f = 3 5 0 -1 2
dfdx = 12 15 0 -
Example 18.5.
f=[1 0 3 2] % f(x) = x^3 + 3x + g=[1 -4 5] % g(x) = x^2 - 4x + fg=conv(f,g) % Find f(x) ⋅⋅ g(x) coefficient vector
d_fg_dx=polyder(fg) % Differentiate the product f(x) ⋅⋅ g(x)
u1=conv(f,polyder(g)); % Find f ⋅⋅ dg/dx u2=conv(g,polyder(f)); % Find g ⋅⋅ df/dx
d_fg_dx=polyadd(u1,u2) % Differentiate f(x) ⋅⋅ g(x) using product rule
f = 1 0 3 2
g = 1 -4 5
fg = 1 -4 8 -10 7 10
d_fg_dx = 5 -16 24 -20 7
d_fg_dx = 5 -16 24 -20 7
0 5 1 0 1 5 2 0 2 5
0
2
4
6
8
1 0
1 2
1 4
1 6
1 8
F u t u r e W o r t h o f 1 6 A n n u a l $ 1 2 0 0 D e p o s i t s V s. A n n u a l I n t e r e s t R a t e i
I n t e r e s t R a t e p e r A n n u m ( % )
Future Worth ($ x 10,000)
Example 18.6.
a=conv([1 0 -1],[1 -2]) % Coefficient vector of f(x)=(x^2 -1) ⋅⋅ (x-2) f_roots=roots(a) xi=-2:0.01:3; % Values of x to evaluate f(x) at yi=polyval(a,xi); % yi=f(xi) plot(xi,yi) % Plot y=f(x) vs x hold on plot([-2,3],[0,0],'k') % Plot x-axis xlabel('x') ylabel('y') title('f(x) vs x')
a = 1 -2 -1 2
f_roots =
-1.
- 2 - 1. 5 - 1 - 0. 5 0 0. 5 1 1. 5 2 2. 5 3
- 1 2
- 1 0
- 8
- 6
- 4
- 2
0
2
4
6
8
x
y
f ( x ) v s x
In this case,
3 2
x x x
x x
x x
The 'residue' function can also be used to sum the individual terms of the
partial fraction expansion. There must be 3 input vector arguments r, p, and k (an empty
vector in the absence of a non-fractional term).
Example 18.7.
r,p [n,d]=residue(r, p, [ ]) % Combine partial fraction terms into n/d
r =
-1.
p = -2. -2. -0.
n = -0.0000 1.0000 3. d = 1.0000 5.0000 7.0000 2.
r1,p1,k [n1,d1]= residue(r1, p1, k1) % Combine r1/p1 + k1 polynomial
r1 = - p1 = - k1 = 1 2 1
n1 = 1 5 7 2 d1 = 1 3
18.8 Curve Fitting
There are a variety of ways to represent a set of empirical data. One method,
Least Squares Regression, is based on the use of low order polynomials to fit the data in
the least squares sense. A least squares fit is a unique curve with the property that the
sum of the squares of the differences between the fitted curve and the given data at the
data points is a minimum. MATLAB generates least squares polynomials for a set of
data by using the function 'polyfit'. It accepts vectors of x and y data and a scalar n
as the order of the polynomial approximation.
The target heart rate during exercise for different ages is tabulated below
Age 20 25 30 35 40 45 50 55 60 65 70
Target Heart Rate 150 146 142 139 135 131 127 124 120 116 112
The data is plotted and a first order least square polynomial is found and
presented on the same graph in the next example.
Example 18.8.
x=20:5:70 % Create vector of age values y=[150 146 142 139 135 131 127 124 120 116 112] % Create vector of % target heart rate values plot(x,y,'*') % Plot raw data a1=polyfit(x,y,1) % Find coefficients of least square line xi=linspace(20,70,250); % Create vector of x values for plotting yi=polyval(a1,xi); % Evaluate points on least square line for plotting hold on v=[20 70 100 150]; plot(xi,yi) axis(v) xlabel('Age (years)') ylabel('Target Heart Rate')
x = 20 25 30 35 40 45 50 55 60 65 70 y = 150 146 142 139 135 131 127 124 120 116 112
a1 = -0.7527 164.
hold on plot(h,p,'') % Plot p vs h table data in left side plot(hi,pi) % Plot quadratic fit hi vs pi in left side xlabel('Altitude, h (ft x 1000)') ylabel ('Air Pressure, p (psi)') title('Air Pressure Vs. Altitude') subplot(1,2,2) % Make right side active hold on plot(h,v,'') % Plot v vs h table data in right side plot(hi,vi) % Plot quadratic fit vi vs hi in right side xlabel('Altitude, h (ft x 1000)') ylabel ('Speed of Sound, v (ft/sec)') title('Speed of Sound Vs. Altitude')
0 1 0 2 0 3 0 4 0
2
4
6
8
1 0
1 2
1 4
1 6
A l t i t u d e , h ( f t x 1 0 0 0 )
Air Pressure, p (psi)
A i r P r e s s u r e V s. A l t i t u d e
0 1 0 2 0 3 0 4 0
9 4 0
9 6 0
9 8 0
1 0 0 0
1 0 2 0
1 0 4 0
1 0 6 0
1 0 8 0
1 1 0 0
1 1 2 0
A l t i t u d e , h ( f t x 1 0 0 0 )
Speed of Sound, v (ft/sec)
S p e e d o f S o u n d V s. A l t i t u d e
