Fourier-sarjat

24.11.05 KP3, L/fouriersarjat.mws ja html

>    restart:

>    with(plots):read("e:ns05.mpl");

Warning, the name changecoords has been redefined

Jaksollinen jatkaminen

Jokainen jollain välillä [a,b] määritelty funktio voidaan jatkaa jaksollisena (p=b-a) koko R:ään.

Esim:

>    f:=x->x^2: # välillä [0,1], jatka 1-jaksoisena!

>    Jf:=JJ(f,0..1);

Jf := proc (x::algebraic) local y; y := floor(x); f(x-y) end proc

>    plot(Jf(x),x=-5..5);

[Maple Plot]

>    f:=x->x^2: # välillä [-1,1], jatka 2-jaksoisena!

>    plot(f,-1..1);

>    Jf:=JJ(f,-1..1);

>    plot(Jf,-4..4);

  Kanttiaalto

>    f:=x->piecewise(x>0,1,x<0,-1); plot(f,-1..1);

f := proc (x) options operator, arrow; piecewise(0 < x,1,x < 0,-1) end proc

[Maple Plot]

>    Jf:=JJ(f,-1..1): plot(Jf,-4..4);

[Maple Plot]

>   

Fourier-kertoimet, esimerkkejä Fourier-sarjoista

Fourier-kertoimien kaavat, jotka saatiin ortogonaalisuuden perusteella. Kaavat johti jo kauan ennen Fourieria , kukas muu kuin Euler  v. 1777.

>    f:='f':

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

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

