Mathe für Physiker III Wintersemester 02/03 Walter Oevel 25. 11. 2002

Ü b u n g s b l a t t 6

Abgabe von ∗ Aufgaben bis zum 2.12.2002, 11^{oo} Uhr, im Zettelkasten auf dem D1.

Aufgabe 23*: (Die DFT als Lösung einer Interpolationsaufgabe. 10 Bonuspunkte) Gegeben sei die L–periodische Funktion f(x). Durch äquidistantes ”Sampling“ im Abstand L/N werden die Werte f_{j} = f(x_{j}) an den Punkten x_{j} = x_{0} + j · L/N (j ∈ ℤ) gemessen. Betrachte die Fourier–Approximationen

S_{N}(x) = ∑_{k=−N/2}^{N/2−1} d_{k} · e^{i·k·2·π·(x−x_{0})/L}

für gerades N bzw.

S_{N}(x) = ∑_{k=−(N−1)/2}^{(N−1)/2} d_{k} · e^{i·k·2·π·(x−x_{0})/L}

für ungerades N. Hierbei seien die d_{k} die diskreten Fourier–Koeffizienten des Datensatzes (f_{0}, …, f_{N−1}), welche mittels d_{N+k} = d_{k} für alle Frequenzen k ∈ ℤ definiert seien. Zeige, dass die Ausgangsfunktion f(x) in folgendem Sinne von S_{N}(x) approximiert wird:

An den Stellen x_{j} gilt S_{N}(x_{j}) = f(x_{j}).

Anleitung: Schreibe die Summen für S_{N} in die Form ∑_{k=0}^{N−1} um.

Musterlösung:

S_{N}(x) = ∑_{k=−N/2}^{N/2−1} d_{k} · e^{i·k·2·π·(x−x_{0})/L} = ∑_{k=0}^{N/2−1} d_{k} · e^{i·k·2·π·(x−x_{0})/L} + ∑_{k=−N/2}^{−1} d_{k} · e^{i·k·2·π·(x−x_{0})/L}

=^{(k′=k+N)} ∑_{k=0}^{N/2−1} d_{k} · e^{i·k·2·π·(x−x_{0})/L} + ∑_{k′=N/2}^{N−1} \underbrace{d_{k′−N}}_{=d_{k′}} ·e^{i·(k′−N)·2·π·(x−x_{0})/L}

=^{(k′=k)} ∑_{k=0}^{N/2−1} d_{k} · e^{i·k·2·π·(x−x_{0})/L} + e^{−i·2·π·N·(x−x_{0})/L} · ∑_{k=N/2}^{N−1} d_{k} · e^{i·k·2·π·(x−x_{0})/L}.

S_{N}(x_{j}) = ∑_{k=0}^{N/2−1} d_{k} · e^{i·k·2·π·j/N} + \underbrace{e^{−i·2·π·j}}_{=1} · ∑_{k=N/2}^{N−1} d_{k} · e^{i·k·2·π·j/N} = ∑_{k=0}^{N−1} d_{k} · e^{i·j·k·2·π/N}.

Mit d_{N+k} = d_{k} gilt für gerades N:

An den Stellen x = x_{j} = x_{0} + j · L/N ergibt sich

Dies ist aber genau die Formel der inversen DFT, also S_{N}(x_{j}) = f_{j} = f(x_{j}).

Für ungerades N geht das analog:

S_{N}(x) = ∑_{k=−(N−1)/2}^{(N−1)/2} d_{k} ·e^{i·k·2·π·(x−x_{0})/L} = ∑_{k=0}^{(N−1)/2} d_{k} ·e^{i·k·2·π·(x−x_{0})/L} + ∑_{k=−(N−1)/2}^{−1} d_{k} ·e^{i·k·2·π·(x−x_{0})/L}

=^{(k′=k+N)} ∑_{k=0}^{(N−1)/2} d_{k} · e^{i·k·2·π·(x−x_{0})/L} + ∑_{k′=(N−1)/2+1}^{N−1} \underbrace{d_{k′−N}}_{=d_{k′}} ·e^{i·(k′−N)·2·π·(x−x_{0})/L}

