>

Approximation by Trigonometric Polynomials.

Ingnore this stuff:

> restart; with(plots): with(LinearAlgebra):

Warning, the name changecoords has been redefined

> ip := proc(f,g) int(f*g, x=-Pi..Pi);

> end;

ip := proc (f, g) int(g*f, x = -Pi .. Pi) end proc

> nrm2 := proc(f) ip(f,f); end;

nrm2 := proc (f) ip(f, f) end proc

> nrm := proc(f) sqrt(nrm2(f)); end;

nrm := proc (f) sqrt(nrm2(f)) end proc

> append := proc(x, llist) [op(llist), x]; end;

append := proc (x, llist) [op(llist), x] end proc

We're working in the space of continuse functions on the interval [-Pi, Pi] .  The inner product is

> ip(f(x),g(x));

int(g(x)*f(x), x = -Pi .. Pi)

n and m will denote postive integers

> assume(n,posint);

> assume(m, posint);

>

Consider the functions 1, cos(nx), sin(nx).  Here's an illustration of what changing n does.

> animate(plot,[sin(n*x), x=-Pi..Pi], n=1..10);

[Plot]

>

Compute the inner product

> ip(sin(n*x), sin(m*x));

0

L:et's see why this is true:

> Int(sin(n*x)*sin(m*x), x) = int(sin(n*x)*sin(m*x), x);

Int(sin(m*x)*sin(n*x), x) = 1/2*sin((-m+n)*x)/(-m+n)-1/2*sin((m+n)*x)/(m+n)

> rh := rhs(%);

rh := 1/2*sin((-m+n)*x)/(-m+n)-1/2*sin((m+n)*x)/(m+n)

> Int(sin(n*x)*sin(m*x), x=-Pi..Pi) = subs(x=Pi, rh)-subs(x=-Pi, rh);

Int(sin(m*x)*sin(n*x), x = -Pi .. Pi) = 1/2*sin(Pi*(-m+n))/(-m+n)-1/2*sin(Pi*(m+n))/(m+n)-1/2*sin(-Pi*(-m+n))/(-m+n)+1/2*sin(-Pi*(m+n))/(m+n)

> value(rhs(%));

0

> ip(cos(n*x), cos(m*x));

0

> ip(cos(n*x), sin(m*x));

0

> ip(1, cos(n*x));

0

> ip(1, sin(n*x));

0

Thus, these functions are othogonal.  To make them ON, we normalize them.

> nrm(1);

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

> u1 := 1/nrm(1);

u1 := 1/2*2^(1/2)/Pi^(1/2)

> nrm(cos(n*x));

Pi^(1/2)

> nrm(sin(n*x));

Pi^(1/2)

> basis := [u1];

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

Here is an ON set of functions:

> for j from 1 to 5 do basis := append(cos(j*x)/sqrt(Pi), basis); basis := append(sin(j*x)/sqrt(Pi), basis); od:

> basis;

[1/2*2^(1/2)/Pi^(1/2), cos(x)/Pi^(1/2), sin(x)/Pi^(1/2), cos(2*x)/Pi^(1/2), sin(2*x)/Pi^(1/2), cos(3*x)/Pi^(1/2), sin(3*x)/Pi^(1/2), cos(4*x)/Pi^(1/2), sin(4*x)/Pi^(1/2), cos(5*x)/Pi^(1/2), sin(5*x)/P...

> for i from 1 to nops(basis) do

> for j from i to nops(basis) do print(i,j, ip(basis[i], basis[j])); od;

> od;

1, 1, 1

1, 2, 0

1, 3, 0

1, 4, 0

1, 5, 0

1, 6, 0

1, 7, 0

1, 8, 0

1, 9, 0

1, 10, 0

1, 11, 0

2, 2, 1

2, 3, 0

2, 4, 0

2, 5, 0

2, 6, 0

2, 7, 0

2, 8, 0

2, 9, 0

2, 10, 0

2, 11, 0

3, 3, 1

3, 4, 0

3, 5, 0

3, 6, 0

3, 7, 0

3, 8, 0

3, 9, 0

3, 10, 0

3, 11, 0

4, 4, 1

4, 5, 0

4, 6, 0

4, 7, 0

4, 8, 0

4, 9, 0

4, 10, 0

4, 11, 0

5, 5, 1

5, 6, 0

5, 7, 0

5, 8, 0

5, 9, 0

5, 10, 0

5, 11, 0

6, 6, 1

6, 7, 0

6, 8, 0

6, 9, 0

6, 10, 0

6, 11, 0

7, 7, 1

7, 8, 0

7, 9, 0

7, 10, 0

7, 11, 0

8, 8, 1

8, 9, 0

8, 10, 0

8, 11, 0

9, 9, 1

9, 10, 0

9, 11, 0

10, 10, 1

10, 11, 0

11, 11, 1

> proj := proc(f, k)

> add(ip(f, basis[i])*basis[i], i=1..k);

> end;

proj := proc (f, k) add(ip(f, basis[i])*basis[i], i = 1 .. k) end proc

proj(f, k) computes the projection of f onto the space spanned by the first k basis elements.   

Consider the following function:

> myfun := expand(x*(x-2)*(x+1));

myfun := x^3-x^2-2*x

> plot1 := plot(myfun, x=-Pi..Pi): plot1;

[Plot]

> pr[1] := proj(myfun, 1);

pr[1] := -1/3*Pi^2

> pr[2] := proj(myfun, 2);

pr[2] := -1/3*Pi^2+4*cos(x)

> pr[3] := proj(myfun, 3);

pr[3] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)

> for j from 4 to nops(basis) do pr[j]:= proj(myfun, j); od;

pr[4] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)-cos(2*x)

pr[5] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)-cos(2*x)-1/2*(-7+2*Pi^2)*sin(2*x)

