Fourierreeksen

Versie 2, 9 oktober 2000.
Dit Mathematica notebook was ooit bedoeld voor wiskunde 2 voor chemici.
Bij het opstellen is gebruik gemaakt van ideeën uit het notebook blok.nb van W. van der Kallen.
Dit is versie 2, waarin de ervaringen verwerkt zijn met versie 1 tijdens toenmalige werkcolleges.
\

Dit notebook is bedoeld om mee te experimenteren. Het is verstandig eerst een werkkopie te maken en de oorspronkelijke versie onveranderd te bewaren.

Voorbereidende opdrachten:

Needs["Graphics`Graphics`"]; Off[NIntegrate::ploss];
Off[NIntegrate::ncvb];
Off[NIntegrate::inum];
Off[NIntegrate::slwcon];

Periodieke functies maken

De functies met periode 2π in de opgaven 10.1.7-9 zijn alle vastgelegd door tussen -π en π iets voor te schrijven en dat periodiek voort te zetten.
Dat mechanisme kunnen we ook in Mathematica gebruiken. We zetten het algemener op: periode 2L, waarbij L positief moet zijn.

Clear[f, voorbeeld]
voorbeeld[x_] := x^2
L = Pi;
f[x_] := 
      Module[{xx}, xx = x - 2  L Round[x/(2  L)]; 
        voorbeeld[xx]]

In de definitie van f vervang je x door xx tussen -L en L, waarbij xx en x een, geschikt gekozen, geheel getal is, verkregen als  Round[x/(2L)].
Hier hebben we de functie de functie in opgave 10.1.8 gemaakt. De grafiek ziet er zo uit:

Plot[f[x], {x, -10, 10}]

Als je een andere functie met periode 2π wilt definiëren, behoef je alleen de definitie van voorbeeld te veranderen. Voor een andere periode moet je ook L aanpassen.

Opgave F1. Plot een grafiek voor de functie in 1.10.7 en ook voor de functie z in E49.

Fourierreeksen van even periodieke functies

Een functie f wordt even genoemd als f(-x)=f(x); de grafiek is dan symmetrisch ten opzichte van de verticale as.

Opgave F2. Welke van de periodieke functies bekeken in de opgaven 10.1.7-9 en E49  zijn even?

We bekijken nu een eenvoudige even functie met periode 2π:

Clear[f, voorbeeld]
sprongpunt = Pi/3;
voorbeeld[x_] := If[Abs[x] ≤ sprongpunt, 1, 0]
L = Pi;
f[x_]:=Module[{xx}, xx = x - 2  L Round[x/(2  L)]; 
        voorbeeld[xx]]

fplot = 
      Plot[f[x], {x, -6, 6}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[ .005]}}]

Opgave F3. Wat is de wiskundige betekenis van de verticale rode lijnen?

We berekenen voor de huidige betekenis van de functie f de 
Fouriercoëfficiënten:
	(Subscript[a, 0]) = (1/(2  L)) (∫_(-L)^(L ))f(
            x) dx

	(Subscript[a, n ]) = (1/L) (∫_(-L)^L) f(
            x) cos(π n x/ L) dx   voor n≥1
Ofwel:

Clear[a]
a[0] = 
        1/(2 L) 
          NIntegrate[f[x], {x, -L, L}];
a[n_] := a[n] = (1/L) 
          NIntegrate[f[x] Cos[Pi n x/L], {x, -L, L}]

We gebruiken NIntegrate om de integralen numeriek te berekenen.
Die integralen willen we niet vaker dan één keer berekenen. Daar zorgt de constructie a[n_]:=a[n]=...  voor.
Laten we de Fouriercoëfficiënten tot n=25 bekijken:

ListPlot[Table[{n, a[n]}, {n, 0, 25}], PlotRange -> All, 
      PlotStyle -> {RGBColor[1, 0, 0], PointSize[ .02]}]

De berekening duurt even. Er moeten 25 integralen berekend worden.

Met de opdracht ?a kun je bekijken welke informatie Mathematica nu omtrent a heeft.

Nu proberen we de functie f te reconstrueren uit de Fourierreeks 
  	(∑_(n = 0)^∞) (Subscript[a, n]) cos (πnx/L)
Die oneindige reeks is de limiet van de beginstukken. Het beginstuk tot en met n=grens beschrijven we als volgt:

Clear[cosreeks]
cosreeks[x_, grens_] := cosreeks[x, grens] = 
        Sum[a[n] Cos[Pi n x/L], {n, 0, grens}]

