Gram-Schmidt orthogonalization and approximation

We want to find the polynomial of degree 4 that is the best approximation to cos(x) on the interval [-Pi, Pi] .  To do this we first take the standard basis for polynomials, {1, x, x^2, x^3, x^4 } and use Gram-Schmidt to construct an orthonormal basis. While we could do this more efficiently using built-in commands of maple, we will just use maple as a calculator.

First, we define our inner product as a function on two functions.  Our inner product is the usual one, that is
`<,>`(p, q) = int(p(x)*q(x), x = -Pi .. Pi) .   Note that ||p-q|| is small precisely when the values of p(x) are close to those of q(x) for all x in the interval [-Pi, Pi] .

For simplicity, we'll operate on expressions in x rather than general functions.

> InnerProd:=(p,q) -> int(p*q, x=-Pi..Pi);

InnerProd := proc (p, q) options operator, arrow; int(p*q, x = -Pi .. Pi) end proc

We'll also want to compute the norm of our functions.  However, Norm is a maple built-in function, so we'll be a little more formal and call Norm by his full name.

> Norman:=p->sqrt(InnerProd(p,p));

Norman := proc (p) options operator, arrow; sqrt(InnerProd(p, p)) end proc

Now we do Gram-Schmidt to construct an orthonormal basis (with respect to this inner product) for polynomials of degree 4.

The first element of this basis is the constant function f[0](x) = 1 .  However, we need to normalize it to be of unit length in our norm.

> f[0] := 1/Norman(1);

f[0] := 1/2*2^(1/2)/Pi^(1/2)

Next, we confirm that f[1](x) = x is orthogonal to f[0] .  Since we want an orthonormal basis, we still want to divide f[1] by its length

> InnerProd(x,f[0]);

0

> f[1]:= x/Norman(x);

f[1] := 1/2*x*6^(1/2)/Pi^(3/2)

Notice that while  f[1](x) is orthogonal to x^2 , f[0] is not.

> InnerProd(x^2,f[0]);

1/3*2^(1/2)*Pi^(5/2)

> InnerProd(x^2,f[1]);

0

This means we have to use Gram-Schmidt to find a function in the span of {f[0], f[1], x^2} which is orthogonal to both f[0] and f[1] .  Note that since f[0] and f[1] are already normalized, we don't have to divide by their lengths.

> ftemp:= x^2 - InnerProd(x^2,f[0])*f[0] - InnerProd(x^2,f[1])*f[1];
f[2]  := collect(ftemp/Norman(ftemp),x);

ftemp := x^2-1/3*Pi^2

f[2] := 3/4*10^(1/2)*x^2/Pi^(5/2)-1/4*10^(1/2)/Pi^(1/2)

Now we play the same game with x^3 and x^4 .  Note that if you were working these integrals out by hand, you'd quickly realize that the even powers are orthogonal to the odd ones.  

> ftemp:= x^3 - InnerProd(x^3,f[0])*f[0] - InnerProd(x^3,f[1])*f[1] - InnerProd(x^3,f[2])*f[2]:
f[3]  := collect(ftemp/Norman(ftemp),x);

f[3] := 5/4*14^(1/2)*x^3/Pi^(7/2)-3/4*14^(1/2)*x/Pi^(3/2)

> ftemp:= x^4 - InnerProd(x^4,f[0])*f[0] - InnerProd(x^4,f[1])*f[1] - InnerProd(x^4,f[2])*f[2] - InnerProd(x^4,f[3])*f[3]:
f[4]  := collect(ftemp/Norman(ftemp),x);

f[4] := 105/16*2^(1/2)*x^4/Pi^(9/2)-45/8*2^(1/2)*x^2/Pi^(5/2)+9/16*2^(1/2)/Pi^(1/2)

So, we now have an orthonormal basis for fourth degree polynomials on [-Pi, Pi] .  Just in case you are disturbed by all the Pi s and square roots, let's also write their approximate values, and also draw their graphs.

> OrthoBasis:=[f[0],f[1],f[2],f[3],f[4]];

