CHAPTER 3: NUMERICAL ALGORITHMS WITH MATLAB

Lecture 3.1: Polynomials and polynomial interpolation

General properties of polynomials:

·         A general polynomial of order n:

y = Pn(x) = c1 xn + c2 xn-1 + … + cn-1 x2 + cn x + cn+1

·         Factorization of polynomials of n-th order by n roots (zeros) x1,x2,…,xn (possibly, complex or multiple):

y = Pn(x) = c1 (x – x1)*(x – x2)*…*(x – xn)

·         Horner's method for nested multiplications of polynomials:

y = Pn(x) = ( … ( ( ( c1*x + c2 )*x + c3  )*x + c4) … )*x + cn+1

Example:

y = x4 + 2 x3 – 7 x2 – 8 x + 12

y = (x-1)*(x-2)*(x+2)*(x+3)

y = ( ( (x + 2)*x – 7)*x – 8)*x + 12

Polynomials with MATLAB:

·         representation of polynomials by a row-vector of coefficients

p1 =     1     2    -7    -8    12

p2 =     2     3     5     9     5

·         roots: find all zeros of a polynomial from its coefficients

r1 =   -3.0000   -2.0000    2.0000    1.0000

r2 =   0.2500 - 1.5612i   0.2500 + 1.5612i  -1.0000            -1.0000

·         poly: find coefficients of a polynomial from its zeros

q1 = poly(r1), q2 = poly(r2)

% NB: coefficients q2 are different from coefficients of p2 by a factor of 2

% The leading-order coefficient c1 is always normalized to 1.

q1 =    1.0000    2.0000   -7.0000   -8.0000   12.0000

q2 =    1.0000    1.5000    2.5000    4.5000    2.5000

% rounding errors may occur in computations of roots

% Wilkinson's example:

roots(poly(1:20))'

ans =

Columns 1 through 11

20.0003   18.9970   18.0118   16.9695   16.0509   14.9319   14.0684   12.9472   12.0345   10.9836   10.0063

Columns 12 through 20

8.9983    8.0003    7.0000    6.0000    5.0000    4.0000    3.0000    2.0000    1.0000

% large rounding error occurs when highly multiple roots are computed

% y = (x-1)^6

r = [ 1 1 1 1 1 1 ]; p = poly(r), rr = roots(p)'

p =     1    -6    15   -20    15    -6     1

rr =   1.0042 - 0.0025i   1.0042 + 0.0025i   1.0000 - 0.0049i   1.0000 + 0.0049i   0.9958 - 0.0024i   0.9958 + 0.0024i

·         polyval: evaluate polynomials at a given point

y1 = polyval(p1,2.5)

x = [ 0 : 0.2 : 1]; y = polyval(p2,x)

y1 =   18.5625

y =    5.0000    7.0272    9.6432   13.1072   17.7552   24.0000

% implementation of the Horner algorithm for nested multiplication:

function [y] = HornerMultiplication(p,x)

[n,m] = size(x);

y = p(1)*ones(n,m);

for k = 2:length(p)

y = y.*x + p(k);

end

y1 = HornerMultiplication(p1,2.5)

y = HornerMultiplication(p2,x)

y1 =

18.5625

y =

5.0000    7.0272    9.6432   13.1072   17.7552   24.0000

·         polyder: computes coefficients of the derivative of a given polynomial

% P(x) = c(1) x^n + c(2) x^(n-1) + … + c(n) x + c(n+1)

% P'(x) = n c(1) x^(n-1) + (n-1) c(2) x^(n-2) + … + c(n)

Pder1 = polyder(p1), Pder2 = polyder(p2)

Pder1 =     4     6   -14    -8

Pder2 =     8     9    10     9

·         polyint:  computes coefficients of the integral of a given polynomial

% P(x) = c(1) x^n + c(2) x^(n-1) + … + c(n) x + c(n+1)

% c(n+2) are constant of integration (to be defined)

Pint1 = polyint(p1,10) % the constant of integration is 10

Pint2 = polyint(p2)  % the constant of integration is 0 (default)

Pint1 =    0.2000    0.5000   -2.3333   -4.0000   12.0000   10.0000

Pint2 =    0.4000    0.7500    1.6667    4.5000    5.0000         0

·         conv: computes a product of two polynomials

% Let p1 be polynomial of order n, p2 be polynomial of order m,

% conv(p1,p2) is polynomial of order (n+m)

p = conv(p1,p2)

·         deconv: computes a division of two polynomials

% Let p1 be polynomial of order n, p2 be polynomial of order m,

% [pq,pr] = deconv(p1,p2) computes the quotient polynomial pq of order (n-m)

% and the remainder polynomial pr of order (m-1) such that p1 = p2*pq + pr

[pq,pr] = deconv(p1,p2)

pq =    0.5000

pr =         0    0.5000   -9.5000  -12.5000    9.5000

·         residue:  computes partial-fraction expansion (residues)

% Let p1 be polynomial of order n and p2 be polynomial of order m

% [C,X,R] = residue(p1,p2) finds coefficients C of the residue terms,

% locations of poles X and the remainder term of a partial fraction expansion

% of the ratio of two polynomials p1(x)/p2(x).

% Example: no multiple roots,

%    p1(x)       C(1)       C(2)             C(n)

%    ----  =  -------- + -------- + ... + -------- + R(x)