Clear[cosplot]
cosplot[grens_] := cosplot[grens] = 
        Plot[cosreeks[x, grens], {x, -6, 6}, 
          PlotStyle -> {{RGBColor[0, 1, 0], 
                Thickness[ .005]}}]

Als je dit combineert met de grafiek van de functie f zie je dat het beginstuk van de cosinusreeks tot en met n=25 de functie redelijk benadert:

Show[fplot, cosplot[25]]

De volgende opdrachten vertonen een serie van benaderingen met beginstukken van de cosinusreeks:

Clear[plottot]
plottot[grens_, opties___] := (approx[x_] = cosreeks[x, grens]; 
        Plot[{f[x], approx[x]}, {x, -6, 6}, 
          opties])
Table[
        Show[GraphicsArray[
            plottot[n, PlotPoints -> 200, PlotRange -> {- .2, 1.2}, 
              PlotStyle -> {{Hue[0.99], Thickness[0.005]}, {Hue[0.1], 
                    Thickness[0.005]}}]]], {n, 0,  16}];

Sluit de reeks grafieken door te twee keer te klikken op de op één na meest links gelegen blauwe lijn ernaast. Nu blijft alleen de begingrafiek over. Klik daarin twee keer om een filmpje te starten.
Je stopt het filmpje door nog een keer te klikken.

De conclusie lijkt te zijn dat de cosinusreeks inderdaad de oorspronkelijke functie als limiet heeft. Een uitzondering treedt op aan de randpunten. Daar schiet de benadering door. Dat heet het Gibbsverschijnsel.
Om dat nader te bekijken nemen we een een benadering met langere beginstukken:

plot40 = 
      Plot[{f[x], cosreeks[x, 40]}, {x, -2, 2}, 
        PlotPoints -> 200]

plot60 = 
      Plot[{f[x], cosreeks[x, 60]}, {x, -2, 2}, 
        PlotPoints -> 200]

Een uitvergroting:

Show[plot60, PlotRange -> {{0, 1.2}, { .5, 1.2}}]

Andere voorbeelden

Om al deze stappen uit te voeren voor andere even periodieke functies, hoeft alleen in de definitie van f de functie voorbeeld veranderd te worden. Daarna kun je de rest van de opdrachten opnieuw uitvoeren.

Voor het gemak van de gebruiker, en om de resultaten voor verschillende functies te kunnen vergelijken, herhalen we hieronder de opdrachten.

Opgave F4. Kies zelf een even periodieke continue functie. (Continu betekent: geen sprongen.)
a) Voer hiervoor de bovenstaande opdrachten uit.
b) Gaan de Fouriercoëfficiënten sneller naar 0 in vergelijking met het eerste voorbeeld? 
c) Treedt er een Gibbsverschijnsel op?

Clear[f, voorbeeld]
voorbeeld[x_] := vulhierzelfwatin
L = Pi;
f[x_]:=Module[{xx}, xx = x - 2  L Round[x/(2  L)]; 
        voorbeeld[xx]]

fplot = 
      Plot[f[x], {x, -6, 6}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[ .005]}}]

Berekening Fouriercoëfficiënten:

Clear[a]
a[0] = 
        1/(2 L) 
          NIntegrate[f[x], {x, -L, L}];
a[n_] := a[n] = (1/L) 
          NIntegrate[f[x] Cos[Pi n x/L], {x, -L, L}]

ListPlot[Table[{n, a[n]}, {n, 0, 25}], PlotRange -> All, 
      PlotStyle -> {RGBColor[1, 0, 0], PointSize[ .02]}]

Beginstukken:

Clear[cosreeks]
cosreeks[x_, grens_] := cosreeks[x, grens] = 
        Sum[a[n] Cos[Pi n x/L], {n, 0, grens}]

Clear[cosplot]
cosplot[grens_] := cosplot[grens] = 
        Plot[cosreeks[x, grens], {x, -6, 6}, 
          PlotStyle -> {{RGBColor[0, 1, 0], 
                Thickness[ .005]}}]

Gecombineerde grafiek:

Show[fplot, cosplot[25]]

Filmpje van een aantal benaderingen:

Clear[plottot]
plottot[grens_, opties___] := (approx[x_] = cosreeks[x, grens]; 
        Plot[{f[x], approx[x]}, {x, -6, 6}, 
          opties])
