1+1 Vlnová rovnice na úsečce

bázi tvoří sinusovky s uzly na koncích

□[f_] := Laplacian[f,{x,y,z}]-D[f,t,t]

fk=Sin[ k x ]Cos[k t]
□[%]

Plot[ Evaluate@Table[Sin[k x](6-k)^2,{k,1,5}],{x,0,Pi},PlotTheme->"Detailed"]

fce0 = Exp[-8(x-2)^2]
yMax = Table[Abs@fce0,{x,0,Pi,0.01}]//Max

Vezmeme počáteční data
 u(x,t=0)=f0(x)
Přepokládáme nulovou časovou derivaci v t=0, proto cos(kt) výše.
Spočteme rozklad do Fourierovy řady

nMax=12;
ak = Table[ (2/Pi) NIntegrate[ fce0  Sin[k x],{x,0,Pi}],{k,1,nMax}]

Sečteme ∑ (Subscript[a, k]) sin(k x) cos(k t)
a ověříme, že

tmp = ak . Table[ fk,{k,1,nMax}]  ;
kontrola==□[tmp]
Plot[ {tmp/.t->0,fce0},{x,0,Pi},PlotTheme->"Detailed",PlotLegends->None,AspectRatio->1/3]

Animate[Plot[ tmp/.t->T,{x,0,Pi},PlotRange->{-yMax,yMax},PlotTheme->"Detailed",PlotLegends->None,AspectRatio->1/3,ImageSize->700],{T,0,2Pi,Pi/80},DisplayAllSteps->True]

2+1 Vlnová rovnice na kruhu

radiální bázi tvoří Besselovy funkce s uzly na okraji kruhu
úhlovou bázi tvoří siny a kosiny

BesselJ[n, k R]Cos[n ϕ]
Laplacian[ BesselJ[n-1, k R]Cos[(n-1) ϕ], {R,ϕ,z},"Cylindrical"]//FullSimplify

Table[BesselJ[n,x],{n,0,4}];
Plot[%,{x,0,10},PlotTheme->"Detailed"]

Známe-li řešení, vybereme několik nejnižších módů jak konečnou bázi

f[n_Integer,m_Integer]:=BesselJ[n-1, R BesselJZero[n-1,m] ]Cos[(n-1) ϕ]
ω[n_Integer,m_Integer]:=BesselJZero[n-1,m]

baze = Table[ {n->n0,m->m0},{n0,1,7},{m0,1,7}];
baze = Flatten[baze,1]

Malování bázových funkcí (chvíli trvá)

