Fourierovy_rady_a_ortogonalni_polynomy.mws

Nekonecne rady s programem Maple

Fourierovy rady a ortogonalni polynomy

Karel Srot, xsrot@math.muni.cz

Zapisnik ilustruje vybrane typy ortogonalnich polynomu a vypocet Fourierovych rad vzhledem k systemu techto polynomu.

>   

>    restart;

Procedura automatizuje Grammuv-Schmidtuv ortogonalizacni proces.

>    GramSchmidt:=proc(base, inner_product) local ortonormbase, norm, i, proj, new_vector;

>      ortonormbase:=NULL;

>      norm:=a->sqrt(inner_product(a,a));  

>      for i from 1 to nops(base) do

>        if i = 1 then

>          ortonormbase:=[simplify(base[i]/norm(base[i]))];  

>        else

>          proj:=sum(inner_product(base[i], ortonormbase[j])/norm(ortonormbase[j])*ortonormbase[j], j=1..(i-1));

>          new_vector:=base[i]-proj;

>          ortonormbase:=[op(ortonormbase), simplify(new_vector/norm(new_vector))];

>        end if;

>      end do;

>      RETURN(ortonormbase);

>    end:

>   

>   

Priklad:

Pomici procedury Gramschmidt  spocitejte prvnich pet Legendrovych polynomu a nakreslete jejich grafy.

>    f:=(a,b)->int(a*b,x=-1..1);

f := (a, b) -> int(a*b,x = -1 .. 1)

>    GramSchmidt([1, x, x^2, x^3, x^4], f);

[1/2*2^(1/2), 1/2*x*6^(1/2), 1/4*(3*x^2-1)*10^(1/2), 1/4*x*(5*x^2-3)*14^(1/2), 3/16*(35*x^4+3-30*x^2)*2^(1/2)]

>    L:=map(pol->expand(pol/eval(pol,x=1)), %);

L := [1, x, 3/2*x^2-1/2, 5/2*x^3-3/2*x, 35/8*x^4+3/8-15/4*x^2]

>    pol:=L[3];

pol := 3/2*x^2-1/2

Overime ortogonalitu ziskanych polynomu

>    for pol2 in L do

>     'f'(pol, pol2) = f(pol,pol2);

>    end do;

f(3/2*x^2-1/2,1) = 0

f(3/2*x^2-1/2,x) = 0

f(3/2*x^2-1/2,3/2*x^2-1/2) = 2/5

f(3/2*x^2-1/2,5/2*x^3-3/2*x) = 0

f(3/2*x^2-1/2,35/8*x^4+3/8-15/4*x^2) = 0

Vykreslime graf

>    plot(L,x=-1..1);

[Maple Plot]

>   

>   

Priklad:

Napiste procedury pocitajici Legendrovy polynomy pomoci vztahu

   P[n](x) = 1/(2^n*n!)*Diff((x^2-1)^n,`$`(x,n))   

  a rekurentni formule

   P[n+1](x) = (2*n+1)/(n+1)*x*P[n](x)-n/(n+1)*p[n-1](x)

a overte jejich spravnost pomoci procedur LegendreP  a orthopoly/P .

>    Leg1:=proc(n::nonnegint,x);

>      if n=0 then

>        1

>      else

>        expand(1/(2^n*n!)*diff((x^2-1)^n,x$n));

>      end if;

>    end proc:

>   

Upravime rekurentni vzorec ze zadani na vhodnejsi tvar

>    P(n+1,x)=(2*n+1)/(n+1)*x*P(n,x)-n/(n+1)*P(n-1,x);

P(n+1,x) = (2*n+1)/(n+1)*x*P(n,x)-n/(n+1)*P(n-1,x)

>    algsubs(n+1=k,%);

P(k,x) = x*P(-1+k,x)*(-1+2*k)/k-P(-2+k,x)*(-1+k)/k

>    Leg2:=proc(k::nonnegint,x) option remember;

>      if k=0 then

>        1

>      elif k=1 then

>        x

>      else

>        expand((2*k-1)/k*x*Leg2(k-1,x)-(k-1)/k*Leg2(k-2,x));

>      end if;

>    end proc:

>   

>    Leg1(8,x);

6435/128*x^8-3003/32*x^6+3465/64*x^4-315/32*x^2+35/128

>    Leg2(9,x);

12155/128*x^9-6435/32*x^7+9009/64*x^5-1155/32*x^3+315/128*x

overime spravnost predchozich vysledku pomoci procedur LegendreP  a orthopoly/P .

>    expand(LegendreP(8,x));

6435/128*x^8-3003/32*x^6+3465/64*x^4-315/32*x^2+35/128

>    with(orthopoly);

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

>    P(9,x);

12155/128*x^9-6435/32*x^7+9009/64*x^5-1155/32*x^3+315/128*x

>   

>   

Priklad:

Spocitejte prvnich pet  koeficientu Fourierovy rady funkce sin(Pi*x)  vzhledem k systemu tvoreneho Legendrovymi polynomy na intervalu [-1,1]. Dale nakreslete graf takto ziskaneho castecneho souctu spolu s grafem funkce sin(Pi*x) .

