CHAPTER 4 PROJECT 1

Function Spaces-Orthogonal Polynomials and Fourier Series

===============================================================

HINTS/ANSWERS TO PROBLEMS

================================================================

Project Part 1: Find the first 4 orthogonal polynomials as defined. Use them to approximate exp(x) from -1 to 1. Graph the result of the Legendre approximation with the function y = exp(x). Since the polynomials were not created here to be orthonormal, you will have to divide by their length (using the inner product given) as you construct the function. On the same graph, include the graph created by the use of a Taylor series. You can look up some of these routines with the command "with(orthopoly)".

We will create an orthogonal, not orthonormal set.

> restart:f1:=1;

f1 := 1

>

> f2:=x-int(x*f1,x=-1..1)/int(f1*f1,x=-1..1)*f1;

f2 := x

>

> f3:=x^2-int(x^2*f1,x=-1..1)/int(f1*f1,x=-1..1)*f1-int(x^2*f2,x=-1..1)/int(f2*f2,x=-1..1)*f2;

f3 := x^2-1/3

>

Now we must create an orthogonal polynomial from x^3. We subtract its projections on 1, x, x^2.

> f4:=x^3- int(f1*x^3,x=-1..1)/int(f1*f1,x=-1..1)*f1 - int(f2*x^3,x=-1..1)/int(f2*f2,x=-1..1)*f2 - int(f3*x^3,x=-1..1)/int(f3*f3,x=-1..1)*f3;

f4 := x^3-3/5*x

>

Now we approximate exp(x)

> f:=int(exp(x)*f1 ,x=-1..1)/int(f1*f1,x=-1..1)*f1 + int(exp(x)*f2,x=-1..1)/int(f2*f2,x=-1..1)*f2+int(exp(x)*f3,x=-1..1)/int(f3*f3,x=-1..1)*f3 +int(exp(x)*f4,x=-1..1)/int(f4*f4,x=-1..1)*f4;

> evalf(f);

f := 1/2*exp(1)-1/2*exp(-1)+3*exp(-1)*x+45/8*(2/3*exp(1)-14/3*exp(-1))*(x^2-1/3)+175/8*(-2*exp(1)+74/5*exp(-1))*(x^3-3/5*x)

.9962940199+.9979548527*x+.5367215194*x^2+.1761391188*x^3

> with(plots):m1:=plot(f,x=-1..1,color=red):m2:=plot(exp(x),x=-1..1,color=blue):m3:=plot(1+x+x^2/2! + x^3/3!,x=-1..1,color=green):display(m1,m2,m3);

The Legendre approximation, which is "f", appears to be the better fit.

=========================================

Project Part 2:

Minimize the function below, by differentiating.

> F:=int((exp(x)-a-b*x-c*x^2-d*x^3)^2,x=-1..1);

F := -2*exp(1)*a+4*d*exp(1)-2*c*exp(1)-1/2*exp(-1)^2+2*exp(-1)*a-32*d*exp(-1)-4*b*exp(-1)+1/2*exp(1)^2+10*c*exp(-1)+2/3*b^2+4/5*b*d+2/5*c^2+2/7*d^2+4/3*c*a+2*a^2

> diff(F,a);

-2*exp(1)+2*exp(-1)+4/3*c+4*a

> diff(F,b);

-4*exp(-1)+4/3*b+4/5*d

> diff(F,c);

-2*exp(1)+10*exp(-1)+4/5*c+4/3*a

> diff(F,d);

4*exp(1)-32*exp(-1)+4/5*b+4/7*d

> M:=matrix( [[4,0,4/3, 0, 2*exp(1)-2*exp(-1)],[0,4/3, 0,4/5, 4*exp(-1)],[4/3, 0, 4/5, 0,2*exp(1)-10*exp(-1)], [0,4/5, 4/7,0, -4*exp(1)+32*exp(-1)]]);

M := matrix([[4, 0, 4/3, 0, 2*exp(1)-2*exp(-1)], [0, 4/3, 0, 4/5, 4*exp(-1)], [4/3, 0, 4/5, 0, 2*exp(1)-10*exp(-1)], [0, 4/5, 4/7, 0, -4*exp(1)+32*exp(-1)]])

> with(linalg):rref(M);

Warning, new definition for norm

Warning, new definition for trace

matrix([[1, 0, 0, 0, -3/4*exp(1)+33/4*exp(-1)], [0, 1, 0, 0, 235/4*exp(-1)-215/28*exp(1)], [0, 0, 1, 0, 15/4*exp(1)-105/4*exp(-1)], [0, 0, 0, 1, 1075/84*exp(1)-1115/12*exp(-1)]])