%    p2(x)     x - X(1)   x - X(2)         x - X(n)

[C,X,R] = residue(p2,p1)

C =   -5.2000

1.2500

4.9500

-2.0000

X =   -3.0000

-2.0000

2.0000

1.0000

R =     2

% the inverse operation:

% from partial fraction expansion to the ratio of two polynomials

[q1,q2] = residue(C,X,R)

q1 =    2.0000    3.0000    5.0000    9.0000    5.0000

q2 =    1.0000    2.0000   -7.0000   -8.0000   12.0000

·         polyfit: computes coefficients of interpolating polynomials of order n that passes through a set of (n+1) data points

x = [ 1,2,-2,-3,0]; y = [0,0,0,0,12];

p = polyfit(x,y,4) % 4 is the order of the polynomial through 5 data points

pp = poly([1,2,-2,-3]) % the same operation if roots of polynomial are known

p =    1.0000    2.0000   -7.0000   -8.0000   12.0000

pp =     1     2    -7    -8    12

Polynomial interpolation:

Problem: Given a set of (n+1) data points:

(x1,y1); (x2,y2); …; (xn,yn); (xn+1,yn+1)

Find a polynomial of degree n, y = Pn(x),  that passes through all (n+1) data points. Assume that the set x1,x2,…,xn,xn+1 is ordered in ascendental order, i.e. x1<x2<…<xn<xn+1. The polynomial y = Pn(x) is called the interpolating polynomial for x1<x<xn+1 and extrapolating polynomial for x < x1 and x > xn+1.

Example: Given a set of five data points for a voltage-current characteristic of a zener diode

 Voltage -1 -0.75 -0.5 -0.25 -0 Current -14.58 -6.15 -1.82 -0.23 0

The five points are connected by a polynomial y = P4(x) of order 4. The data points are shown by blue stars, the polynomial is shown by green solid line:

Methods:

·         Power series (Vandermonde interpolation)

·         Lagrange interpolation polynomials

·         Newton forward difference interpolation

·         Newton backward difference interpolation

Vandermonde interpolation:

The (n+1) coefficients ck of the interpolating polynomial y = Pn(x) are to be matched with (n+1) data points:

Pn(xk) = yk.

The matching conditions produce a system of  (n+1) linear equations for (n+1) unknown coefficients of the interpolating polynomial:

c1 (xk)n +  c2 (xk)n-1 + … + cn-1 (xk)2 + cn xk + cn+1 = yk

The linear system can be solved with the use of MATLAB linear algebra solver.

x = [ -1,-0.75,-0.5,-0.25,-0]; A = vander(x)

y = [ -14.58;-6.15;-1.82;-0.23;-0.00]; c = A\y; c'

xInt = -1 : 0.01 : 0; yInt = HornerMultiplication(c,xInt);

plot(xInt,yInt,'g',x,y,'b*');

A =

1.0000   -1.0000    1.0000   -1.0000    1.0000

0.3164   -0.4219    0.5625   -0.7500    1.0000

0.0625   -0.1250    0.2500   -0.5000    1.0000

0.0039   -0.0156    0.0625   -0.2500    1.0000

0         0         0         0    1.0000

ans =

0.2133   15.0400    0.3067    0.0600         0

The Vandermond matrices are ill-conditioned for large n and the linear algerba solvers produce inaccurate numerical results. Lagrange and Newton's interpolations are direct methods of finding the interpolating polynomials, without the necessity to solve linear systems of equations.

Lagrange interpolation:

Lagrange polynomials are defined in the factorization form:

Ln,j(x) =

·         Ln,j(x) is a polynomial of order n

·         Ln,j(xi) = 0 for any i j

·         Ln,j(xj) = 1

The interpolating polynomial y = Pn(x) has a simple form in terms of the Lagrange polynomials Ln,j(x):

Pn(x) = = y1Ln,1(x) + … +yn+1Ln,n+1(x)

The interpolating polynomial y = Pn(x) has order n and it passes through all (n+1) data points. It can be proved that interpolating polynomials are unique, i.e. the polynomials y = Pn(x) obtained by the Vandermonde method and by the Lagrange method are the same polynomials.

function [yi] = LagrangeInter(x,y,xi)

% Lagrange interpolation algorithm

% x,y - row-vectors of (n+1) data values (x,y)

% xi - a row-vector of values of x, where the polynomial y = Pn(x) is evaluated

% yi - a row-vector of values of y, evaluated with y = Pn(x)

n = length(x) - 1; % order of interpolation polynomial y = Pn(x)

ni = length(xi); % number of points where the interpolation is to be evaluated

L = ones(ni,n+1); % the matrix for Lagrange interpolating polynomials L_(n,j)(x)

% L has (n+1) columns for each point j = 1,2,...,n+1

% L has ni rows for each point of xi

for j = 1 : (n+1)

for i = 1 : (n+1)

if (i ~= j)

L(:,j) = L(:,j).*(xi' - x(i))/(x(j)-x(i));

end

end

end

yi = y*L';

Example of Lagrange interpolation:

x = [ -1,-0.75,-0.5,-0.25,-0];

y = [ -14.58,-6.15,-1.82,-0.23,-0.00];

xInt = -1 : 0.01 : 0;

yInt = LagrangeInter(x,y,xInt);

plot(xInt,yInt,'g',x,y,'b*');