OrthoBasis := [1/2*2^(1/2)/Pi^(1/2), 1/2*x*6^(1/2)/Pi^(3/2), 3/4*10^(1/2)*x^2/Pi^(5/2)-1/4*10^(1/2)/Pi^(1/2), 5/4*14^(1/2)*x^3/Pi^(7/2)-3/4*14^(1/2)*x/Pi^(3/2), 105/16*2^(1/2)*x^4/Pi^(9/2)-45/8*2^(1/2...

> evalf(OrthoBasis,4);

[.3988, .2197*x, .1354*x^2-.4458, 0.8494e-1*x^3-.5035*x, 0.5359e-1*x^4-.4540*x^2+.4486]

> plot(OrthoBasis,x=-Pi..Pi, color=[black,red,blue,green,brown],thickness=2);

[Plot]

The orthonormal basis we have constructed is closely related to the Legendre polynomials, one of several families of orthogonal polynomials.

Now we can turn to our minimization problem, which becomes easy with an orthonormal basis.

Recall that the minimum with respect to our norm is just the projection of the target function onto the subspace spanned by the basis.  With an orthogonal basis, this is trivial since it is just a sum of inner products times each basis element.
While we know that the inner product of
cos(x) with f[0] , f[1] , and f[3] is zero, we write them anyway.

> approx := collect(
 InnerProd(cos(x),f[0])*f[0] + InnerProd(cos(x),f[1])*f[1] + InnerProd(cos(x),f[2])*f[2] + InnerProd(cos(x),f[3])*f[3] + InnerProd(cos(x),f[4])*f[4], x);

evalf(approx);

approx := -1575/8*(2*Pi^2-21)*x^4/Pi^8+(-45/2/Pi^4+675/4*(2*Pi^2-21)/Pi^6)*x^2+15/2/Pi^2-135/8*(2*Pi^2-21)/Pi^4

0.2615982020e-1*x^4-.4522878091*x^2+.9783263892

The fit is pretty good, differing by less than .07 on the whole interval.

> plot([cos(x),approx], x=-Pi..Pi, color=[red, green], thickness=2);

[Plot]

For comparison, let's look at the Taylor polynomial for the cosine on the same interval.  Note that while the Taylor polynomial fits much better for small values, over the whole interval it is much worse.

> tailor:=convert(taylor(cos(x),x,5),polynom);

> plot([cos(x), tailor], x=-Pi..Pi, color=[red, blue], thickness=2);

tailor := 1-1/2*x^2+1/24*x^4

[Plot]

More explictly, we can compare the norms of the differences, and see the dramatic difference.  We also show a plot of the two differences on the same axes.

> evalf(Norman(cos(x)-approx));
evalf(Norman(cos(x)-tailor));

0.4601613986e-1

.7995084396

> plot([cos(x)-approx, cos(x)-tailor], x=-Pi..Pi, color=[green,blue], thickness=2);

[Plot]

Just to emphasize the ease of solving the minimization problem once the basis is in hand, here are a couple more examples.

Naturally, since our basis was constructed specifically for the interval [-Pi, Pi] , we must work there. However, it would be not hard to adapt this same basis to work on any interval.

In order to avoid lots of typing, we write a Maple function that computes the inner product for us, and another that shows us the graphs.

> MakeApprox:=(f,basis)->collect(add(InnerProd(f,basis[i])*basis[i], i=1..nops(basis)),x):

ShowApprox:=proc(f,basis)

 local a;

 a:=MakeApprox(f,basis);

 print(a);

 print(evalf(a), " norm of difference=",evalf(Norman(a-f)));

 plot([f, a], x=-Pi..Pi, color=[red, green], thickness=2);

end:

First, just to check, we apply it to the same problem we just solved.

> ShowApprox(cos(x),OrthoBasis);

-1575/8*(2*Pi^2-21)*x^4/Pi^8+(-45/2/Pi^4+675/4*(2*Pi^2-21)/Pi^6)*x^2+15/2/Pi^2-135/8*(2*Pi^2-21)/Pi^4

0.2615982020e-1*x^4-.4522878091*x^2+.9783263892,

[Plot]

Now we will try several other functions, just to show off.  First, sin(x), which doesn't do as well (we'd need to add a fifth degree polynomial to our basis for a good fit.)

> ShowApprox(sin(x),OrthoBasis);

35/2*(-15+Pi^2)*x^3/Pi^6+(3/Pi^2-21/2*(-15+Pi^2)/Pi^4)*x

-0.9338769716e-1*x^3+.8569833271*x,

[Plot]

> ShowApprox(cos(1+x)-x,OrthoBasis);

-1575/8*cos(1)*(2*Pi^2-21)*x^4/Pi^8-35/2*sin(1)*(-15+Pi^2)*x^3/Pi^6+(-45/2*cos(1)/Pi^4+675/4*cos(1)*(2*Pi^2-21)/Pi^6)*x^2+(-(3*sin(1)+Pi^2)/Pi^2+21/2*sin(1)*(-15+Pi^2)/Pi^4)*x+15/2*cos(1)/Pi^2-135/8*c...-1575/8*cos(1)*(2*Pi^2-21)*x^4/Pi^8-35/2*sin(1)*(-15+Pi^2)*x^3/Pi^6+(-45/2*cos(1)/Pi^4+675/4*cos(1)*(2*Pi^2-21)/Pi^6)*x^2+(-(3*sin(1)+Pi^2)/Pi^2+21/2*sin(1)*(-15+Pi^2)/Pi^4)*x+15/2*cos(1)/Pi^2-135/8*c...

0.1413421117e-1*x^4+0.7858303749e-1*x^3-.2443721461*x^2-1.721126604*x+.5285920041,

[Plot]

> ShowApprox(exp(x)*cos(x),OrthoBasis);

-315/128*(90*Pi^2+210*Pi-4*Pi^4+105-90*exp(2*Pi)*Pi^2+210*exp(2*Pi)*Pi+4*exp(2*Pi)*Pi^4-105*exp(2*Pi))*exp(-Pi)*x^4/Pi^9-35/16*(-15-15*Pi+2*Pi^3+15*exp(2*Pi)-15*exp(2*Pi)*Pi+2*exp(2*Pi)*Pi^3)*exp(-Pi)...-315/128*(90*Pi^2+210*Pi-4*Pi^4+105-90*exp(2*Pi)*Pi^2+210*exp(2*Pi)*Pi+4*exp(2*Pi)*Pi^4-105*exp(2*Pi))*exp(-Pi)*x^4/Pi^9-35/16*(-15-15*Pi+2*Pi^3+15*exp(2*Pi)-15*exp(2*Pi)*Pi+2*exp(2*Pi)*Pi^3)*exp(-Pi)...-315/128*(90*Pi^2+210*Pi-4*Pi^4+105-90*exp(2*Pi)*Pi^2+210*exp(2*Pi)*Pi+4*exp(2*Pi)*Pi^4-105*exp(2*Pi))*exp(-Pi)*x^4/Pi^9-35/16*(-15-15*Pi+2*Pi^3+15*exp(2*Pi)-15*exp(2*Pi)*Pi+2*exp(2*Pi)*Pi^3)*exp(-Pi)...-315/128*(90*Pi^2+210*Pi-4*Pi^4+105-90*exp(2*Pi)*Pi^2+210*exp(2*Pi)*Pi+4*exp(2*Pi)*Pi^4-105*exp(2*Pi))*exp(-Pi)*x^4/Pi^9-35/16*(-15-15*Pi+2*Pi^3+15*exp(2*Pi)-15*exp(2*Pi)*Pi+2*exp(2*Pi)*Pi^3)*exp(-Pi)...-315/128*(90*Pi^2+210*Pi-4*Pi^4+105-90*exp(2*Pi)*Pi^2+210*exp(2*Pi)*Pi+4*exp(2*Pi)*Pi^4-105*exp(2*Pi))*exp(-Pi)*x^4/Pi^9-35/16*(-15-15*Pi+2*Pi^3+15*exp(2*Pi)-15*exp(2*Pi)*Pi+2*exp(2*Pi)*Pi^3)*exp(-Pi)...-315/128*(90*Pi^2+210*Pi-4*Pi^4+105-90*exp(2*Pi)*Pi^2+210*exp(2*Pi)*Pi+4*exp(2*Pi)*Pi^4-105*exp(2*Pi))*exp(-Pi)*x^4/Pi^9-35/16*(-15-15*Pi+2*Pi^3+15*exp(2*Pi)-15*exp(2*Pi)*Pi+2*exp(2*Pi)*Pi^3)*exp(-Pi)...

-.1116938107*x^4-.5009320068*x^3-.2395710070*x^2+1.204634804*x+1.126116582,

[Plot]



Scott Sutherland May 9, 2005