>    f:=x->sin(Pi*x);

f := x -> sin(Pi*x)

>    g:=(a,b)->int(a*b,x=-1..1);

g := (a, b) -> int(a*b,x = -1 .. 1)

>    cn:=n->g(LegendreP(n,x),f(x))/g(LegendreP(n,x),LegendreP(n,x));

cn := n -> g(LegendreP(n,x),f(x))/g(LegendreP(n,x),LegendreP(n,x))

>    N:=6:

>    for i from 0 to (N-1) do

>      print(c[i]=cn(i));

>    end do;

c[0] = 0

c[1] = 3/Pi

c[2] = 0

c[3] = 7*(-15+Pi^2)/Pi^3

c[4] = 0

c[5] = 11*(Pi^4+945-105*Pi^2)/Pi^5

>    S:=sum(cn(k)*LegendreP(k,x),k=0..N-1);

S := 3/Pi*x+7*(-15+Pi^2)/Pi^3*LegendreP(3,x)+11*(Pi^4+945-105*Pi^2)/Pi^5*LegendreP(5,x)

>    pol:=expand(S);

pol := 105/8/Pi*x+39375/4/Pi^3*x^3-16065/8/Pi^3*x-315/4/Pi*x^3+693/8/Pi*x^5+654885/8/Pi^5*x^5+155925/8/Pi^5*x-363825/4/Pi^5*x^3-72765/8/Pi^3*x^5

>    evalf(pol);

3.10346042*x-4.81438833*x^3+1.7269051*x^5

>    g1:=plot(f(x),x=-1..1,thickness=4,color=green,scaling=constrained):

>    g2:=plot(pol,x=-1..1,thickness=2,color=red,scaling=constrained):

>    plots[display](g2,g1);

[Maple Plot]

Pro srovnani nakreslime take graf Maclaurinova polynomu 3. stupne.

>    tayl:=convert(taylor(f(x),x=0,N),polynom);

tayl := Pi*x-1/6*Pi^3*x^3+1/120*Pi^5*x^5

>    g3:=plot(tayl,x=-1..1,thickness=2,color=blue,scaling=constrained):

>    plots[display](g3,g1);

[Maple Plot]

>   

>   

Priklad:

Nakreslete grafy prvnich peti Cebysevovych polynomu prvniho a druheho druhu.

Cebysevovy polynomy prvniho druhu:

definujeme skalarni soucin

>    g:=(a,b)->int((1-x^2)^(-1/2)*a*b,x=-1..1);

g := (a, b) -> int(1/(1-x^2)^(1/2)*a*b,x = -1 .. 1)

>    N:=5:

>    GramSchmidt([seq(x^n, n=0..N-1)], g);

[1/(Pi^(1/2)), x*2^(1/2)/Pi^(1/2), (2*x^2-1)*2^(1/2)/Pi^(1/2), x*(4*x^2-3)*2^(1/2)/Pi^(1/2), (8*x^4+1-8*x^2)*2^(1/2)/Pi^(1/2)]

Podobne jako v pripade Legendrovych polynomu ziskane normovane polynomy normalizujeme pomoci vztahu   T[n](1) = 1  

>    L:=map(pol->expand(pol/eval(pol,x=1)), %);

L := [1, x, 2*x^2-1, 4*x^3-3*x, 8*x^4+1-8*x^2]

>    plot(L,x=-1..1,thickness=2);

[Maple Plot]

Cebysevovy polynomy druheho druhu:

definujeme skalarni soucin

>    h:=(a,b)->int((1-x^2)^(1/2)*a*b,x=-1..1);

h := (a, b) -> int((1-x^2)^(1/2)*a*b,x = -1 .. 1)

>    N:=5:

>    GramSchmidt([seq(x^n, n=0..N-1)], h);

[2^(1/2)/Pi^(1/2), 2*x*2^(1/2)/Pi^(1/2), (4*x^2-1)*2^(1/2)/Pi^(1/2), 4*x*(2*x^2-1)*2^(1/2)/Pi^(1/2), (16*x^4+1-12*x^2)*2^(1/2)/Pi^(1/2)]

Cebysevovy polynomy druheho druhu normalizujeme pomoci vztahu U[n](1) = n+1  

>    L2:=expand([seq(%[i]/eval(%[i],x=1)*i,i=1..N)]);

L2 := [1, 2*x, 4*x^2-1, 8*x^3-4*x, 16*x^4+1-12*x^2]

>    plot(L2,x=-1..1,thickness=2);

[Maple Plot]

>   

>   

Priklad:

Aproximujte funkci abs(x) na intervalu [-1,1] polynomem pateho stupne, jez je castecnym souctem jeji Fourierovy rady vzhledem k systemu tvoreneho

Cebysevovymi polynomy prvniho druhu.

>    f:=x->abs(x);

f := x -> abs(x)

>    g:=(a,b)->int(a*b/sqrt(1-x^2),x=-1..1);

g := (a, b) -> int(a*b/sqrt(1-x^2),x = -1 .. 1)