> f:=-3/4*exp(1)+33/4*exp(-1) + x *(235/4*exp(-1)-215/28*exp(1)) + x^2*(15/4*exp(1)-105/4*exp(-1))+x^3 *(1075/84*exp(1)-1115/12*exp(-1));evalf(f);

f := -3/4*exp(1)+33/4*exp(-1)+x*(235/4*exp(-1)-215/28*exp(1))+x^2*(15/4*exp(1)-105/4*exp(-1))+x^3*(1075/84*exp(1)-1115/12*exp(-1))

.996294019+.74039599*x+.536721528*x^2+.60540390*x^3

Comparing results with the Legendre polynomials provided by Maple, and we obtain exactly the same result!

> with(orthopoly);

[G, H, L, P, T, U]

> P(1,x);

x

> P(0,x);

1

> P(2,x);

3/2*x^2-1/2

> P(3,x);

5/2*x^3-3/2*x

> f:=int(exp(x)*P(0,x),x=-1..1)/int(P(0,x)*P(0,x),x=-1..1)*P(0,x)+int(exp(x)*P(1,x),x=-1..1)/int(P(1,x)*P(1,x),x=-1..1)*P(1,x)+int(exp(x)*P(2,x),x=-1..1)/int(P(2,x)*P(2,x),x=-1..1)*P(2,x)+int(exp(x)*P(3,x),x=-1..1)/int(P(3,x)*P(3,x),x=-1..1)*P(3,x);

f := 1/2*exp(1)-1/2*exp(-1)+3*exp(-1)*x+5/2*(exp(1)-7*exp(-1))*(3/2*x^2-1/2)+7/2*(-5*exp(1)+37*exp(-1))*(5/2*x^3-3/2*x)

> evalf(f);

> .9962940180+.9979548790*x+.5367215250*x^2+.1761390750*x^3;

.9962940180+.9979548790*x+.5367215250*x^2+.1761390750*x^3

This is the same result as the quantity found with the method of our orthogonal polynomials.

.9962940173+.9979548182*x+.5367215271*x^2+.1761391763*x^3

================================================

Project Part 3:

Find the Fourier series for y = x^3 on [-Pi..Pi].

We'll start out this way, for the sake of pedagogy, but we can short-cut these steps.

> restart;

> a0:=1/(sqrt(2)*Pi)*Int(x^3,x=-Pi..Pi);

a0 := 1/2*2^(1/2)/Pi*Int(x^3,x = -Pi .. Pi)

> a0:=1/(sqrt(2)*Pi)*int(x^3,x=-Pi..Pi);

a0 := 0

> a1:=1/(Pi)*Int(x^3*cos(x),x=-Pi..Pi);

a1 := 1/Pi*Int(x^3*cos(x),x = -Pi .. Pi)

> a1:=1/(Pi)*int(x^3*cos(x),x=-Pi..Pi);

a1 := 0

>

> a2:= 1/(Pi)*Int(x^3*cos(2*x),x=-Pi..Pi);

a2 := 1/Pi*Int(x^3*cos(2*x),x = -Pi .. Pi)

> a2:= 1/(Pi)*int(x^3*cos(2*x),x=-Pi..Pi);

a2 := 0

> a3:=(1/Pi)*int(x^3*cos(3*x),x=-Pi.. Pi);

a3 := 0

> #Notice: all even functions (cosines) have coeffs of 0.

> b1:=(1/Pi)*Int(x^3*sin(x),x=-Pi..Pi);

b1 := 1/Pi*Int(x^3*sin(x),x = -Pi .. Pi)

> b1:=(1/Pi)*int(x^3*sin( x),x=-Pi..Pi);

b1 := 1/Pi*(2*Pi^3-12*Pi)

> b2:=(1/Pi)*int(x^3*sin(2*x),x=-Pi..Pi);

b2 := 1/Pi*(-Pi^3+3/2*Pi)

> b3:=(1/Pi)*int(x^3*sin(3*x),x=-Pi..Pi);

b3 := 1/Pi*(2/3*Pi^3-4/9*Pi)

> #In general,

> b(n):=int(x^3*sin(n*x),x=-Pi..Pi); a(n):=int(x^3*cos(n*x),x=-Pi..Pi);

>

b(n) := -2*(n^3*Pi^3*cos(Pi*n)-3*n^2*Pi^2*sin(Pi*n)+6*sin(Pi*n)-6*Pi*n*cos(Pi*n))/n^4

a(n) := 0

