Fourier-sarjojen perustyöarkki

22.11.01  HA

                                                                                                  14.8.02

Apufunktioita: JJ, trigsiev, jana

>    restart:

Warning, the name changecoords has been redefined

>    with(plots):

>    #read("/p/edu/mat-1.414/maple/fourier.mpl"): # muuta v302:ksi

>    #read("c:\\opetus\\maple\\ohjelmat.mpl"):

>    read("/home/apiola/opetus/peruskurssi/v2-3/maple/v302.mpl"):

Määritellään uudet koodit varmuuden vuoksi myös tässä, vaikka ovatkin  fourier.mpl:ssä (ja ohjelmat.mpl:ssä).

>    # Jaksollinen jatko: MapleTech Vol 3 No. 3 1996
#
JJ:=proc(f,d::range)
  subs({'F'=f,'L'=lhs(d),'D'=rhs(d)-lhs(d)},
  proc(x::algebraic) local y;
    y:=floor((x-L)/D);
    F(x-y*D);
  end)
end:

# Esim:
#sw:=JJ(signum,-1..1);
#plot(sw,-4..4);
#
 parillinenjatko:=f->unapply(piecewise(x>0,f(x),x<0,f(-x)),x):
 paritonjatko:=f->unapply(piecewise(x>0,f(x),x<0,-f(-x)),x):

>    trigsiev:=proc(lauseke,n)
local tulos;
tulos:=subs(sin(n*Pi)=0,lauseke);
tulos:=subs(cos(n*Pi)=(-1)^n,tulos);
tulos:=algsubs(exp(-2*I*n*Pi)=1,tulos);
tulos:=algsubs(exp(2*I*n*Pi)=1,tulos);
tulos:=subs((-1)^(2*n)=1,tulos);simplify(tulos);
end:

>    jana:=(x,y)->plot([[x,0],[x,y]],thickness=3,color=blue):

>    #display([seq(jana(n,abs(an)),n=1..6)],axes=box);

>    op(trigsiev):

Fourier-kertoimet

Reaalinen muoto, cos ja sin

Nämä on syytä ottaa pohjaksi työskentelyyn ja olla valmis modifioimaan esim. integrointirajoja,

integrointia paloittain ym. (Kokemus opettaa välttämään liian pitkälle vietyä automatiikkaa.)

>    f:='f': T:='T': omega:='omega': #omega:=2*Pi/T:

>    a0:=(2/T)*int(f(t),t=-T/2..T/2):
an:=(2/T)*int(f(t)*cos(n*omega*t),t=-T/2..T/2):
bn:=(2/T)*int(f(t)*sin(n*omega*t),t=-T/2..T/2):

>    #sarja:=a0/2+sum(an*cos(n*omega*t)+bn*sin(n*omega*t),n=1..infinity):

>    osasumma:=N->a0/2+sum(an*cos(n*omega*t)+bn*sin(n*omega*t),n=1..N):

>   

Sinisarja ja kosinisarja

Jäljempänä esitämme parillisen ja parittoman laajennusoperaattorin. Käytännössä parillista/paritonta laajennusta vastaa kosini/sinisarja, jossa lähdetään liikkeelle välillä [0,d] määritellystä  funktiosta f, joka laajennetaan parilliseksi/parittomaksi. Tällöin yllä olevat kaavat tulevat muotoon:

1) kosinisarja:

>    f:='f': T:='T': omega:='omega': #omega:=2*Pi/T:

>    a0:=(1/d)*int(f(t),t=-d..d);an:=(1/d)*int(f(t)*cos(n*omega*t),t=-d..d);bn:=0;

>   

a0 := 1/d*int(f(t),t = -d .. d)

an := 1/d*int(f(t)*cos(n*omega*t),t = -d .. d)

bn := 0

Koska yllä olevat integroitavat ovat parilllisia, voidaan integroida [0..d], kunhan tulos kerrotaan 2:lla.

Tällöin tullaan toimeen alkuperäisellä, välillä [0,d] määritellyllä funktiolla, eikä tarvitse kohdistaa integrointia parilliseen laajennukseen. Siispä kosinisarjakaavat ovat:

>    f:='f':T:='T':omega:='omega':#omega:=2*Pi/T:
a0:=(2/d)*int(f(t),t=0..d);
an:=(2/d)*int(f(t)*cos(n*omega*t),t=0..d);bn:=0;
sarja:='a0'/2+sum('an'*cos(n*omega*t),n=1..infinity);

a0 := 2/d*int(f(t),t = 0 .. d)

an := 2/d*int(f(t)*cos(n*omega*t),t = 0 .. d)

bn := 0

sarja := 1/2*a0+sum(an*cos(n*omega*t),n = 1 .. infinity)

Sinisarja saadaan aivan vastaavasti:

>    f:='f':T:='T':omega:='omega':#omega:=2*Pi/T:
bn:=(2/d)*int(f(t)*sin(n*omega*t),t=0..d);
sarja:=sum('bn'*cos(n*omega*t),n=1..infinity);

bn := 2/d*int(f(t)*sin(n*omega*t),t = 0 .. d)

sarja := sum(bn*cos(n*omega*t),n = 1 .. infinity)

Kompleksinen muoto

Kokeillaan vaihteeksi sellaista tyyliä, että c on n:n funktio, tässä se on mukavaa, kun kaavatkin ovat yhtenäiset.  c-funktio on syytä määritellä unapply-tavalla. Laskentatehossa nimittäin on dramaattinen ero. (unapply evaluoi (esim. integroi) määrittelyaikana, kun taas -> määrittely tekee sen kutsuaikana.)

>    T:='T': f:='f': omega:='omega': omega:=2*Pi/T:

>    c:=unapply(1/T*int(f(t)*exp(-I*n*omega*t),t=0..T),n);c(0):=1/T*int(f(t),t=0..T);

c := proc (n) options operator, arrow; 1/T*int(f(t)*exp(-2*I*n*Pi/T*t),t = 0 .. T) end proc

c(0) := 1/T*int(f(t),t = 0 .. T)

Vaikka seuraava toimii, niin ei kannata tehdä, kts. yllä.

>    #c:=n->1/T*int(f(t)*exp(-I*n*omega*t),t=0..T);

>    # sarja:=sum('c(n)'*exp(I*n*'omega'*t),n=-infinity..infinity);

Esim.

Muodosta alla määriteltävän funktion Fourier-sarja, piirrä jaksollinen laajennus ja Fourier-osasummien kuvaajia. Piirrä myös taajuusspektrit.  Käsittele sekä reaalimuotoista sin/cos-sarjaa että kompleksimuotoista.

>    f:=t->piecewise(-Pi<t and t < 0,t+Pi,0<t and t<Pi,Pi);

f := proc (t) options operator, arrow; piecewise(-Pi < t and t < 0,t+Pi,0 < t and t < Pi,Pi) end proc

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

[Maple Plot]

>    jaksf:=JJ(f,-Pi..Pi):

>    plot(jaksf,-2*Pi..2*Pi);

[Maple Plot]

>    T:=2*Pi:omega:=2*Pi/T:

>    a0:=(2/T)*int(f(t),t=-Pi..Pi);
an:=(2/T)*int(f(t)*cos(n*omega*t),t=-Pi..Pi):an:=trigsiev(an,n);
bn:=(2/T)*int(f(t)*sin(n*omega*t),t=-Pi..Pi):bn:=trigsiev(bn,n);

a0 := 3/2*Pi

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

bn := -(-1)^n/n

Koska halusimme välttyä assume -komennon aiheuttamilta sivuvaikutuksilta, sievensimme omalla trigsiev- komennolla. Huomaamme, että int   osasi käsitellä piecewise -määriteltyä funktiota ihan nätisti. Luultavaa on, että se osaa myös epäjatkuvalle, mutta kannattaa tarkkailla.

>    #sarja:=a0/2+sum(an*cos(n*omega*t)+bn*sin(n*omega*t),n=1..infinity);

>    osasumma:=N->a0/2+sum(an*cos(n*omega*t)+bn*sin(n*omega*t),n=1..N);

osasumma := proc (N) options operator, arrow; 1/2*a0+sum(an*cos(n*omega*t)+bn*sin(n*omega*t),n = 1 .. N) end proc

>    osasumma(5);

3/4*Pi+2/Pi*cos(t)+sin(t)-1/2*sin(2*t)+2/9*1/Pi*cos(3*t)+1/3*sin(3*t)-1/4*sin(4*t)+2/25*1/Pi*cos(5*t)+1/5*sin(5*t)

Tässä tyylissä (jossa an ja bn integroitiin symbolisesti valmiiksi), on aivan sopivaa (ja suositeltavaa) käyttää -> funktiomääritystapaa.

>    plot([jaksf(t),osasumma(10)],t=-2*Pi..2*Pi);

[Maple Plot]

>    plot([jaksf(t),seq(osasumma(k),k=1..10)],t=-2*Pi..2*Pi);

[Maple Plot]

Taajuusspekrit:

Piirtäminen hoituu oikein näppärästi yllä määrittelemällämme jana -funktiolla.

>    A0:=a0/2; An:=abs(an-I*bn);

A0 := 3/4*Pi

An := abs(-1/Pi*((-1)^n-1)/n^2+(-1)^n/n*I)

>    phi0:=0: phin:=argument(an-I*bn);

phin := argument(-1/Pi*((-1)^n-1)/n^2+(-1)^n/n*I)