pr[6] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)-cos(2*x)-1/2*(-7+2*Pi^2)*sin(2*x)+4/9*cos(3*x)

pr[7] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)-cos(2*x)-1/2*(-7+2*Pi^2)*sin(2*x)+4/9*cos(3*x)+2/9*(-8+3*Pi^2)*sin(3*x)

pr[8] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)-cos(2*x)-1/2*(-7+2*Pi^2)*sin(2*x)+4/9*cos(3*x)+2/9*(-8+3*Pi^2)*sin(3*x)-1/4*cos(4*x)

pr[9] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)-cos(2*x)-1/2*(-7+2*Pi^2)*sin(2*x)+4/9*cos(3*x)+2/9*(-8+3*Pi^2)*sin(3*x)-1/4*cos(4*x)-1/16*(-19+8*Pi^2)*sin(4*x)

pr[10] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)-cos(2*x)-1/2*(-7+2*Pi^2)*sin(2*x)+4/9*cos(3*x)+2/9*(-8+3*Pi^2)*sin(3*x)-1/4*cos(4*x)-1/16*(-19+8*Pi^2)*sin(4*x)+4/25*cos(5*x)

pr[11] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)-cos(2*x)-1/2*(-7+2*Pi^2)*sin(2*x)+4/9*cos(3*x)+2/9*(-8+3*Pi^2)*sin(3*x)-1/4*cos(4*x)-1/16*(-19+8*Pi^2)*sin(4*x)+4/25*cos(5*x)+2/125*(-56+25*Pi^2)*sin(5*...pr[11] := -1/3*Pi^2+4*cos(x)+2*(-8+Pi^2)*sin(x)-cos(2*x)-1/2*(-7+2*Pi^2)*sin(2*x)+4/9*cos(3*x)+2/9*(-8+3*Pi^2)*sin(3*x)-1/4*cos(4*x)-1/16*(-19+8*Pi^2)*sin(4*x)+4/25*cos(5*x)+2/125*(-56+25*Pi^2)*sin(5*...

Plot these projections to see how well they approximate the function.

> for j from 1 to nops(basis) do plt[j]:= plot(pr[j], x=-Pi..Pi, color=blue); od:

> for j from 1 to nops(basis) do plt[j] := display([plot1,plt[j]]); od:

> display([plot1, seq(plt[j], j=1..nops(basis))], insequence=true);

[Plot]

Try it with more basis elements!

> basis := [u1];

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

> for j from 1 to 50 do basis := append(cos(j*x)/sqrt(Pi), basis); basis := append(sin(j*x)/sqrt(Pi), basis); od:

> pr[1] := ip(myfun, basis[1])*basis[1];

pr[1] := -1/3*Pi^2

> for j from 2 to nops(basis) do pr[j] := pr[j-1]+ip(myfun, basis[j])*basis[j]; od:

> for j from 1 to nops(basis) do plt[j]:= plot(pr[j], x=-Pi..Pi, color=blue); od:

> for j from 1 to nops(basis) do plt[j] := display([plot1,plt[j]]); od:

> display([plot1, seq(plt[j], j=1..nops(basis))], insequence=true);

[Plot]

Evaluate the norm of the difference between f and the projections:

> for j from 1 to nops(basis) do

> nn2 := evalf(Int((myfun-pr[j])^2,x=-Pi..Pi)); nn := sqrt(nn2); norms[j] := nn; od:

>

> for j from 1 to nops(basis) do print(j, norms[j]); od;

1, 22.59195278

2, 21.45066078

3, 20.40112980

4, 20.32398840

5, 16.89982978

6, 16.88145980

7, 14.57883094

8, 14.57209532

9, 12.97038388

10, 12.96728318

11, 11.78517975

12, 11.78353413

13, 10.86921190

14, 10.86824881

15, 10.13541890

16, 10.13481349

17, 9.531127066

18, 9.530725150

19, 9.022571454

20, 9.022292896

21, 8.587057098

22, 8.586857190

23, 8.208723891

24, 8.208576238

25, 7.876124106

26, 7.876012379

27, 7.580767399

28, 7.580681098

29, 7.316209758

30, 7.316141901

31, 7.077462288

32, 7.077408103

33, 6.860595557

34, 6.860551696

35, 6.662467578

36, 6.662431643

37, 6.480532359

38, 6.480502599

39, 6.312702304

40, 6.312677421

41, 6.157247485

42, 6.157226497

43, 6.012720651

44, 6.012702807

45, 5.877900578

46, 5.877885298

47, 5.751748680

48, 5.751735510

49, 5.633375377

50, 5.633363957

51, 5.522013719

52, 5.522003760

53, 5.416998496

54, 5.416989766

55, 5.317749540

56, 5.317741851

57, 5.223758249

58, 5.223751446

59, 5.134576632

60, 5.134570590

61, 5.049808339

62, 5.049802950

63, 4.969101252

64, 4.969096429

65, 4.892141343

66, 4.892137012

67, 4.818647537

68, 4.818643634

69, 4.748367396

70, 4.748363870

71, 4.681073482

72, 4.681070286

73, 4.616560255

74, 4.616557350

75, 4.554641444

76, 4.554638797

77, 4.495147778

78, 4.495145361

79, 4.437925053

80, 4.437922842

81, 4.382832444

82, 4.382830414

83, 4.329741045

84, 4.329739180

85, 4.278532619

86, 4.278530900

87, 4.229098470

88, 4.229096885

89, 4.181338491

90, 4.181337025

91, 4.135160304

92, 4.135158947

93, 4.090478510

94, 4.090477251

95, 4.047214024

96, 4.047212855

97, 4.005293495

98, 4.005292406

99, 3.964648773

100, 3.964647759

101, 3.925216450

>