>    cn:=n->g(ChebyshevT(n,x),f(x))/g(ChebyshevT(n,x),ChebyshevT(n,x));

cn := n -> g(ChebyshevT(n,x),f(x))/g(ChebyshevT(n,x),ChebyshevT(n,x))

>    N:=6:

>    for i from 0 to (N-1) do

>      print(value(c[i]=cn(i)));

>    end do;

c[0] = 2/Pi

c[1] = 0

c[2] = 4/3/Pi

c[3] = 0

c[4] = -4/15/Pi

c[5] = 0

>    S:=sum(cn(k)*ChebyshevT(k,x),k=0..N-1);

S := 2/Pi*ChebyshevT(0,x)+4/3*1/Pi*ChebyshevT(2,x)-4/15*1/Pi*ChebyshevT(4,x)

>    pol:=expand(S);

pol := 2/5/Pi+24/5/Pi*x^2-32/15/Pi*x^4

>    evalf(pol);

.1273239544+1.527887453*x^2-.6790610902*x^4

>   

>    g1:=plot(f(x),x=-1..1,color=green,thickness=3):

>    g2:=plot(pol,x=-1..1,color=red,thickness=2):

>    plots[display](g2,g1);

[Maple Plot]

>   

>    minmax_f:=numapprox[minimax](f, -1..1,[N-1,0]);

minmax_f := x -> .679162859e-1+(-.2472105295e-7+(1.928563325+(.4216266422e-7+(-1.063933613-.1744161126e-7*x)*x)*x)*x)*x

>    minmaxpol:=expand(minmax_f(x));

minmaxpol := .679162859e-1-.2472105295e-7*x+1.928563325*x^2+.4216266422e-7*x^3-1.063933613*x^4-.1744161126e-7*x^5

>    g3:=plot(minmaxpol,x=-1..1,color=blue,thickness=2):

>    plots[display](g3,g2,g1);

[Maple Plot]

>   

>    g4:=plot(minmaxpol-f(x),x=-1..1,color=blue,thickness=2):

>    g5:=plot(pol-f(x),x=-1..1,color=red,thickness=2):

>    plots[display](g4,g5,scaling=constrained);

[Maple Plot]

>   

>   

Priklad:

Aproximujte funkci exp(x) na intervalu [-1,1] polynomem tretiho stupne, jez je castecnym souctem jeji Fourierovy rady vzhledem k systemu tvoreneho

Cebysevovymi polynomy prvniho druhu.

>    f:=x->exp(x);

f := x -> exp(x)

>    g:=(a,b)->int(a*b/sqrt(1-x^2),x=-1..1);

g := (a, b) -> int(a*b/sqrt(1-x^2),x = -1 .. 1)

>    cn:=n->g(ChebyshevT(n,x),f(x))/g(ChebyshevT(n,x),ChebyshevT(n,x));

cn := n -> g(ChebyshevT(n,x),f(x))/g(ChebyshevT(n,x),ChebyshevT(n,x))

>    N:=4:

>    for i from 0 to (N-1) do

>      print(c[i]=cn(i));

>    end do;

c[0] = 1/Pi*int(ChebyshevT(0,x)*exp(x)/(1-x^2)^(1/2),x = -1 .. 1)

c[1] = 2/Pi*int(ChebyshevT(1,x)*exp(x)/(1-x^2)^(1/2),x = -1 .. 1)

c[2] = 2/Pi*int(ChebyshevT(2,x)*exp(x)/(1-x^2)^(1/2),x = -1 .. 1)

c[3] = 2/Pi*int(ChebyshevT(3,x)*exp(x)/(1-x^2)^(1/2),x = -1 .. 1)

>    S:=sum(cn(k)*ChebyshevT(k,x),k=0..N-1):

>    pol:=evalf(expand(S));

pol := .9945705384+.9973076585*x+.5429906792*x^2+.1773473994*x^3

>   

>    g1:=plot(f(x),x=-1..1,color=green,thickness=3):

>    g2:=plot(pol,x=-1..1,color=red,thickness=2):

>    plots[display](g2,g1);

[Maple Plot]

>   

>    minmax_f:=numapprox[minimax](f, -1..1,[N-1,0]);

minmax_f := x -> .9945718793+(.9956674823+(.5429879654+.1795337116*x)*x)*x

>    minmaxpol:=expand(minmax_f(x));

minmaxpol := .9945718793+.9956674823*x+.5429879654*x^2+.1795337116*x^3

>    g3:=plot(minmaxpol,x=-1..1,color=blue,thickness=2):

>    plots[display](g3,g2,g1);

[Maple Plot]

>   

>    g4:=plot(abs(minmaxpol-f(x)),x=-1..1,color=blue,thickness=2):

>    g5:=plot(abs(pol-f(x)),x=-1..1,color=red,thickness=2):

>    plots[display](g4,g5);

[Maple Plot]

>    evalf(maximize(abs(minmaxpol-f(x)),x=-1..1));

.55359407e-2

>    evalf(maximize(abs(pol-f(x)),x=-1..1));

.6065553e-2

>   

>