>    an:=(1/Pi)*Int(f(x)*cos(n*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);

bn := 1/Pi*Int(f(x)*sin(n*x),x = -Pi .. Pi)

>    sarja:='a0'+Sum('an'*cos(n*x)+'bn'*sin(n*x),n=1..infinity);

sarja := a0+Sum(an*cos(n*x)+bn*sin(n*x),n = 1 .. infinity)

>    osasumma:=(x,N)->a0+add(an*cos(n*x)+bn*sin(n*x),n=1..N):

Fourier -sarjojen muodostaminen on nyt helppoa. Toki voi tulla integrointivaivaa, mutta Maplea voi käyttää

nautinnollisesti hyödyksi. Voi tietysti syntyä myös tilanteita, joissa on turvauduttava numeeriseen integrointiin.

Siihenkin Maplessa on välineitä.

Esim 1, kanttiaalto

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

f := proc (x) options operator, arrow; piecewise(0 < x,1,x < 0,-1) end proc

>    Jf:=JJ(f,-Pi..Pi); #plot(Jf,-2*Pi..2*Pi);

Jf := proc (x::algebraic) local y; y := floor(1/2*(x+Pi)/Pi); f(x-2*y*Pi) end proc

>   

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

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

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

a0:=0;

>    an:=(1/Pi)*Int('f'(x)*cos(n*x),x=-Pi..Pi); an:=0;

an := 1/Pi*Int(f(x)*cos(n*x),x = -Pi .. Pi)

>   

>    an:=0;

an := 0

>    bn:=(1/Pi)*Int('f'(x)*sin(n*x),x=-Pi..Pi);

bn := 1/Pi*Int(f(x)*sin(n*x),x = -Pi .. Pi)

>    bn:=2*(1/Pi)*int(1*sin(n*x),x=0..Pi);

bn := -2/Pi*(cos(Pi*n)-1)/n

>    bn:=trigsiev(bn,n);

bn := -2/Pi*((-1)^n-1)/n

>    seq(bn,n=1..7);

4/Pi, 0, 4/3/Pi, 0, 4/5/Pi, 0, 4/7/Pi

>    sarja:=Sum('b[n]'*sin(n*x),n=1..infinity);

sarja := Sum(b[n]*sin(n*x),n = 1 .. infinity)

>    osasumma:=(x,N)->add(bn*sin(n*x),n=1..N);

osasumma := proc (x, N) options operator, arrow; add(bn*sin(n*x),n = 1 .. N) end proc

>    osasumma(x,9);

4/Pi*sin(x)+4/3*1/Pi*sin(3*x)+4/5*1/Pi*sin(5*x)+4/7*1/Pi*sin(7*x)+4/9*1/Pi*sin(9*x)

>    sarja:=(4/Pi)*Sum(sin(2*k-1)*x/(2*k-1),k=1..infinity); # Muista iso Sum !

sarja := 4/Pi*Sum(sin(2*k-1)*x/(2*k-1),k = 1 .. infinity)

Kuvaajia:

>    plot([Jf(x),osasumma(x,1)],x=-2*Pi..2*Pi);

[Maple Plot]

>    plot([Jf(x),osasumma(x,1),osasumma(x,3),osasumma(x,5)],x=-2*Pi..2*Pi);

[Maple Plot]

>    plot([Jf(x),seq(osasumma(x,N),N=1..7)],x=-Pi..Pi);

[Maple Plot]

>    plot([Jf(x),seq(osasumma(x,N),N=[1,3,5,31,67])],x=-Pi..Pi,color=[gray,black,green,blue,red]);

[Maple Plot]

>    plot([Jf(x),osasumma(x,101)],x=-0.1..Pi+.1,color=[gray,blue],numpoints=100);

[Maple Plot]

>    display(%,view=[0..Pi,0.7..1.3]);

[Maple Plot]

Miten voi kanttiaallon F-sarja supeta, vaikka kertoimet ovat harmonisen sarjan kertoimet. sin(n*x)-termien on pakko vaikuttaa asiaan. Tutkitaan:

>    g:=(x,k)->sin((2*k-1)*x);

g := proc (x, k) options operator, arrow; sin((2*k-1)*x) end proc

>    plot3d(g(x,k),x=-Pi..Pi,k=1..10):

>    x:=Pi/2:plot([seq([k,g(x,k)],k=1..20)]):

>    kuva:=(x,N)->plot([seq([k,g(x,k)],k=1..N)],style=line);

kuva := proc (x, N) options operator, arrow; plot([seq([k, g(x,k)],k = 1 .. N)],style = line) end proc

>    kuva(0.1,100):

>    display(seq(kuva(x,150),x=[seq(0.1*k,k=0..10*3)]),insequence=true):

>   

Yleinen jakso p=2L

>    restart: with(plots): read("e:ns05.mpl");

Warning, the name changecoords has been redefined

>    a0:=1/(2*L)*Int('f'(x),x=-L..L);

a0 := 1/2*1/L*Int(f(x),x = -L .. L)

>    an:=1/L*Int('f'(x)*cos(n*Pi*x/L),x=-L..L);

an := 1/L*Int(f(x)*cos(n*Pi*x/L),x = -L .. L)

>    bn:=1/L*Int('f'(x)*sin(n*Pi*x/L),x=-L..L);

bn := 1/L*Int(f(x)*sin(n*Pi*x/L),x = -L .. L)

Esim: Puoliaaltotasasuunnattu siniaalto

>    u:= piecewise(t>-L and t < 0,0,t > 0 and t < L,E*sin(omega*t));

u := PIECEWISE([0, -1/2*Pi-t < 0 and t < 0],[E*sin(2*t), -t < 0 and t-1/2*Pi < 0])

>    omega:=2: E:=1:

>    L:=Pi/omega:plot(u(t),t=-1..1,axes=box); omega='omega': E:='E':L:='L':L:=Pi/omega:

>   

[Maple Plot]

>    a0:=1/(2*L)*int(E*sin(omega*t),t=0..L);

a0 := 1/Pi*E

>    an:=1/L*int(E*sin(omega*t)*cos(n*omega*t),t=0..L);  # Huom! n > 1

an := -1/Pi*(cos(n*Pi)+1)*E/(1+n)/(-1+n)

>    an:=trigsiev(an,n);

an := -1/Pi*((-1)^n+1)*E/(-1+n^2)

>    a1:=1/L*int(E*sin(omega*t)*cos(1*omega*t),t=0..L);

a1 := 0

>    bn:=1/L*int(E*sin(omega*t)*sin(n*omega*t),t=0..L);  # Huom! n > 1

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

>    bn:=trigsiev(bn,n);  

bn := 0

>    b1:=1/L*int(E*sin(omega*t)*sin(1*omega*t),t=0..L);

b1 := 1/2*E

>    sarja:=a0+Sum(an*cos(n*omega*t),n=2..infinity)+b1*sin(omega*t);

>   

sarja := 1/Pi*E+Sum(-1/Pi*((-1)^n+1)*E/(-1+n^2)*cos(2*n*t),n = 2 .. infinity)+1/2*E*sin(2*t)

>    osasumma:=unapply(a0+sum(an*cos(n*omega*t),n=2..N)+b1*sin(omega*t),t,N);

osasumma := proc (t, N) options operator, arrow; 2/Pi-1/2*(N+2)/(N+1)/Pi/(-1+(N+1)^2)*(-1)^(N+1)-1/2*(1+2*N)/(N+1)/N/Pi+sum(-2/Pi/(-1+n^2)*(-1)^n*cos(n*t)^2-2/Pi/(-1+n^2)*cos(n*t)^2,n = 2 .. N)+1/2*sin...
osasumma := proc (t, N) options operator, arrow; 2/Pi-1/2*(N+2)/(N+1)/Pi/(-1+(N+1)^2)*(-1)^(N+1)-1/2*(1+2*N)/(N+1)/N/Pi+sum(-2/Pi/(-1+n^2)*(-1)^n*cos(n*t)^2-2/Pi/(-1+n^2)*cos(n*t)^2,n = 2 .. N)+1/2*sin...

>    E:=1: omega:=2:

>    plot([u,osasumma(t,10)],t=-L..L);

[Maple Plot]

>    plot([seq(osasumma(t,N),N=2..10)],t=-L..L);

[Maple Plot]

>    plot([osasumma(t,50)],t=-L..L,numpoints=100,axes=boxed);

[Maple Plot]

>