Notice that all the cos terms are zero. (a(n) = 0). This is because y = x^3 is ODD, and the Fourier series decides it doesn't need any even terms like cos, in its approximation. Here's the shortcut.

> ourseries:=a0/(sqrt(2)*Pi)+sum(a('n')*cos('n'*x)+b('n')*sin('n'*x),'n'=1..10);

ourseries := -2*(-Pi^3+6*Pi)*sin(x)-1/8*(8*Pi^3-12*Pi)*sin(2*x)-2/81*(-27*Pi^3+18*Pi)*sin(3*x)-1/128*(64*Pi^3-24*Pi)*sin(4*x)-2/625*(-125*Pi^3+30*Pi)*sin(5*x)-1/648*(216*Pi^3-36*Pi)*sin(6*x)-2/2401*(-343*Pi^3+42*Pi)*sin(7*x)-1/2048*(512*Pi^3-48*Pi)*sin(8*x)-2/6561*(-729*Pi^3+54*Pi)*sin(9*x)-1/5000*(1000*Pi^3-60*Pi)*sin(10*x)

> with(plots):m1:=plot(ourseries, x = -2*Pi..2*Pi,color=red, thickness=2):m2:=plot(x^3,x=- 2*Pi.. 2*Pi,color=black):display(m1,m2);

Notice how the Fourier series (red) gives a wavey-approximation to x^3 on the interval [-Pi, Pi], then repeats and repeats itself. Let's see what happens if we use more terms in the Fourier series.

>

> ourseries2:=sum(b('n')*sin('n'*x),'n'=1..15);

ourseries2 := -2*(-Pi^3+6*Pi)*sin(x)-1/8*(8*Pi^3-12*Pi)*sin(2*x)-2/81*(-27*Pi^3+18*Pi)*sin(3*x)-1/128*(64*Pi^3-24*Pi)*sin(4*x)-2/625*(-125*Pi^3+30*Pi)*sin(5*x)-1/648*(216*Pi^3-36*Pi)*sin(6*x)-2/2401*(-343*Pi^3+42*Pi)*sin(7*x)-1/2048*(512*Pi^3-48*Pi)*sin(8*x)-2/6561*(-729*Pi^3+54*Pi)*sin(9*x)-1/5000*(1000*Pi^3-60*Pi)*sin(10*x)-2/14641*(-1331*Pi^3+66*Pi)*sin(11*x)-1/10368*(1728*Pi^3-72*Pi)*sin(12*x)-2/28561*(-2197*Pi^3+78*Pi)*sin(13*x)-1/19208*(2744*Pi^3-84*Pi)*sin(14*x)-2/50625*(-3375*Pi^3+90*Pi)*sin(15*x)

> with(plots):ourseries2plot:=plot(-2*(-Pi^3+6*Pi)*sin(x)-1/8*(8*Pi^3-12*Pi)*sin(2*x)-2/81*(-27*Pi^3+18*Pi)*sin(3*x)-1/128*(64*Pi^3-24*Pi)*sin(4*x)-2/625*(-125*Pi^3+30*Pi)*sin(5*x)-1/648*(216*Pi^3-36*Pi)*sin(6*x)-2/2401*(-343*Pi^3+42*Pi)*sin(7*x)-1/2048*(512*Pi^3-48*Pi)*sin(8*x)-2/6561*(-729*Pi^3+54*Pi)*sin(9*x)-1/5000*(1000*Pi^3-60*Pi)*sin(10*x)-2/14641*(-1331*Pi^3+66*Pi)*sin(11*x)-1/10368*(1728*Pi^3-72*Pi)*sin(12*x)-2/28561*(-2197*Pi^3+78*Pi)*sin(13*x)-1/19208*(2744*Pi^3-84*Pi)*sin(14*x)-2/50625*(-3375*Pi^3+90*Pi)*sin(15*x),x=-2*Pi..2*Pi,color=blue):plot(m2,ourseries2,x=-2*Pi..2*Pi):display(m1,ourseries2plot,m2);

Error, (in plot) invalid arguments

The approximation improves, but to the eye, not by too much. And the convergence (of the red graph to the black) is NOT GOOD at the endpoints, -Pi and Pi. This is the problem with Fourier Series. When you are interpreting waves, it is not good if the boundaries of the waves are bad. This leads to bad "edges" when you're viewing an x-ray, a CT scan, a photograph, or trying to hear consonants over the phone. Wavelets are creating better "edges."

Now find the Fourier series for y = exp(x).

> restart:with(plots):

> a0:=1/(sqrt(2)*Pi)*int(exp(x),x=-Pi..Pi);

a0 := 1/2*2^(1/2)/Pi*(exp(Pi)-exp(-Pi))