=^{(k′=k)} ∑_{k=0}^{(N−1)/2−1} d_{k} · e^{i·k·2·π·(x−x_{0})/L} + e^{−i·2·π·N·(x−x_{0})/L} · ∑_{k=(N−1)/2+1}^{N−1} d_{k} · e^{i·k·2·π·(x−x_{0})/L}.

An den Stellen x = x_{j} = x_{0} + j · L/N ergibt sich

S_{N}(x_{j}) = ∑_{k=0}^{(N−1)/2} d_{k} · e^{i·k·2·π·j/N} + \underbrace{e^{−i·2·π·j}}_{=1} · ∑_{k=(N−1)/2+1}^{N−1} d_{k} · e^{i·k·2·π·j/N} = ∑_{k=0}^{N−1} d_{k} · e^{i·j·k·2·π/N}.

Dies ist aber genau die Formel der inversen DFT, also S_{N}(x_{j}) = f_{j} = f(x_{j}).

Aufgabe 24*: (FFT mit MuPAD. 10 Bonuspunkte) Lies die Hilfeseite zu MuPADs numeric::fft. Betrachte die Funktion

f(x) = |x| = \frac{π}{2} − \frac{4}{π} · (\frac{cos(x)}{1^{2}} + \frac{cos(3 · x)}{3^{2}} + \frac{cos(5 · x)}{5^{2}} + ⋯)

auf dem Intervall [−π, π]. Wähle N = 16. Lasse MuPAD die diskreten Fourier–Koeffizienten zu f_{j} = f(−π + j · 2 · π/N) berechnen. Stelle f(x) und die in Aufgabe 23 betrachtete Approximation S_{N}(x) in einer gemeinsamen Graphik dar.

Relevante MuPAD–Funktion: abs, numeric::fft, map, float, Re, plotfunc2d.

Achtung: FFT_{MuPAD} = N · FFT_{Skript}! Achtung: Die numerische Auswertung enthält durch Rundungsfehler kleine Imaginärteile. Der Plotter plotfunc2d ignoriert nicht-reelle Eingabedaten, d.h., die Imaginärteile müssen vor dem Plotten entfernt werden!

Musterlösung:

>> x0:= -PI:              // Linkes Intervallende
>> L:= 2*PI:              // Periode
>> f:= x -> abs(x):       // die Funktion auf dem Intervall [x0, x0 + L]
>> N:= 16:                // Anzahl der Stuetzstellen
>> ff:= [f(x0 + j*L/N) $ j = 0..N-1]: // Sampling
>> d:= numeric::fft(ff):  // DFT der Daten ff

Um an Aufgabe 23 anknüpfen zu können, muss die FFT MuPADs auf die FFT des Skripts umgebogen werden:

>> d:= map(d, x -> x/N):

Die Rekonstruktion nach Aufgabe 23 braucht die diskreten Fourier–Koeffizienten d_{k} mit Frequenzen k von −N/2 bis N/2 − 1. In der MuPAD–Liste d = [d_{0}, …, d_{N−1}] müssen die negativen Frequenzen k daher per d_{k} = d_{N+k} addressiert werden. Weiterhin ist zu beachten, dass MuPAD-Listen mit 1 beginnend indiziert werden, also d_{k} = d[1 + k]:

>> S:= x -> _plus(d[1 + N + k]*exp(I*k*(x-x0)) $ k=-N/2 .. -1)
          + _plus(d[1 +   k  ]*exp(I*k*(x-x0)) $ k = 0 .. N/2 - 1):

Man beachte, dass ein Aufruf von S durch numerische Rundungsfehler kleine Imaginärteile enthält:

>> S(0.3)
               0.2605531595 + 1.137599129e-16 I

Daher wird nur der Realteil von S(x) an plotfunc2d übergeben:

>> plotfunc2d(f(x), Re(S(x)), x = -PI..PI, Grid = 200):

Aufgabe 25*: (Eine graphische FFT–Spielerei. 10 Bonuspunkte) Zu Abschnitt 1.3.11 (”Digitales Filtern“) des Skripts. Die im Internet erhältliche Datei DieKleine.txt enthält ein MuPAD–Array A = array(1..64, 1..32) mit reellen Daten, die als Grauwerte eines verrauschten Bildes zu interpretieren sind. Das Bild kann in einer MuPAD–Sitzung mittels

