Fourierovy_trigonometricke

Nekonecne rady s program em Maple

Fourierov y trigonometricke rady

Karel Srot, xsrot@math.muni.cz

Zapisnik ilustruje vypocty Fourierovych rad a jejich aplikace pri reseni diferencialnich rovnich ci vysetrovani kmitajici struny.

>

vypocet metodou krok za krokem

Priklad:

Spocitejte Fourierovu radu funkce definovane na intervalu [ , ] a pomoci animace znazornete konvergenci jejich castecnych souctu.

>

> restart;

> assume(n,integer);

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

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

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

> a0/2+Sum(aN*cos(n*x)+bN*sin(n*x), n=1..infinity);

> four:=(m,x)->a0/2+sum(aN*cos(n*x)+bN*sin(n*x), n=1..m):

> T[3](x)=four(3,x);

> with(plots):

> graf1:=plot(x^2,x=-Pi..Pi,color=aquamarine,thickness=2):

> graf2:=plot(four(3,x),x=-Pi..Pi,color=red):

> display(graf1,graf2);

> clenu:=10:

> anim:=animate(four(m,x),x=-3*Pi..3*Pi,m=0..clenu,color=red,frames=clenu+1,
numpoints=300):

> display(anim, graf1);

> display([graf1,anim],view=[-1.2..1.2,-0.8..1.6]);

> animate(x^2-four(m,x), x=-Pi..Pi, m=1..clenu, frames=clenu, numpoints=300,
view=[-Pi..Pi,-1..1], scaling=constrained, thickness=2);

>

Priklad:

Urcete kvadratickou odchylku prvnich deseti castecnych souctu Fourierovy rady funkce z predchazejiciho prikladu.

> qdev:= (f, g, var) -> evalf(sqrt(int((f-g)^2, var))):

> for i from 0 to 9 do

> qdev(x^2, four(i,x), x=-Pi..Pi);

> end do;

> qdev(x^2, four(50,x), x=-Pi..Pi);

> qdev(x^2, four(100,x), x=-Pi..Pi);

> qdev(x^2, four(200,x), x=-Pi..Pi);

>

Priklad:

Spocitejte Fourierovu radu funkce f(x)=|sin(x)| na intervalu [ , ].

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

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

> aN:=1/Pi*int(f(x)*cos(n*x), x=-Pi..Pi);

> bN:=1/Pi*int(f(x)*sin(n*x), x=-Pi..Pi);

> a0/2+Sum(aN*cos(n*x)+bN*sin(n*x),n=1..infinity);

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

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

> a0/2+Sum(aN*cos(n*x)+bN*sin(n*x),n=2..infinity);

>

Priklad:

Spocitejte Fourierovu radu funkce definovane predpisem

0 -Pi < x < 0