>

> b(n):=1/Pi*int(exp(x)*sin(n*x),x=-Pi..Pi); a(n):=1/Pi*int(exp(x)*cos(n*x),x=-Pi..Pi);

>

b(n) := 1/Pi*(-n*exp(Pi)*cos(Pi*n)+exp(Pi)*sin(Pi*n)+n*exp(-Pi)*cos(Pi*n)+exp(-Pi)*sin(Pi*n))/(1+n^2)

a(n) := -1/Pi*(-exp(Pi)*cos(Pi*n)-n*exp(Pi)*sin(Pi*n)+exp(-Pi)*cos(Pi*n)-n*exp(-Pi)*sin(Pi*n))/(1+n^2)

> ourseries:=a0/(sqrt(2)*Pi)+sum(a('n')*cos('n'*x)+b('n')*sin('n'*x),'n'=1..10);

ourseries := 1/2/Pi^2*(exp(Pi)-exp(-Pi))-1/2/Pi*(exp(Pi)-exp(-Pi))*cos(x)+1/2/Pi*(exp(Pi)-exp(-Pi))*sin(x)-1/5/Pi*(-exp(Pi)+exp(-Pi))*cos(2*x)+1/5/Pi*(-2*exp(Pi)+2*exp(-Pi))*sin(2*x)-1/10/Pi*(exp(Pi)-exp(-Pi))*cos(3*x)+1/10/Pi*(3*exp(Pi)-3*exp(-Pi))*sin(3*x)-1/17/Pi*(-exp(Pi)+exp(-Pi))*cos(4*x)+1/17/Pi*(-4*exp(Pi)+4*exp(-Pi))*sin(4*x)-1/26/Pi*(exp(Pi)-exp(-Pi))*cos(5*x)+1/26/Pi*(5*exp(Pi)-5*exp(-Pi))*sin(5*x)-1/37/Pi*(-exp(Pi)+exp(-Pi))*cos(6*x)+1/37/Pi*(-6*exp(Pi)+6*exp(-Pi))*sin(6*x)-1/50/Pi*(exp(Pi)-exp(-Pi))*cos(7*x)+1/50/Pi*(7*exp(Pi)-7*exp(-Pi))*sin(7*x)-1/65/Pi*(-exp(Pi)+exp(-Pi))*cos(8*x)+1/65/Pi*(-8*exp(Pi)+8*exp(-Pi))*sin(8*x)-1/82/Pi*(exp(Pi)-exp(-Pi))*cos(9*x)+1/82/Pi*(9*exp(Pi)-9*exp(-Pi))*sin(9*x)-1/101/Pi*(-exp(Pi)+exp(-Pi))*cos(10*x)+1/101/Pi*(-10*exp(Pi)+10*exp(-Pi))*sin(10*x)

> with(plots): ourseries:=a0/sqrt(2)+sum(a('n')*cos('n'*x)+b('n')*sin('n'*x),'n'=1..10);

ourseries := 1/2/Pi*(exp(Pi)-exp(-Pi))-1/2/Pi*(exp(Pi)-exp(-Pi))*cos(x)+1/2/Pi*(exp(Pi)-exp(-Pi))*sin(x)-1/5/Pi*(-exp(Pi)+exp(-Pi))*cos(2*x)+1/5/Pi*(-2*exp(Pi)+2*exp(-Pi))*sin(2*x)-1/10/Pi*(exp(Pi)-exp(-Pi))*cos(3*x)+1/10/Pi*(3*exp(Pi)-3*exp(-Pi))*sin(3*x)-1/17/Pi*(-exp(Pi)+exp(-Pi))*cos(4*x)+1/17/Pi*(-4*exp(Pi)+4*exp(-Pi))*sin(4*x)-1/26/Pi*(exp(Pi)-exp(-Pi))*cos(5*x)+1/26/Pi*(5*exp(Pi)-5*exp(-Pi))*sin(5*x)-1/37/Pi*(-exp(Pi)+exp(-Pi))*cos(6*x)+1/37/Pi*(-6*exp(Pi)+6*exp(-Pi))*sin(6*x)-1/50/Pi*(exp(Pi)-exp(-Pi))*cos(7*x)+1/50/Pi*(7*exp(Pi)-7*exp(-Pi))*sin(7*x)-1/65/Pi*(-exp(Pi)+exp(-Pi))*cos(8*x)+1/65/Pi*(-8*exp(Pi)+8*exp(-Pi))*sin(8*x)-1/82/Pi*(exp(Pi)-exp(-Pi))*cos(9*x)+1/82/Pi*(9*exp(Pi)-9*exp(-Pi))*sin(9*x)-1/101/Pi*(-exp(Pi)+exp(-Pi))*cos(10*x)+1/101/Pi*(-10*exp(Pi)+10*exp(-Pi))*sin(10*x)

