Propagation eines Wellenpakets mit Dispersion

Parameter :

(* Halbe Breite und Zentrumsprosition des Wellenpakets zur Zeit t=0  *)
σ0=0.2;
x0=0;
 (* Phasen- und Gruppengeschwindigkeit (Dispersionsfreie Welle) *)
c=2;      
(* Breite des Rechenfensters (aus praktischen Gründen abhängig von σ0) *)
w=50*σ0;
(* Rechengenauigkeit: In wie viele Punkte (im Ortsraum) soll die x-Achse zerlegt werden? *)
n=100;
(* Dispersionsrelation *)
ω[k_]:=c*(k)^1.03
(* ============================================================ *)
(* Das Folgende sind abgeleitete Größen.                       *)
(* Das Sampling Intervall im Ortsraum ergibt sich aus w und n *)
(* Die maximale Frequenz im k-Raum ist die sogenannte          *)
(* "Nyquist-Frequenz". Sie ist gerade die Hälfte der Sampling- *) 
(* Frequenz n/w. Dies folgt aus dem "Abtast-Theorem" das besagt *)
(* dass ein Signal mindestens mit der doppelten Frequenz        *) 
(* abgetatstet werden muss. Sind in dem Orginal-Signal         *)
(* Frequenzen größer als die halbe Abtastfrequenz enthalten, so *)
(* gehen diese Komponenten bereits beim Abtasten verloren.     *)
(* Bei der folgenden Fourier-Zerlegung "lohnt" es sich daher   *)
(* nicht, nach Freqenzen größer als der Nyquist-Frequenz zu    *)
(* suchen.                                                      *)
dx=w/n;
fnyq=0.5/dx;

Dies ist unser Ausgangswellenpaket im Orts- (x-) Raum. Es kann z.B. ein Wellenberg auf einer Wasseroberfläche sein. Oder eine lokale Auslenkung einer gespannten Saite. Hier kann man zwischen einer Gauß-Kurve oder einer Kastenfunktion wählen (vrgl. Aufgabe 2, Blatt 2) Nur eine der beiden Zellen ausführen!

ψ0[x_]:=PDF[NormalDistribution[x0,σ0],x]

ψ0[x_]:=UnitStep[x0+σ0-x]*UnitStep[x-(x0-σ0)]

Plot[ψ0[x], {x,-w/2,w/2},PlotRange->All, Axes->False, Frame->True, PlotStyle->{Thick},Filling->Axis]

Nun müssen wir bestimmen, aus welchen Schwingungen (im Ortsraum, wir reden noch nicht von der Zeit) ψ0 besteht. Wir müssen es also in ein Spektrum von harmonischen Schwingungen mit Ortsfrequenzen k zerlegen. Das ist mathematisch eine Fouriertransformation. Wir approximieren das Integral numerisch durch eine endliche Summe.

ϕ[k_]=Sum[ψ0[x]*Exp[-I*k*x]*dx, {x,-w/2,w/2, dx}]/Sqrt[2*π];

Das quadrierte ϕ ist das "Spektrum" von ψ0 im inversen Orts- oder k-Raum. Da die Fouriertarnsformation auf einem diskretisierten (Orts-) Signal von n Datenpunkten basiert, die eine Breite w im Ortsraum überspannen, kann das Fouriertarsnformierte Signal keine Frequenzen größer als die Nyquist-Frequenz fnyq = n/2w enthalten.

Plot[Abs[ϕ[k]]^2,{k, -2*π*fnyq,2*π*fnyq},PlotRange->All, Axes->False, Frame->True, PlotStyle->{Thick},Filling->Axis]

Nun bauen wir aus dem bekannten Spektrum das zeitabhängige Wellenpaket ψ(x,t) als (prinzipiell unendliche, aber hier numerisch diskretisierte) Summe aus ebenen Wellen auf. 
Dazu wird jede Frequenzkomponente von ϕ mit einem Faktor Exp(-iωt) multipliziert, wobei ω(k)  gerade die Dispersionsrelation der Welle ist, also die Funktion, die Orts - und Zeitfrequenz miteinander verbindet. Anschließend wird das ganze Fourier-rücktransformiert.
 Später, in der Quantenmechanik werden wir den Faktor Exp(-iωt) treffenderweiseals "Zeitentwicklungsoperator" bezeichnen.

ψ[x_,t_]=Sum[ϕ[k]*Exp[-I*ω[k]*t]*Exp[I*k*x]*2*π*fnyq/n, {k,-2*π*fnyq,2*π*fnyq,2*π*fnyq/n}]/Sqrt[2*π];

Nun können wir das zeitentwickelte Wellenpaket plotten. Durch Rundungsfehler enthält die aus dem Zeitentwickelten Spektrum rücktransformierte Funktion kleine imaginäre Anteile. Diese "bügeln" wir weg, indem wir nur den Realteil betrachten.

Plot[{ψ0[x],Re[ψ[x,1]]},{x,-w/2,w/2}, PlotRange->All, Axes->False, Frame->True, PlotStyle->{Thick},Filling->Axis]

Oder als hübsche Animation ...

ListAnimate[Table[Plot[{ψ0[x],Re[ψ[x,time]]},{x,-w/2,w/2}, PlotRange->All, Axes->False, Frame->True, PlotStyle->{Thick},Filling->Axis], {time, 0, 2.5, 0.33}], AnimationRunning->False]