>> read("D:\\Dokumente und Einstellungen\\meinPfad\\DieKleine.txt"):
>> plot(plot::matrixplot(A))

eingelesen und visualisiert werden (wenn die Datei lokal unter dem angegebenen Pfad gespeichert wurde). Mit FourierData:= numeric::fft(A) kann die 2-dimensionale DFT von A durchgeführt werden. Die Filter–Routine

>> filter:= proc(eps, FourierData)
   local M,N, i, j;
   begin
     M:= op(FourierData, [0, 2, 2]):
     N:= op(FourierData, [0, 3, 2]):
     for i from 1 to M do
       for j from 1 to N do
         if abs(i/M - 1/2) < eps or abs(j/N - 1/2) < eps then
             FourierData[i,j]:= 0;
         end_if:

       end_for:
     end_for:
     return(FourierData);
   end_proc:

wird in der Datei DieKleine.txt mitgeliefert. Sie entfernt aus einer Datenmatrix alle Daten, die in der ”Mitte“ des Datenfeldes sitzen (bei diskreten Fourier–Daten entspricht dies hochfrequenten Anteilen). Finde interaktiv einen Wert des Parameters eps ∈ (0, 1), bei dem das im Bild versteckte ”Nutzsignal“ sichtbar wird: es handelt sich um den Namen und den Geburtstag meiner jüngsten Tochter. Anleitung: die Fourier–Rücktransformation der gefilterten Fourier–Daten ist ein komplexer Datensatz. Man visualiere mittels plot::matrixplot den Realteil oder den Absolutbetrag des gefilterten Datensatzes.

Musterlösung: Folgende Schleife erzeugt gefilterte Bilder, in denen die Werte eps = 0.1, 0.2, …, 0.9 etc. ausprobiert werden:

>> read("DieKleine.txt"):
>> FourierData:= numeric::fft(A):
>> for eps in [j*0.1 $ j = 1..9] do
     FilteredFourierData:= filter(eps, FourierData):
     FilteredData:= map(numeric::invfft(FilteredFourierData), abs):
     plot(plot::matrixplot(FilteredData)):
  end_for:

Das beste Ergebnis zeigt sich für eps = 0.3:

                                                                  eps = 0.2

                                                              eps = 0.3

                                                                  eps = 0.4

Aufgabe 26*: (Die Wellengleichung. 10 Bonuspunkte) Betrachte eine schwingende Saite, die an den Punkten x = 0 und x = L eingespannt ist. Ihre Auslenkung wird durch u(x, t) beschrieben, wobei u(0, t) = u(L, t) = 0 ∀t ∈ ℝ gilt. Die zeitliche Dynamik wird durch die partielle DGL (”Wellengleichung“)

\frac{∂u}{∂t^{2}} = c^{2} · \frac{∂^{2}u}{∂x^{2}}

beschrieben, wobei c > 0 eine Materialkonstante ist.

a) Setze für die Lösung u(x, t) = ∑_{k=1}^{∞} b_{k}(t) · sin(k · π · x/L) mit (zeitabhängigen) Fourier– Koeffizienten b_{k}(t) an. Die Randbedingungen u(0, t) = u(L, t) = 0 sind damit automatisch erfüllt. Setze den Fourier–Ansatz in die Wellengleichung ein. Für die einzelnen Fourier–Koeffizienten b_{k}(t) ergibt sich jeweils eine einfach zu lösende gewöhnliche DGL. Zeige, dass die Wellengleichung durch

u(x, t) = ∑_{k=1}^{∞} b_{k}(0) · cos(c · k · π · t/L) · sin(k · π · x/L)

+ ∑_{k=1}^{∞} \frac{L · ḃ_{k}(0)}{c · k · π} · sin(c · k · π · t/L) · sin(k · π · x/L)

mit beliebigen Konstanten b_{k}(0), ḃ_{k}(0) gelöst wird.

