Fall EE103 2007 Lecture Notes (SEJ) Section 8 SECTION 8: INTRODUCTION TO ORTHOGONAL POLYNOMIALS............. 122 Orthogonal Polynomials........................................................................................... 124 Trigonometric Polynomials...................................................................................... 125 Background: .......................................................................................................... 125 Trigonometric Interpolation ................................................................................ 127 Examples of Trigonometric Interpolation:......................................................... 129
EE103 Fall 2007 Lecture Notes (SEJ) Section 8 SECTION 8: INTRODUCTION TO ORTHOGONAL POLYNOMIALS Ex: Consider the following data for which we wish to compute a least-squares approximating polynomial.
[x y] 10.0000 10.2000 10.4000 10.6000 10.8000 11.0000
0 0.0040 0.0160 0.0360 0.0640 0.1000
The X matrix, for the approximating quadratic, is
X= 100.0000 104.0400 108.1600 112.3600 116.6400 121.0000 10.0000 10.2000 10.4000 10.6000 10.8000 11.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
The condition number of X T X is given by MatLab to be cond(X'*X) = 1.5598e+010. Keep in mind that the "basis polynomials" that we've been using are, for a quadratic leastsquares approximating polynomial,
2 ( x) = x2 1 ( x ) = x 0 ( x) 1
However, suppose we choose other "basis polynomials"; e.g.,
2 ( x ) = ( x - 10.5) 2 - 0.1166667 1 ( x ) = x - 10.5 0 ( x) 1
Do not, at this point, be concerned with the derivation of these polynomials; however, as before, the X matrix becomes
122
EE103 Fall 2007 Lecture Notes (SEJ) Section 8
2 ( x1 ) 1 ( x1 ) ( x ) ( x ) 2 2 1 2 X = 2 ( x6 ) 1 ( x6 ) 1
1 1
and we wish to compute c2 , c1 , c0 so that P2 ( x ) = c2 2 ( x ) + c11 ( x ) + c0 0 ( x ) is the least-squares approximating polynomial. In this case, the X matrix is 0.1333 -0.0267 -0.1067 -0.1067 -0.0267 0.1333 -0.5000 -0.3000 -0.1000 0.1000 0.3000 0.5000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
and the condition number is cond(X'*X) = 100.4464, a remarkable reduction. However, this is not all: The matrix X'X = 0.0597 0 0.0000 0 0.7000 0 0.0000 0 6.0000
a diagonal matrix! Therefore, the system of equations, X T Xc = X T y can be trivially solved. Of course, we need to see why this has occurred. The functions x m , x m -1 , , x, 1
comprise a set of "basis polynomials" for the vector space of m th degree polynomials. That is, any m th degree polynomial is a linear combination of these "basis polynomials"; e.g.,
am x m + am-1 x m -1 +
a1 x + a0
123
EE103 Fall 2007 Lecture Notes (SEJ) Section 8 Orthogonal Polynomials Suppose we can find, or construct, another set of "basis polynomials" for the vector space of m th degree polynomials; e.g.,
m ( x ), m -1 ( x ),
and, suppose further, that
, 1 ( x ), 0 ( x )
( x ) ( x ) = 0
k =1 i k j k
n
for all i , j with i j
Definition: When this latter condition holds, we say the "basis polynomials" are orthogonal. [Note: The definition depends upon the x values.]
HW: Show that the "basis polynomials" x 2 , x , 1 are not orthogonal with respect to the data of the above example; then show that the functions
2 ( x ) = ( x - 10.5) 2 - 0.1166667 1 ( x ) = x - 10.5 0 ( x) 1
are orthogonal. Recall the following about "orthogonal" vectors. The vectors { v1 , v 2 , orthogonal if for all i , j (i j ) v i v j = 0 Therefore, the matrix V TV is a diagonal matrix. To see this, consider the (i , j )th , i j , element of the matrix V TV . That element is the inner product of the i th row of V T with the j th column of V . But the i th row of V T is the i th column of V . But we've assumed orthogonality and, therefore, the (i , j )th , i j , element is zero. If the "basis functions" m ( x ), m -1 ( x ), data set ( xk , yk ), k = 1, , 1 ( x ), 0 ( x ) are orthogonal (with respect to the
T
, v k } are said to be
, n ) then the X T X matrix is a diagonal matrix, where
124
EE103 Fall 2007 Lecture Notes (SEJ) Section 8
m ( x1 ) m ( x2 ) X = m ( x3 ) ( x ) m n
m -1 ( x1 ) m-2 ( x1 ) 1 ( x1 ) 0 ( x1 ) m -1 ( x2 ) m-2 ( x2 ) 1 ( x2 ) 0 ( x2 ) m -1 ( x3 ) m -2 ( x3 ) 1 ( x3 ) 0 ( x3 ) m-1 ( xn ) m -2 ( xn )
( X T X )c = X T y 1 ( xn ) 0 ( x n )
and, hence, the solution for the coefficients, via,
is simply computed. To see this, note that the definition of "orthogonal basis functions" is
( x )
k =1 i k
n
j
( xk ) = 0 for all i , j with i j
and, further note, this is precisely the "inner product" of the i th and j th columns of X . There is a computational implementation of the Gram-Schmidt process for computing a set of orthogonal "basis polynomials", given the data ( xk , yk ), k = 1, However, we will not develop that method in these notes.
, n for the problem.
Rather, we introduce
approximation/interpolation based upon "trigonometric polynomials", indeed, orthogonal
trigonometric polynomials. This topic is at the heart of Fourier analysis and engineering
applications.
Trigonometric Polynomials
Perhaps the most important engineering applications of polynomial interpolation and approximation involve the use of trigonometric polynomials.
Background:
Let f be a function of period 2 ; i.e., f ( x + 2 ) = f ( x ) , x . Fourier's great insight is that virtually any such function can be expressed as
f ( x) =
a0 2
+ [a j cos( jx ) + bk sin( jx )]
j =0
(1)
Of course, cos( jx ), sin( jx ) , both have period 2 . If one is concerned with a function f of period T, then (1) can be rewritten as
125
EE103 Fall 2007 Lecture Notes (SEJ) Section 8
f ( x) =
a0 2
+ [a j cos( jw0 x ) + b j sin( jw0 x )]
j =0
(2)
where w0 = 2 / T , the so-called "fundamental frequency".
Ex: Show that cos( jw0 x ) and sin( jw0 x ) have periods equal to T .
The leading coefficient,
a0 2
, is generally included in various discussions merely to enable
closed-form solutions for the a j , b j coefficients. Those solutions are:
1 aj =
-
f ( x ) cos( jx )dx , j = 0,1, 2, j = 0,1, 2, (3)
bj =
1
f ( x)sin( jx)dx ,
-
These solutions are the consequence of the fact that the "trigonometric basis functions" are orthogonal, in the following sense:
cos( kx ) sin( rx )dx = 0
-
cos( kx ) cos( rx )dx = 0 , k r
-
sin( kx ) sin( rx )dx = 0 , k r
-
HW: Show that typical examples of the above are true; you may use Matlab's symbolic toolbox, if you wish. Note: One could also choose the interval of integration to be
[0, 2 ] , or
[0, T ] or [ - T , T ] when using the fundamental frequency w0 = 2 / T . 2 2
The solutions given in (3) can be derived, although a good deal of somewhat technical mathematics is involved (of course, deriving the values of the coefficients depends upon the creativity of Fourier: i.e., equation (1) ). The simplest derivation that is for a0 . Integrating both sides of equation (1) yields
-
a f ( x )dx = 0 + [a j cos( jx ) + b j sin( jx )] dx - 2 j =0
=
a0 2
-
dx + a j cos( jx )dx + b j sin( jx )dx
j =0 - j =0 -
126
EE103 Fall 2007 Lecture Notes (SEJ) Section 8 assuming the interchanges of integration and summation are valid. summations are equal to zero and, therefore, we have
a 0 = -1
-
The last two
f ( x )dx
Note: We'd have exactly the same result if the interval of integration is [0, 2 ] .
The usefulness of Equations (1) or (2) arises by utilizing a modification when given a finite sample for f (e.g., we've sampled the generally unknown function f , assumed to be of period T )
T xk = k N , k = 0,1, 2,
,N
f ( xk ) , f ( x0 ) = f ( xN ) on the interval [0, T ] . Equation (2) suggests we use a finite sum
f ( x ) Fm ( x )
a0 2
+ [a j cos( jw0 x ) + b j sin( jw0 x )]
j =0
m
(4)
where m should be chosen, relative to N , so that 2m + 1 N (note: b0 doesn't appear since sin(0) = 0 ; therefore m should be chosen so that N is large enough to provide enough information to determine the 2m + 1 coefficients a0 , a1 ,
Trigonometric Interpolation1
, am ; b1 , b2 ,
, bm ).
Rather than use (4), we "equivalently" use in its place, dropping the leading term a0 / 2 ,
f ( x ) Fm ( x )
[a
j =0
m
j
cos( jw0 x ) + b j sin( jw0 x )]
(5)
Assume f to be a T periodic function; we have "sampled" the function at N + 1 points, equally spaced on the interval [0, T ] , and we have the following data at hand:
( x0 , f 0 ), ( xN , f N ) , f N = f 0 ,N
xk = kT / N , k = 0, 1,
1
This discussion follows that of Charles Van Loan, Introduction to Scientific Computing.
127
EE103 Fall 2007 Lecture Notes (SEJ) Section 8 where, of course, f k = f ( xk ) . Assume that N is even and that m = N / 2 . We seek coefficients a0 , points ( x0 , f 0 ),
, am , b0 , , bm so that the function in (5), Fm ( x ) , interpolates the data , am , b0 , , bm so that
, ( x N , f N ) . That is, we seek coefficients a0 ,
f ( xk ) = Fm ( xk ) , k = 0, , N (6) As noted above, b0 does not actually appear in Fm ( x ) because sin( 0) = 0 . It's also the
case that bm does not appear because sin( jw0 x N ) = sin( jw0T ) = sin( 2 j ) = 0 . Also, the first and last equations of (6) are identical, due to periodicity of f , i.e., f n = f 0 . Therefore, we delete the last equation.
We now have an NxN system of linear equations to determine the coefficients of the
interpolating trigonometric polynomial:
f k = a0 + [a j cos( jw0 xk ) + b j sin( jw0 xk )] + am cos( mw0 xk ) , k = 0,
j =1 m -1
, N -1
(7)
and we see that we seek m + 1 " a " coefficients and m - 1 " b " coefficients. The a ' s are associated with the cosine function and the b ' s are associated with the sine function. Since w0 = 2 / T , m = N / 2 , and each xk = kT / N , we have that each jw0 xk = jk / m . Hence, (7) becomes
f k = a0 + [a j cos( jk / m ) + b j sin( jk / m )] + am cos( k ) , k = 0,
j =1 m -1
, N -1
(8)
128
EE103 Fall 2007 Lecture Notes (SEJ) Section 8 The following, as an example, is a resulting 6x6 matrix ( N = 6, m = 3 ):
1 1 1 1 1 1
1 1 1 1 1 1
cos 0 cos 3 cos 23 cos 33 cos 43 cos 53
1.0000 0.5000 -0.5000 -1.0000 -0.5000 0.5000
cos 0 cos 23 cos 43 cos 63 cos 83 cos 10 3
1.0000 -0.5000 -0.5000 1.0000 -0.5000 -0.5000
1 -1 1 -1 1 -1
1 -1 1 -1 1 -1
sin 0 sin 3 sin 23 sin 33 sin 43 sin 53
0 0.8660 0.8660 0.0000 -0.8660 -0.8660
sin 0 a0 sin 23 a1 sin 43 a2 = sin 63 a3 sin 83 b1 sin 10 b2 3
0 0.8660 -0.8660 -0.0000 0.8660 -0.8660
f0 f1 f2 f3 f4 f5
Examples of Trigonometric Interpolation:
Figure 1: function with a jump discontinuity at 1
129
EE103 Fall 2007 Lecture Notes (SEJ) Section 8
Figure 2: 12 point trigonometric polynomial interpolation
Figure 3: 24 point trigonometric polynomial interpolation
130
EE103 Fall 2007 Lecture Notes (SEJ) Section 8
Figure 4: 48 point trigonometric polynomial interpolation
Figure 5: Bell shaped function
131
EE103 Fall 2007 Lecture Notes (SEJ) Section 8
Figure 6: 8 point trigonometric polynomial interpolation
Figure 7: 16 point trigonometric polynomial interpolation
132
EE103 Fall 2007 Lecture Notes (SEJ) Section 8
The following MatLab function, due to Charles Van Loan, produces the solution for the interpolation equations (8): the coefficients, a j ' s for the cosine functions, and the b j ' s
for the sine functions, and, of course, the corresponding orthogonal matrix.
function [F,P] = CSInterp(f) % F = CSInterp(f): Adopted from Charles Van Loan % f is a column n vector and n = 2m. % F.a is a column m+1 vector and F.b is a column m-1 % vector so that if % tau = (0:n-1)'*pi/m, then % % f = F.a(1)*cos(0*tau) +...+ F.a(m+1)*cos(m*tau) + % F.b(1)*sin(tau) +...+ F.b(m-1)*sin((m-1)*tau) n = length(f); m = n/2; tau = (pi/m)*(0:n-1)'; P = []; for j=0:m, P = [P cos(j*tau)]; end for j=1:m-1, P = [P sin(j*tau)]; end y = P\f; %note: this code makes no use of P'P, diagonal. F = struct('a',y(1:m+1),'b',y(m+2:n)); %F.a are the weights on the cosines; F.b for the sines.
Chart 1: Matlab Function for the solution of the Trigonometric Interpolation Eqs
Of course, we're generally interested in evaluating the interpolating trigonometric polynomial at points other than x0 , x1 , interpolating polynomial.
, xN , and especially if we desire to plot the
The following is a succinct MatLab function, also due to Charles Van Loan, that produces the values of the trigonometric interpolating polynomial at each of the points in tvals (a selection of domain values chosen by the user); the other two inputs are (1) the output "structure" F from CSInterp.m, and (2) the period T.
133
EE103 Fall 2007 Lecture Notes (SEJ) Section 8
function Fvals = CSeval(F,T,tvals) % Due to Charles Van Loan % F.a is a length m+1 column vector, F.b is a length m-1 column vector, % T is a positive scalar, and tvals is a column vector. % If % F(t) = F.a(1) + F.a(2)*cos((2*pi/T)*t) +...+ F.a(m+1)*cos((2*m*pi/T)*t) + % F.b(1)*sin((2*pi/T)*t) +...+ F.b(m-1)*sin((2*m*pi/T)*t) % % then Fvals = F(tvals). Fvals = zeros(length(tvals),1); tau = (2*pi/T)*tvals; for j=0:length(F.a)-1, Fvals = Fvals + F.a(j+1)*cos(j*tau); end for j=1:length(F.b), Fvals = Fvals + F.b(j)*sin(j*tau); end
Chart 2: Matlab Function for Evaluating the Trigonometric Interpolation Polynomial at a selection of domain values, tvals.
%clear all;clc; n=16; x=linspace(-1,1,n+1)'; f=1./(1+25*x.^2); [F,P]=CSinterp(f(1:n)); xd=linspace(0,2,100)'; Fvals=CSeval(F,2,xd); plot(xd-1,Fvals,'r',x(1:n),f(1:n),'bd'); hold on; f1=inline('1./(1+25*x.^2)'); fplot(f1,[-1,1]); hold on; label=strcat({'familiar (?) function 16 sample points ,'},{' f=1./(1+25*x.^2) '}); title(label);grid on; hold off;
Chart 3 (Tutorial): Matlab Script for Solving for a,b coefficients, and Plotting the sampled points, the actual sampled function, and the Trig Interpolating Polynomial (note selection of xd and xd-1 (why?))
134