cys={R->Sqrt[x^2+y^2],ϕ->ArcTan[x,y]};(*převod kartézských do polárních souřadnic*)
barvicky =(Blend[{Yellow,Red,Darker@Red,Black,Darker@Blue,Blue,Cyan},-#+1/2]&);
uvnitrKruhu=Function[{x,y},x^2+y^2<1];
bazoveFunkce = Table[ f[n,m]/.q,{q,baze}];
nahodnyVyber=Take[DeleteDuplicates@RandomInteger[{1,Length@baze},100],UpTo[9] ]//Sort;
ParallelTable[ 
Rasterize@ContourPlot[f/.cys,{x,-1,1},{y,-1,1},
RegionFunction->uvnitrKruhu,PlotRange->Full,MaxRecursion->0,PlotPoints->40,
ColorFunction->barvicky,ColorFunctionScaling->False,ImageSize->480,Contours->40,ContourStyle->None],
{f,bazoveFunkce[[nahodnyVyber]]}]

Počáteční data

f0 =Exp[-8((R Cos[ϕ]-0.6)^2+(R Sin[ϕ])^2)](1-R^2) (* faktor (1-R^2) zaručuje splnění okr. podmínek *)

grafF0=ContourPlot[f0/.cys,{x,-1,1},{y,-1,1},RegionFunction->uvnitrKruhu,PlotRange->Full,MaxRecursion->0,PlotPoints->40,ColorFunction->barvicky,ColorFunctionScaling->False,ImageSize->280,Contours->40,ContourStyle->None];
Plot3D[f0/.cys,{x,-1,1},{y,-1,1},RegionFunction->uvnitrKruhu,PlotRange->{-1,1},MaxRecursion->0,PlotPoints->40,ColorFunction->(barvicky[#3]&),ColorFunctionScaling->False,
MeshStyle->Gray,ImageSize->380];
Column[{Style["Počáteční data u(x,y,t=0)",{20}], Row[{%%,%}]},Alignment->Center]

Nyní je třeba počáteční data zprojektovat na bázi
dot[. , n,m] představuje  skalární součin s (Subscript[f, nm])

dot[h_,n_Integer,m_Integer]:=NIntegrate[h  f[n,m] R,{R,0,1},{ϕ,0, 2Pi},MaxRecursion->0, Method->{"CartesianRule",Method->{{"GaussKronrodRule","Points"->16},{"TrapezoidalRule","Points"->12}}}]

ParallelEvaluate[Off[NIntegrate::ncvb]];

PrintTemporary["Počítám počáteční hodnoty amplitud ..."];
amplitudy = ParallelTable[ dot[f0,n,m]/dot[f[n,m],n,m] /.q,{q,baze}]

Nyní máme řešení jako funkci R,ϕ a t

funkcet= Table[ f[n,m] Cos[ ω[n,m] t] /.q,{q,baze}].amplitudy//N;
funkce0= funkcet/.t->0//N;

Zkontrolujeme, že řeší vlnovou rovnici

Laplacian[funkcet,{R,ϕ,z},"Cylindrical"]-D[funkcet,t,t];
%/.{R->0.7,ϕ->0,t->0}//N

Vizuálně zkontrolujeme, že báze stačí k popisu počátečních dat
(malá báze by nestačila a obě funkce by se lišily)

ContourPlot[funkce0/.cys,{x,-1,1},{y,-1,1},RegionFunction->uvnitrKruhu,PlotRange->Full,MaxRecursion->0,PlotPoints->40,ColorFunction->barvicky,ColorFunctionScaling->False,ImageSize->280,Contours->40,ContourStyle->None];
Row[%,grafF0]

Animace vývoje pole 
(trvá to ...)

SetSharedVariable[progress];progress=Graphics[];
Monitor[
frames3D = ParallelTable[ progress=
Rasterize@Plot3D[funkcet/.cys,{x,-1,1},{y,-1,1},RegionFunction->uvnitrKruhu,PlotRange->{-1,1},MaxRecursion->0,PlotPoints->40,
ColorFunction->(barvicky[#3]&),ColorFunctionScaling->False,ImageSize->580],{t,0,5,1/20}],
Show[progress,ImageSize->333]];
ListAnimate[frames3D,DisplayAllSteps->True]

Ještě jedna animace

SetSharedVariable[progress];progress=Graphics[];
Monitor[ 
frames =ParallelTable[
 progress=
Rasterize@ContourPlot[funkcet/.cys,{x,-1,1},{y,-1,1},RegionFunction->uvnitrKruhu,PlotRange->Full,MaxRecursion->0,PlotPoints->40,ColorFunction->barvicky,ColorFunctionScaling->False,ImageSize->480,Contours->Range[-1,1,1/60],ContourStyle->None],
{t,0,5,1/20}],

Show[progress,ImageSize->333]
]

ListAnimate[frames,DisplayAllSteps->True]

Grupová a fázová rychlost

(rychle a názorně a bez Fourierových transformací )

Animate[
Plot[ Sin[10x -t]/(1+((7x-t)/12)^4),{x,-10,10},PlotRange->{-1,1},PlotTheme->"Scientific",PlotPoints->200],{t,-50,50,0.5},DisplayAllSteps->True]

Vlnová rovnice s disperzí
Vlnový balík s obálkou

(Subscript[□, p])[f_] := Laplacian[f,{x,y,z}]-(1/(c^2))D[f,t,t]-(((Subscript[ω, p])^2)/(c^2))f
g[x-(Subscript[v, g]) t]Exp[I k(x- (Subscript[v, f]) t)]
rovnice = (Subscript[□, p])[%]//FullSimplify

Předpokládáme, že k je velké
členy u (k^2) a k musejí vymizet vhodnou volbou rychlostí (Subscript[v, f]) a (Subscript[v, g])

Solve[ (k^2) ((c^2)-(Subscript[v, f]^2))+(Subscript[ω, p]^2)==0 &&(c^2)-(Subscript[v, f]) (Subscript[v, g])==0,{(Subscript[v, f]) ,(Subscript[v, g])}][[2]]//Simplify

Ještě ukážeme, že (Subscript[v, g])= ((dω)/dk)

D[(Sqrt[(c^2) (k^2)+(Subscript[ω, p]^2)]),k]
%/.Solve[(c^2) (k^2)+(Subscript[ω, p]^2)==ω^2,k][[2]]//Simplify