b) Gegeben seien Anfangsbedingungen u(x, 0) = u_{0}(x) und \frac{∂u}{∂t}(x, 0) = u̇_{0}(x) mit 0 ≤ x ≤ L zum Zeitpunkt t = 0. Wie lassen sich die auf dem Intervall [0, L] gegebenen Funktionen u_{0}(x) und u̇_{0}(x) in Sinus–Reihen entwickeln? Wie ergibt sich mit a) diejenige Lösung der Wellengleichung, die zum Zeitpunkt t = 0 diese Anfangsbedingungen erfüllt?

Von der schwingenden Saite werden periodische Anregungen α_{k} · cos(c · k · π · t/L) an die umgebende Luft abgegeben, die dann als vom Instrument erzeugte Töne wahrgenommen werden. Die zu k = 1 gehörende Schwingung niedrigster Frequenz c · π/L heißt ”Grundton“, die zu k > 1 gehörenden Schwingungen ”Obertöne“ oder auch ”Harmonische Oberschwingungen“. Der Klang eines Instruments wird durch das Amplitudenverhältnis der Obertöne zum Grundton bestimmt.

c) Zeige: der Klang einer Gitarre ändert sich mit dem Ort x_{0} ∈ (0, L), an dem die Saiten angezupft werden. Betrachte dazu die spezielle Anfangsbedingungen

u_{0}(x) = \{ \begin{matrix}\frac{x}{x_{0}} · U_{0} & für 0 ≤ x ≤ x_{0}, \\ \frac{L − x}{L − x_{0}} · U_{0} & für x_{0} ≤ x ≤ L \end{matrix}

und u̇_{0}(x) = 0. Berechne die Fourier–Koeffizienten von u_{0}(x) in Abhängigkeit von x_{0} und stelle die entsprechende Lösung der Wellengleichung über ihre Fourier–Reihe dar. Betrachte für die Fälle x_{0} = L/2 und x_{0} = L/3 die Amplitudenverhältnisse von Grundton und Obertönen, die von einem Punkt x der Saite ausgesandt werden.

d) Zeige: bei einer E-Gitarre hängt der Klang weiterhin von der Plazierung der ”Pick-Ups“ ab. (Unter den Saiten sind an einer Stelle x_{p} ∈ (0, L) kleine Spulen (= pick-ups) angebracht, in denen die schwingende Stahlsaite eine elektrische Spannung U(t) proportional zur Auslenkung u(x_{p}, t) induziert.) Diskutiere wiederum das Amplitudenverhältnis der Obertöne zum Grundton.

Musterlösung:

a) Für den Fourier–Ansatz werden nur sin–Terme gewählt, da dadurch die Randbedingungen

u(0, t) = ∑_{k=1}^{∞} b_{k}(t) · sin(k · 0) = 0, u(L, t) = ∑_{k=1}^{∞} b_{k}(t) · sin(k · π) = 0

automatisch erfüllt sind. Einsetzen des Fourier–Ansatzes in die Wellengleichung liefert

\frac{∂^{2}}{∂t^{2}} u(x, t) = \frac{∂^{2}}{∂t^{2}} ∑_{k=1}^{∞} b_{k}(t) · sin(k · π · x/L) = ∑_{k=1}^{∞} (\frac{d^{2}}{dt^{2}} b_{k}(t)) · sin(k · π · x/L)

= c^{2} · \frac{∂^{2}}{∂x^{2}} u(x, t) = c^{2} · \frac{∂^{2}}{∂x^{2}} ∑_{k=1}^{∞} b_{k}(t) · sin(k · π · x/L)

= c^{2} · ∑_{k=1}^{∞} b_{k}(t) · \frac{∂^{2}}{∂x^{2}} sin(k · π · x/L) = ∑_{k=1}^{∞} b_{k}(t) · ( − \frac{c^{2} · k^{2} · π^{2}}{L^{2}} ) · sin(k · π · x/L).

Vergleich der einzelnen Fourier–Moden in

∑_{k=1}^{∞} (\frac{d^{2}}{dt^{2}} b_{k}(t)) · sin(k · π · x/L) = ∑_{k=1}^{∞} ( − \frac{c^{2} · π^{2} · k^{2}}{L^{2}} · b_{k}(t)) · sin(k · ·π · x/L)

liefert die gewöhnliche DGL