>    display([jana(0,A0),seq(jana(n,An),n=1..6)],axes=box,title="Amplitudispektri");display([jana(0,phi0),seq(jana(n,phin),n=1..6)],title="Vaihespektri");

[Maple Plot]

[Maple Plot]

Kompleksinen muoto

>    c:=unapply(1/T*int(f(t)*exp(-I*n*omega*t),t=0..T),n);c(0):=1/T*int(f(t),t=0..T);

>    sarja:=sum('c(n)'*exp(I*n*'omega'*t),n=-infinity..infinity);

c := proc (n) options operator, arrow; -1/2*I*(exp(n*Pi*I)-1)/n/exp(n*Pi*I) end proc

c(0) := 1/2*Pi

sarja := sum(c(n)*exp(n*omega*t*I),n = -infinity .. infinity)

>    c(0);

1/2*Pi

>    c(1);

-I

>    seq(c(n),n=-5..5);

1/5*I, 0, 1/3*I, 0, I, 1/2*Pi, -I, 0, -1/3*I, 0, -1/5*I

>    display([seq(jana(n,abs(c(n))),n=-6..6)],axes=box,title="Amplitudispektri");

[Maple Plot]

>    display([seq(jana(n,argument(c(n))),n=-6..6)],title="Vaihespektri");

[Maple Plot]

>   

Parilliset ja parittomat laajennukset

Määritellään laajennusoperaattorit " parillinenjatko"  ja " paritonjatko ", jotka jatkavat
positiivisella R-akselilla määritellyn funktion parillisena/parittomana koko R:ään.
(Ovat ohjelmat.mpl:ssä)

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

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

>    f:=t->t;

f := proc (t) options operator, arrow; t end proc

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

[Maple Plot]

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

[Maple Plot]

>    plot(paritonjatko(x->x^2),-1..1);

[Maple Plot]

>    jakspariton:=JJ(paritonjatko(x->x^2),-1..1):

>    plot(jakspariton,-3..3);

[Maple Plot]

Onpas meillä nyt hyvät välineet. Tehtävien ratkominen on pelkkää nautiskelua!

Lisäksi kannattaa ottaa vielä numeerisen integroinnin tapaus.

Fourier-sarjan laskeminen numeerisesti integroiden

Titenkään integrointi ei aina onnistu symbolisesti.  Otetaan esimerkki, jossa samalla tulee sinisarjan kaavan käyttö.

Kehitettävä sinisarjaksi välillä [-1,1] funktio:

>    f:=t->t^sin(4*Pi*t);

f := proc (t) options operator, arrow; t^sin(4*Pi*t) end proc

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

[Maple Plot]

>    d:=1:omega:=Pi:

>    b:=array(1..20):

>    for n from 1 to 20 do
   b[n]:=(2/d)*evalf(Int(f(t)*sin(n*omega*t),t=0..d));
od: n:='n':

Tässä yksi sitkeä Maple "feature" : for-loopin jälkeen n:llä on arvo 21, mikä johtaa jatkossa virhetilanteeseen, siksi n vapautetaan tässä.

>    osasumma:=N->sum(b[n]*sin(n*omega*t),n=1..N);

osasumma := proc (N) options operator, arrow; sum(b[n]*sin(n*omega*t),n = 1 .. N) end proc

>    osasumma(10);

1.494333420*sin(Pi*t)+.1753057331*sin(2*Pi*t)+.1309937131*sin(3*Pi*t)-.7743672430*sin(4*Pi*t)-.8964208140e-1*sin(5*Pi*t)+.7947576352e-1*sin(6*Pi*t)+.3068905782*sin(7*Pi*t)+.3346763050e-1*sin(8*Pi*t)+.3...
1.494333420*sin(Pi*t)+.1753057331*sin(2*Pi*t)+.1309937131*sin(3*Pi*t)-.7743672430*sin(4*Pi*t)-.8964208140e-1*sin(5*Pi*t)+.7947576352e-1*sin(6*Pi*t)+.3068905782*sin(7*Pi*t)+.3346763050e-1*sin(8*Pi*t)+.3...
1.494333420*sin(Pi*t)+.1753057331*sin(2*Pi*t)+.1309937131*sin(3*Pi*t)-.7743672430*sin(4*Pi*t)-.8964208140e-1*sin(5*Pi*t)+.7947576352e-1*sin(6*Pi*t)+.3068905782*sin(7*Pi*t)+.3346763050e-1*sin(8*Pi*t)+.3...

>    plot([paritonjatko(f)(t),op(map(osasumma,[3,5,7,17]))],t=-1..1);

[Maple Plot]

>   

>   

>   

Tehtiin pikku briljeeraus vielä mappaamalla osasumma erinäisille indekseille.