> evalf(ourseries);

3.676077911-3.676077911*cos(x)+3.676077911*sin(x)+1.470431164*cos(2.*x)-2.940862328*sin(2.*x)-.7352155821*cos(3.*x)+2.205646746*sin(3.*x)+.4324797542*cos(4.*x)-1.729919017*sin(4.*x)-.2827752239*cos(5.*x)+1.413876119*sin(5.*x)+.1987069141*cos(6.*x)-1.192241484*sin(6.*x)-.1470431164*cos(7.*x)+1.029301815*sin(7.*x)+.1131100895*cos(8.*x)-.9048807163*sin(8.*x)-.8966043683e-1*cos(9.*x)+.8069439315*sin(9.*x)+.7279362199e-1*cos(10.*x)-.7279362199*sin(10.*x)

>

> with(plots):m1:=plot(ourseries,x = -2*Pi..2*Pi,color=red, thickness=2):m2:=plot(exp(x),x=- 2*Pi.. 2*Pi,color=black):display(m1,m2);

Looking up close:

> m1:=plot(ourseries,x = -Pi..Pi,color=red, thickness=2):m2:=plot(exp(x),x=-Pi.. Pi,color=black):display(m1,m2);

========================================================

Project Part 4:

Maple has trouble finding Fourier series for f(x) = sin^2(x) * cos(x) ; division by 0 is the problem.

> restart:with(plots):

> a0:=1/(sqrt(2)*Pi)*int(sin^2(x)*cos(x),x=-Pi..Pi);

a0 := 0

>

> b(n):=1/Pi*int(sin^2(x)*cos(x)*sin(n*x),x=-Pi..Pi); a(n):=1/Pi*int(sin^2(x)*cos(x)*cos(n*x),x=-Pi..Pi);

>

b(n) := 0

a(n) := -2/Pi/(1+n)/(-1+n)*sin^2*n*sin(Pi*n)

It would appear that a1 is undefined, and a3 = 0.

> with(plots): ourseries:=a0/sqrt(2)+sum(a('n')*cos('n'*x)+b('n')*sin('n'*x),'n'=1..10);

Error, (in sum) division by zero

> ourseries:=a0/sqrt(2)+sum(a('n')*cos('n'*x)+b('n')*sin('n'*x),'n'=2..10);

ourseries := 0

However, when we find a1 and a3 by themselves (outside of the grand formula)

> a1:=1/Pi*int((sin(x))^2*cos(x)*cos(x),x=-Pi..Pi);

> a3:=1/Pi*int((sin(x))^2*cos(x)*cos(3*x),x=-Pi..Pi);

a1 := 1/4

a3 := -1/4

So we see these terms are not 0, as the general formula predicted.

>

> restart:with(plots):m1:=plot(1/4*cos(x)-1/4*cos(3*x),x = -2*Pi..2*Pi,color=red, thickness=2):m2:=plot((sin(x))^2*cos(x),x=- 2*Pi..2*Pi,color=black):display(m1);display(m2);

> with(plots):

=========================================

Project Part 5:

The Taylor series for y = exp(-1/x^2) gives 0 (or Maple says, "undefined"). Use Legendre polynomials (up to 3rd degree) to approximate this function on x = -2..2.

> restart:with(linalg):with(plots):taylor(exp(-1/x^2),x=0,4);

Warning, new definition for norm

Warning, new definition for trace

Error, (in series/exp) unable to compute series

> limit(exp(-1/x^2),x=0);

0

> deriv1:=limit(exp(-1/h^2)/h,h=0);

deriv1 := 0

> g:=exp(-1/x^2):a1:=evalf(int(exp(-1/x^2)*f11,x=-1..1));

a1 := .178147710*f11

> a2:=evalf(int(g*f22,x=-1..1));

a2 := .178147710*f22

> a3:=evalf(int(g*f33,x=-1..1));

a3 := .178147710*f33

> a4:=evalf(int(g*f44,x=-1..1));

a4 := .178147710*f44

> g:=simplify(a1*f11+a2*f22+a3*f33+a4*f44);

g := .1781477100*f11^2+.1781477100*f22^2+.1781477100*f33^2+.1781477100*f44^2

>

> m1:=plot(-.3674848897e-1+.3774670343*x^2,x=-1..1):m2:=plot(exp(-1/x^2),x=-1..1,color=blue):display(m1,m2);