Table[
        Show[GraphicsArray[
            plottot[n, PlotPoints -> 200, PlotRange -> {- .2, 1.2}, 
              PlotStyle -> {{Hue[0.99], Thickness[0.005]}, {Hue[0.1], 
                    Thickness[0.005]}}]]], {n, 0,  16}];

Misschien moet je hierboven de PlotRange veranderen om de hele grafiek te zien te krijgen.

Fourierreeksen van algemene functies

Als de periodieke functies niet even zijn, heb je ook de sinustermen nodig. In 10.3 vind je de Fouriercoëfficiënten (Subscript[a, n]) en (Subscript[b, n]), met (Subscript[a, n]) als hierboven en voor n⩾1
	(Subscript[b, n ])= (1/L)(∫_(-L)^(L  ))f(x)sin(πnx/L) dx
en de Fourierreeks
	(Subscript[a, 0])+ (∑_(n = 1)^∞) ( (Subscript[a, n]) 
                  cos (πnx/L)  +  (Subscript[b, n ]) 
                  sin (πnx/L) )

Opgave F5. Waarom zijn alle  (Subscript[b, n ])  gelijk aan 0  als f een even functie is?

Clear[f, voorbeeld]
voorbeeld[x_] := If[x < 0, 0, x]
L = 1;
f[x_]:=Module[{xx}, xx = x - 2  L Round[x/(2  L)]; 
        voorbeeld[xx]]

fplot = 
      Plot[f[x], {x, -1.7, 1.7}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[ .005]}}]

Berekening Fouriercoëfficiënten:

Clear[a, b]
a[0] = 
        1/(2 L) 
          NIntegrate[f[x], {x, -L, L}];
a[n_] := a[n] = (1/L) 
          NIntegrate[
            f[x] Cos[Pi n x/L], {x, -L, 
              L}]
b[n_] := b[n] = (1/L) 
          NIntegrate[f[x] Sin[Pi n x/L], {x, -L, L}]

a[0]
a[1]
b[1]

Opgave F6. Controleer de waarde van   (Subscript[a, 0]) door de integraal ook met de hand uit te rekenen.

alijst = 
      ListPlot[Table[{n, a[n]}, {n, 0, 25}], PlotRange -> All, 
        PlotStyle -> {RGBColor[1, 0, 0], 
            PointSize[ .02]}]
blijst = 
      ListPlot[Table[{n, b[n]}, {n, 1, 25}], PlotRange -> All, 
        PlotStyle -> {RGBColor[0, 0, 1], 
            PointSize[ .02]}]
Show[alijst, blijst]

Beginstukken:

Clear[fourreeks]
fourreeks[x_, grens_] := fourreeks[x, grens] = 
        a[0] + Sum[
            a[n] Cos[Pi n x/L] + b[n] Sin[Pi n x/L], {n, 1, 
              grens}]

Clear[plotreeks]
plotreeks[grens_] := plotreeks[grens] = 
        Plot[fourreeks[x, grens], {x, -1.7, 1.7}, 
          PlotStyle -> {{RGBColor[0, 1, 0], 
                Thickness[ .005]}}]

Gecombineerde grafiek:

Show[fplot, plotreeks[10]]

Filmpje van een aantal benaderingen:

Clear[plottot]
plottot[grens_, opties___] := (approx[x_] = fourreeks[x, grens]; 
        Plot[{f[x], approx[x]}, {x, -1.7, 1.7}, 
          opties])
Table[
        Show[GraphicsArray[
            plottot[n, PlotPoints -> 200, PlotRange -> {- .2, 1.2}, 
              PlotStyle -> {{Hue[0.99], Thickness[0.005]}, {Hue[0.1], 
                    Thickness[0.005]}}]]], {n, 0,  16}];

Ook hier treedt bij het sprongpunt het Gibbsverschijnsel op.

Opgave F7. Onderzoek op dezelfde manier een andere periodieke functie.

Opgave F8. Aan welke eigenschap moet een functie voldoen willen alle (Subscript[a, n])  gelijk aan 0 zijn?
Zo'n periodieke functie kan dus in een sinusreeks ontwikkeld worden.

Opgave F9. Wat gebeurt er als je een beginstuk van een Fourierreeks weglaat? Probeer het eerst te beredeneren, en probeer het dan uit.
Met een beginstuk bedoelen we bijvoorbeeld:
 	(Subscript[a, 0]) +  (∑_(n = 1)^2)( 
            (Subscript[a, n])cos(πnx/L)+ (Subscript[b, n ])sin(πnx/L) ).