\frac{d^{2}}{dt^{2}} b_{k}(t) = −c^{2} · \frac{π^{2}}{L^{2}} · k^{2} · b_{k}(t)

mit der allgemeinen Lösung

b_{k}(t) = b_{k}(0) · cos(c · k · π · t/L) + \frac{L · ḃ_{k}(0)}{c · k · π} · sin(c · k · π · t/L).

Nach Konstruktion erfüllt damit

u(x, t) = ∑_{k=1}^{∞} b_{k}(0) · cos(c · k · π · t/L) · sin(k · π · x/L)

+ ∑_{k=1}^{∞} \frac{L · ḃ_{k}(0)}{c · k · π} · sin(c · k · π · t/L) · sin(k · π · x/L)

für jede beliebige Wahl der Integrationskonstanten b_{k}(0), ḃ_{k}(0) die Wellengleichung und die Randbedingungen u(0, t) = u(L, 0) = 0.

b) Setzt man in der obigen Lösungsformel t = 0, so ergibt sich

u(x, 0) = ∑_{k=1}^{∞} b_{k}(0) · sin(k · π · x/L), \frac{∂u}{∂t}(x, 0) = ∑_{k=1}^{∞} ḃ_{k}(0) · sin(k · π · x/L).

Damit stellen sich die Integrationskonstanten b_{k}(0) in der Tat als die Fourier–Koeffizienten von u(x, 0) = u_{0}(x) heraus, während die ḃ_{k}(0) die Fourier–Koeffizienten von \frac{∂u}{∂t}(x, 0) sind.

Frage: läßt sich jede Funktion u_{0}(x) mit x ∈ [0, L] in eine Sinus–Reihe entwickeln (wo bleiben die Cosinus–Terme)? Antwort: Die Sinus–Terme reichen! Setze u_{0}(x) (definiert auf [0, L]) \underline{ungerade} auf das Interval [−L, L] fort. Die übliche Fourier–Entwicklung auf [−L, L] liefert dann eine Sinus–Reihe. Analog ist die Funktion u̇(x) ungerade auf [−L, L] fortzusetzen und in eine übliche Fourier–Sinus– Reihe auf [−L, L] zu entwickeln.

0) Gegeben u_{0}(x) und u̇_{0}(x) mit x ∈ [0, L].

1) Berechne die Fourier–Koeffizienten von u_{0}(x) und u̇_{0}(x) über die Formel für ungerade 2 · L– periodische Funktionen über dem Intervall [−L, L]:

b_{k}(0) = \frac{2}{L} · ∫_{0}^{L} u_{0}(x) · sin(k · π · x/L) dx, k ∈ ℕ,

ḃ_{k}(0) = \frac{2}{L} · ∫_{0}^{L} u̇_{0}(x) · sin(k · π · x/L) dx, k ∈ ℕ.

2) Die Lösung der Wellengleichung ist gegeben durch

u(x, t) = ∑_{k=1}^{∞} b_{k}(0) · cos(c · k · π · t/L) · sin(k · π · x/L)

+ ∑_{k=1}^{∞} \frac{L · ḃ_{k}(0)}{c · k · π} · sin(c · k · π · t/L) · sin(k · π · x/L)

mit den in 1) berechneten Werten b_{k}(0) und ḃ_{k}(0).

c) Wir ziehen die in b) aufgestellten Schritte 1), 2) durch. Die Fourier–Koeffizienten ḃ_{k}(0) von u̇_{0}(x) = 0 sind alle 0. Die Fourier–Koeffizienten b_{k}(0) von u_{0}(x) werden mit MuPAD 3.0 ermittelt:

>> assume(L_0 > 0):
>> assume(0 < x0 < L0):
>> assume(L > 0):
>> assume(k <> 0):
>> b[k]:= 2/L*int(x/x0*U0*sin(k*PI*x/L), x = 0 .. x0) +
          2/L*int((L-x)/(L-x0)*U0*sin(k*PI*x/L), x = x0 .. L):
>> assume(k, Type::Integer):
>> b[k]:= simplify(b[k])
                      2       / PI k x0 \
                   2 L  U0 sin| ------- |
                              \    L    /
                   ----------------------
                       2  2
                     PI  k  x0 (L - x0)

