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]

 

p1 = 1 2 -7 -8 12

p2 = 2 3 5 9 5

 

         roots: find all zeros of a polynomial from its coefficients

 

r1 = roots(p1)', r2 = roots(p2)'

 

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)

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

% 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)

 

p = 2 7 -3 -18 -12 -57 -47 68 60

 

         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.00

-0.75

-0.50

-0.25

-0.00

Current

-14.58

-6.15

-1.82

-0.23

0.00

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*');