Wat gebeurt er als je de cosinustermen weglaat?

Warmtevergelijking

Als een toepassing van Fourierreeksen bekijken we de eendimensionale warmtevergelijking 
	(∂_t) T(x, t)  =  (c^2)(∂_(x, 2))T(x, t)
zoals ook besproken wordt in Kreyszig 11.5.
T(x,t) is de temperatuur op plaats x op tijdstip t in een goed ingepakte staaf. De constante c is afhankelijk van de soortelijke warmte en het warmtegeleidingsvermogen van de staaf. Laten we voor het gemak c=1 nemen.

Opgave F10. We hebben de differentiaalvergelijking in Mathematica notatie opgeschreven. Hoe ziet deze vergelijking eruit in "gewone" notatie?

We veronderstellen dat de staaf lengte L=4 heeft en dat de plaats x voldoet aan 2≤x≤2.
Op het begintijdstip t=0 wordt de temperatuur beschreven door
	T(x,0) = 15+Abs[Abs[x]-1]

Opgave F11. Maak een grafiek van de temperatuurverdeling op tijdstip t=0.

We berekenen de cosinus-Fouriercoëfficiënten van deze even functie T(x,0), en noemen die a[n,0]:

Clear[a, T]
L = 2;
T[x_, 0] := 15 + Abs[Abs[x] - 1]
a[0, 0] = 
        1/(2 L) 
          NIntegrate[T[x, 0], {x, -L, L}];
a[n_, 
        0] := a[n, 0] = (1/L) 
          NIntegrate[T[x, 0] Cos[Pi n x/L], {x, -L, L}]

ListPlot[Table[{n, a[n, 0]}, {n, 0, 30}], 
      PlotRange -> All, 
      PlotStyle -> {RGBColor[1, 0, 0], PointSize[ .02]}]

En we geven een bijbehorende benadering van de temperatuursverdeling op tijdstip t=0:

Plot[Sum[a[n, 0] Cos[Pi n x/L], {n, 0, 30}], {x, -2, 2}, 
      PlotStyle -> {{RGBColor[0, 1, 0], Thickness[ .005]}}, 
      PlotRange -> {{-2, 2}, {14.9, 16.1}}]

Niet alleen voor t=0, maar voor alle t≥0 is er een ontwikkeling 
   T(x,t) = (Subscript[a, 0])(t) + (∑_(n = 1)^∞) (Subscript[a, n])(t) cos(πnx/2)
 geldig voor -2≤x≤2.
 Invullen in de differentiaalvergelijking levert door vergelijking van de voorfactoren van de cosinussen:
 	(Subscript[a, 0])'(t)  =  0
 	(Subscript[a, n])'(t)  = -( π(n/2))^2 (Subscript[a, n])(t)

Opgave F12. Controleer dat.

We vinden dus  
	(Subscript[a, 0])(t) = (Subscript[a, 0])(0)
	(Subscript[a, n])(t) = (Subscript[a, n])(0) (e^(-((t(πn/2))^2)))

Opgave F13. Waarom?

Dat vertellen we Mathematica:

a[n_, t_] := a[n, t] = 
        a[n, 0] Exp[-t (Pi n/2)^2]

Nu kunnen we op verschillende momenten de temperatuursverdeling schetsen:

Clear[teken]
teken[t_] := 
      Plot[Sum[a[n, t] Cos[Pi n x/L], {n, 0, 30}], {x, -2, 2}, 
        PlotStyle -> {{RGBColor[0, 1, 0], Thickness[ .005]}}, 
        PlotRange -> {{-2, 2}, {14.9, 16.1}}]

teken[0.3]

En een filmpje:

Clear[plott]
plott[t_, 
        opties___] := (approx[x_] = 
          Sum[a[n, t] Cos[Pi n x/L], {n, 0, 30}]; 
        Plot[approx[x], {x, -2, 2}, opties])
Table[
        Show[GraphicsArray[
            plott[n* .01, PlotPoints -> 200, 
              PlotRange -> {{-2, 2}, {14.9, 16.1}}, 
              PlotStyle -> {{Hue[0.99], Thickness[0.005]}}]]], {n, 0,  
          50}];

Hier wordt de temperatuursverdeling gegeven voor t=0, t=.01, t=0.2, ..., t=0.3.

Opgave F14. Voer hetzelfde uit voor een andere temperatuursverdeling op tijdstip t=0.