Anmerkung: in MuPAD 2.5 reagiert der Integrierer nicht richtig auf die per assume gesetzten Eigenschaften. Ein workaround besteht darin, dem Integrierer durch die Option Continuous mitzuteilen, dass er sich nicht um potentielle Unstetigkeiten der Stammfunktion kümmern soll:

>> b[k]:= 2/L*int(x/x0*U0*sin(k*PI*x/L), x = 0 .. x0, Continuous) +
          2/L*int((L-x)/(L-x0)*U0*sin(k*PI*x/L), x = x0 .. L, Continuous):
>> assume(k, Type::Integer):
>> b[k]:= simplify(b[k])
                        2       / k PI x0 \
                     2 L  U0 sin| ------- |
                                \    L    /
                  - ------------------------
                     2   2   2      2   2
                    k  PI  x0  - L k  PI  x0

Hiermit ergibt sich nach b) die Lösung u(x, t):

u(x, t) = \frac{2 · U_{0}}{π^{2} · \frac{x_{0}}{L} · (1 − \frac{x_{0}}{L} )} · ∑_{k=1}^{∞} \frac{sin(k · π · x_{0}/L)}{k^{2}} · cos(c · k · π · t/L) · sin(k · π · x/L).

Sei

a_{k} = \frac{2 · U_{0}}{π^{2} · \frac{x_{0}}{L} · (1 − \frac{x_{0}}{L} )} · \frac{sin(k · π · x_{0}/L)}{k^{2}} · sin(k · π · x/L)

die Amplitude der k-ten Oberschwingung, die vom Ort x ausgesandt wird. Das Amplitudenverhältnis der k-ten Harmonischen zum Grundton ist damit

\frac{a_{k}}{a_{1}} = \frac{1}{k^{2}} · \frac{sin(k · π · x_{0}/L)}{sin(π · x_{0}/L)} · \frac{sin(k · π · x_{0}/L)}{sin(π · x_{0}/L)}.

\underline{Für x_{0} = L/2:}

\frac{a_{k}}{a_{1}} = \frac{1}{k^{2}} · \frac{sin(k · π/2)}{sin(π/2)} · \frac{sin(k · π · x/L)}{sin(π · x/L)} = \frac{sin(k · π/2)}{k^{2}} · \frac{sin(k · π · x/L)}{sin(π · x/L)}

= \{ \begin{matrix}0 & für gerades k, \\ \frac{(−1)^{n}}{k^{2}} · \frac{sin(k · π · x/L)}{sin(π · x/L)} & für ungerades k = 2 · n + 1. \end{matrix}

x_{0} = L/2: [figure: k; |a_{k}|; 1; 2; 3; 4; 5; 6; 7] Alle geraden Harmonischen verschwinden!

\underline{Für x_{0} = L/3:}

\frac{a_{k}}{a_{1}} = \frac{1}{k^{2}} · \frac{sin(k · π/3)}{sin(π/3)} · \frac{sin(k · π · x/L)}{sin(π · x/L)} = \{ \begin{matrix}0 & für k = 3 · n, n ∈ ℕ, \\ \frac{(−1)^{n}}{k^{2}} · \frac{sin(k · π · x/L)}{sin(π · x/L)} & für k = 3 · n + 1, n ∈ ℕ_{0}, \\ \frac{(−1)^{n}}{k^{2}} · \frac{sin(k · π · x/L)}{sin(π · x/L)} & für k = 3 · n + 2, n ∈ ℕ_{0}. \end{matrix}

x_{0} = L/3: [figure: k; |a_{k}|; 1; 2; 3; 4; 5; 6; 7] Jede dritte Harmonische verschwindet!

d) Am Punkt x = x_{p} ist die Auslenkung der Seite

u(x_{p}, t) ∼ ∑_{k=1}^{∞} \underbrace{(\frac{sin(k · π · x_{0}/L) · sin(k · π · x_{p}/L)}{k^{2}})}_{Amplitude der Zeitschwingung} · cos(c · k · π · t/L).

Die Abhängigkeit des Klangs von x_{p} ist damit die selbe wie die Abhängigkeit vom ”Zupfpunkt“ x_{0}.