f(x) = {

sin(x) 0 <= x < Pi

na intervalu [ , ].

> f:=x->piecewise(x<0,0,sin(x)):

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

> aN:=1/Pi*int(f(x)*cos(n*x), x=-Pi..Pi);

> bN:=1/Pi*int(f(x)*sin(n*x), x=-Pi..Pi);

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

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

> a0/2+b1*sin(x)+Sum(aN*cos(n*x)+bN*sin(n*x),n=2..infinity);

>

Priklad:

Spocitejte hodnoty koeficientu Fourierovy rady funkce zadane predpisem

pro x z intervalu [-1/2, 1].

> restart;

> assume(n, integer):

> f:=x->1-x^2+x:

> a:=-1/2:

> b:=1:

> L:=abs(a-b):

>

> a0:=2/L*int(f(x), x=a..b);

> aN:=2/L*int(f(x)*cos(2*Pi/L*n*x), x=a..b);

> bN:=2/L*int(f(x)*sin(2*Pi/L*n*x), x=a..b);

> a0/2+Sum(aN*cos(2*Pi/L*n*x)+bN*sin(2*Pi/L*n*x),n=1..infinity);

> four:=(m,x)->a0/2+sum(aN*cos(2*Pi/L*n*x)+bN*sin(2*Pi/L*n*x), n=1..m):

> with(plots):

> graf1:=plot(f(x), x=a..b, color=aquamarine, thickness=2):

> for i in [2,3,5,10] do

> display(graf1, plot(four(i,x), x=-2..2, color=red, numpoints=300),
scaling=constrained);

> end do;

>

vypocet pomoci procedur knihovny FourierTrigSeries

> restart;

> with(FourierTrigSeries);

>

Priklad:

Spocitejte Fourierovu radu funkce f(x)=sgn(cos x) na intervalu a vykreslete grafy castecnych souctu , a spolu s periodickym rozsirenim funkce f .

> f:=signum(cos(x)):

> FRada:=GetFourierSeries(f,x=-Pi..Pi);

> gb:=DrawPeriodicExtension(f,x=-Pi..Pi,-7..7,color=aquamarine,thickness=2,drawlimits):

> g1:=DrawPartialSum(FRada,1,-7..7,color=blue):

> g2:=DrawPartialSum(FRada,3,-7..7,color=green):

> g3:=DrawPartialSum(FRada,5,-7..7):

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

>

Priklad:

Spocitejte Fourierovu radu funkce f(x)=abs(sin x) na intervalu a vykreslete grafy castecnych souctu , , a spolu s periodickym rozsirenim funkce f .

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

> FRada:=GetFourierSeries(f,x=-Pi..Pi);

> gb:=DrawPeriodicExtension(f,x=-Pi..Pi,-7..7,color=aquamarine,thickness=2):

> g1:=DrawPartialSum(FRada,1,-7..7,color=blue):

> g2:=DrawPartialSum(FRada,3,-7..7,color=green):

> g3:=DrawPartialSum(FRada,5,-7..7,color=black):

> g4:=DrawPartialSum(FRada,7,-7..7):

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

>

Priklad:

Spocitejte Fourierovu radu funkce f(x)=x sin(x) a sestrojte animaci znazornujici konvergenci Fourierovy rady.

> f:=x*sin(x);

> FRada:=GetFourierSeries(f,x=-Pi..Pi);

> gb:=DrawPeriodicExtension(f,x=-Pi..Pi,-6..6,thickness=2, color=aquamarine):

> anim:=seq(plots[display](gb, DrawPartialSum(FRada,i,-6..6, numpoints=200)), i=0..10):

> plots[display](anim, insequence=true,scaling=constrained);

>

Priklad:

> f:=x->(Pi-x)/2;

> FRada:=GetFourierSeries(f(x),x=0..2*Pi);

> gb:=DrawPeriodicExtension(f(x),x=0..2*Pi,-Pi..3*Pi,drawlimits,color=aquamarine,thickness=2):

> g1:=DrawPartialSum(FRada, 1, -Pi..3*Pi, color=green):

> g2:=DrawPartialSum(FRada, 2, -Pi..3*Pi, color=blue):

> g3:=DrawPartialSum(FRada, 4, -Pi..3*Pi, color=violet):

> g4:=DrawPartialSum(FRada, 8, -Pi..3*Pi, color=red):

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

>

Priklad:

> f:=cos(a*x);

> GetFourierSeries(f,x=-1..1);

>

Priklad:

> f:=piecewise(x<0,a*x,b*x);

> FRada:=GetFourierSeries(f,x=-Pi..Pi);

>

Priklad:

> f:=piecewise(x<0,1/2*x,2*x);

> FRada:=GetFourierSeries(f,x=-Pi..Pi);

> gb:=DrawPeriodicExtension(f,x=-Pi..Pi, -6..6,drawlimits,color=aquamarine,thickness=2):

> anim:=seq(plots[display](gb, DrawPartialSum(FRada,i,-6..6, numpoints=200)), i=0..15):

> plots[display](anim, insequence=true);

>

Priklad:

> f:=x->x*cos(x);

> FRada:=GetFourierSeries(f(x),x=-Pi..Pi);

> gb:=DrawPeriodicExtension(f(x),x=-Pi..Pi,-6..6,drawlimits, color=aquamarine,thickness=2):

> anim:=seq(plots[display](gb, DrawPartialSum(FRada,i,-6..6, numpoints=250)), i=0..20):

> plots[display](anim, insequence=true,scaling=constrained);

>

Priklad:

> f:=(Pi^2-x^2)^2;

> FRada:=GetFourierSeries(f,x=-Pi..Pi);

> gb:=DrawPeriodicExtension(f,x=-Pi..Pi,-6..6,color=aquamarine,thickness=2):

> anim:=seq(plots[display](gb, DrawPartialSum(FRada,i,-6..6, numpoints=200)), i=0..10):

> plots[display](anim, insequence=true);

>

Priklad:

> f:=x->piecewise(x<0,a,b);

> GetFourierSeries(f(x),x=-Pi..Pi);

> SimplifyCoefficients(%, simplify);

>

Priklad:

> f:=x->arcsin(sin(x));

> FRada:=GetFourierSeries(f(x),x=-Pi..Pi);

> g:=x->piecewise(x<-Pi/2,-x-Pi,x<Pi/2,x,Pi-x):

> GetFourierSeries(g(x),x=-Pi..Pi);

> SimplifyCoefficients(%, simplify);

>

Priklad:

Spocitejte Fourierovu ?adu funkce y=signum(x) pro x z intervalu [-Pi, Pi] a urcete funkci hodnotu castecnych souctu

pro n=10, 50, 100, 500, 1000, 5 000 a 10 000.

Dale nakreslete graf castecneho souctu stupne 20 a graf prislusneho souctu ziskaneho metodou -aproximace.

> f:=x->signum(x):

> FS:=GetFourierSeries(f(x), x=-Pi..Pi);

> for i in [ 10,50,100,500,1000, 5000, 10000 ] do

> evalf(Evaluate(FS, x=(Pi-Pi/i), trunc=i));

> end do;

Vsimneme si hodnoty prekmitu, ktera cini temer 18 %.

> plot(GetPartialSum(FS,20), x=-2*Pi..2*Pi, numpoints=200, scaling=constrained);

> sinc:=x->piecewise(x=0,1,sin(x)/x):

> N:=20:

> s:=Coefficients(FS, 0):

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

> s:=s+sinc(i*Pi/N)*Coefficients(FS,i)[2]*sin(i*x):

> end do:

> plot(s,x=-2*Pi..2*Pi,scaling=constrained,numpoints=200);

>

Priklad:

Spocitejte Fourierovu radu funkce f(x)=x, pro x z intervalu [ , ]. S vyuzitim Besselovy identity spocitejte kvadratickou odchylku castecnych souctu stupnu 100, 500, 1000.

> f:=x->x;

> FS:=GetFourierSeries(f(x), x=-Pi..Pi);

Vypocet pomoci funkce qdev (pocita podle definice kvadraticke odchylky) je casove a pametove narocny.

> qdev:= (f, g, var) -> evalf(sqrt(int((f-g)^2, var))):

>

> s1:=GetPartialSum(FS, 100):

> qdev(f(x), s1, x=-Pi..Pi);

Procedura qdevBE pocita kvadratickou odchylku pomoci Besselovy identity.

> qdevBE:=proc(f, var, N) local FS, L, normf, tmpsum, res, i;

> FS:=GetFourierSeries(f,var);

> normf:=int(f^2,var);

> L:=op(2, op(2,var)) - op(1, op(2,var));

> tmpsum:= L*Coefficients(FS,0)^2;

> for i from 1 to N do

> tmpsum:= tmpsum + L/2*( Coefficients(FS,i)[1]^2 + Coefficients(FS,i)[2]^2) ;

> end do:

> res:=evalf(sqrt(normf - tmpsum));

> return res;

> end proc:

>

> qdevBE(f(x),x=-Pi..Pi,100);

> qdevBE(f(x),x=-Pi..Pi,500);

> qdevBE(f(x),x=-Pi..Pi,1000);

>

Priklad:

Ilustrujte pomoci animace kmitani struny, jejiz krajni body jsou [0, 0] a [ , 0] a jejiz vychozi tvar pri napnuti popisuje funkce

x/2 pro x z intervalu [0, )

phi(x) = {

pro x z intervalu [ , ]

Nejdrive nadefinujeme funkci phi(x) a jeji liche rozsireni.

> phi:=x->piecewise(x<Pi/2 ,x/2,Pi/2-x/2):

>

> g1:=plot(phi,0..Pi, scaling=constrained, thickness=2,
color=aquamarine):%;

Spocitame Fourierovu radu licheho rozsireni funkce phi(x).

> FS:=GetFourierSeries(phi(x), x=0..Pi, odd);

> FS:=SimplifyCoefficients(FS, simplify);

Nyni potrebujeme rozsirit koeficienty Fourierovy rady vyrazem cos(nt) . Vytvorime tedy novou radu.

> FS2:=Create({ Coefficients(FS)*cos(n*t)}, SinTrigP(n,x));

Pro aproximaci funkce phi(x) pouzijeme prvnich 15 clenu rady.

> s:=GetPartialSum(FS2, 15);

> with(plots):

> g2:=animate(s,x=0..Pi,t=0..2*Pi,frames=50, scaling=constrained):

> display(g1,g2);

Nasledujici serii prikazu zobrazime animaci privnich tri nenulovych harmonickych slozek.

> H1:= GetPartialSum(FS2,1);

> H3:= GetPartialSum(FS2,3)-GetPartialSum(FS2,2);

> H5:= GetPartialSum(FS2,5)-GetPartialSum(FS2,4);

> animate({ H1,H3,H5} ,x=0..Pi,t=0..2*Pi, frames=50);

>

Priklad:

> phi:=x->1/4*x*(Pi-x):

>

> g1:=plot(phi,0..Pi, scaling=constrained, thickness=2,
color=aquamarine):%;

> FS:=GetFourierSeries(phi(x), x=0..Pi, odd);

> FS2:=Create({Coefficients(FS)*cos(n*t)}, SinTrigP(n,x));

> s:=GetPartialSum(FS2, 15);

> with(plots):

> g2:=animate(s,x=0..Pi,t=0..2*Pi,frames=50, scaling=constrained):

> display(g1,g2);

> H1:=GetPartialSum(FS2,1);

> H3:=GetPartialSum(FS2,3)-GetPartialSum(FS2,2);

> H5:=GetPartialSum(FS2,5)-GetPartialSum(FS2,4);

> animate({H1,H3,H5},x=0..Pi,t=0..2*Pi, frames=50);

>

Priklad:

> phi:=x->piecewise(x<Pi/5,x,x<2*Pi/5,2*Pi/5-x,0):

>

> g1:=plot(phi,0..Pi, scaling=constrained, thickness=2,
color=aquamarine):%;

> FS:=GetFourierSeries(phi(x), x=0..Pi, odd);

> FS:=SimplifyCoefficients(FS, simplify);

> FS2:=Create({Coefficients(FS)*cos(n*t)}, SinTrigP(n,x));

> s:=GetPartialSum(FS2, 20):

> with(plots):

> g2:=animate(s,x=0..Pi,t=0..2*Pi,frames=50, scaling=constrained):

> display(g1,g2);

>

Priklad:

Metodou Fourierovych rad reste diferencialni rovnici

Vytvorime obecnou reprezentaci Fourierovy rady a dosadime ji do diferencialni rovnice. Hodnoty koeficientu pak ziskame porovnanim s radou na prave strane.

Nejdrive vytvorime reprezentaci prave strany rovnice. Vyuzijeme procedury Create.

> RHS:=Create({[0], 'general'=[0,1/n^4]}, CosSinTrigP(n,x));

Vytvorime obecnou reprezentaci Fourierovy rady.

> FS:=Create({[a0], 'general'=[aN,bN]}, CosSinTrigP(n,x));

Radu dvakrat derivujeme a k vysledku pricteme dvojnasobek rady samotne. Ziskame tak levou stranu diferencialni rovnice po dosazeni obecne Fourierovy rady.

> tmp:=Derivate(Derivate(FS, 'x'), 'x');

> LHS:=Add(tmp,FS, 1, 2);

Konecne odecteme pravou stranu rovnice od leve.

> S:=Add(LHS,RHS, 1, -1);

Ziskali jsme nekonecnou radu, jejiz vsechny koeficienty musi byt rovny nule. V nasem pripade staci pouze sestrojit a vyresit soustavu rovnic pro nezname aN a bN.

> Coefficients(S);

> {%[1]=0, %[2]=0};

${2*aN-n^2*aN = 0, -1/(n^4)+2*bN-n^2*bN = 0}$

> solve(%, {aN, bN});

${bN = -1/(n^4*(-2+n^2)), aN = 0}$

Konecne muzeme sestrojit hledane reseni.

> Y:=Create({subs(%, bN)}, SinTrigP(n,x));

Spravnost muzeme overit dosazenim reseni do leve strany rovnice.

> Add(Derivate(Derivate(Y, 'x'), 'x'), Y, 1, 2);

> SimplifyCoefficients(%, simplify);

>