Regge-Wheeler-ova rovnice  
Přepis do prvního řádu v čase

(∂_tt) ψ==(∂_xx)ψ-V ψ//HoldForm
(∂_t)MatrixForm[ {ψ,χ} ]==MatrixForm[ {χ,(∂_xx)ψ-V ψ} ]//HoldForm

želví souřadnice

rStar==r+2M Log[r/(2M)-1]
Solve[%,r][[1]]
rInv[rStar_]=Evaluate@(r/.%/.M->1)
rInv[-10.]

Konkrétní potenciál pro numerické řešení

params = {M->1, s->2,l->2}
(Subscript[V, RW])[r_]=(1-2M/r)(l(l+1)/r^2+2M/r^3 (1-s^2))/.params

Row[{ Plot[ (Subscript[V, RW])[r],{r,2,10},ImageSize->Medium,PlotTheme->"Scientific"],Plot[ (Subscript[V, RW])[rInv[x]],{x,-8,10},ImageSize->Medium,PlotTheme->"Scientific"] }]

Pravá strana rovnice 
(∂_t)(({{ψ}, {χ}}))==(({{χ}, {(∂_xx)ψ-V ψ}})) 
Numerický integrátor pro MOL

rhs := Compile[ {{x,_Real,1},{V,_Real,1},{u,_Real,2}},Module[ { i,j,f,dx,n=Length[u],u0,iψ = 1,iχ = 2},
dx = x[[2]]-x[[1]];
u0 ={0.,0.};
f=u ;
f[[1]] = {u[[1,iχ]],0};
f[[-1]] = {u[[n,iχ]],0};
For[i=2, i≤n-1,i++,
f[[i,1]] =u[[i,2]];
f[[i,2]] = (u[[i-1,1]]-2*u[[i,1]]+u[[i+1,1]])/(dx^2)-V[[i]] u[[i,1]];
];
f
]
(*,CompilationOptions->{"InlineExternalDefinitions"->True}*)];

(*<<"CompiledFunctionTools`"
CompilePrint[rhs]*)

(* Heun's RK3, no explicit time in rhs *)
RK3Step[u_,dt_] :=Block[ {u1,u2,u3,dt3},
dt3=dt/3.;
u1 = u+dt3 rhs[xDomain,Vx, u];
u2 = u+2dt3 rhs[xDomain,Vx,u1];
u3 = u1+dt rhs[xDomain,Vx,u2];
0.75 u3+0.25 u
]

Počáteční prostorový tvar vlny 
Parametry gridu
Předpočítané hodnoty potenciálu

ψ0[x_Real]:=If[x>3&&x<6,{Sin[Pi x/3]^2,0},{0,0}]

nDomain = 5000;
xBdry = 100.;
xDomain = Range[ -xBdry,xBdry,2xBdry/nDomain];
Vx = (Subscript[V, RW])[rInv[#]]&/@ xDomain;
uInit = ψ0/@ xDomain;
iObserver = 11/20*nDomain
rObserver==rInv[xDomain[[iObserver]]]

Inicializace
Hlavní smyčka řešení pomocí MOL

tFin = 94;
dtPrint = 0.5;

cfl = 0.5;
t0 = 0; 
dx0 = xDomain[[2]]-xDomain[[1]];
dt0 = cfl * dx0
u0 = uInit;

tPrint = t0;
results = {};observed={};
(*u0
u0 = RK3Step[ xDomain, u0,dt0]
*)
Do[
If[ t0≥tPrint,AppendTo[ results,{t0,u0}]; tPrint+=dtPrint];
AppendTo[ observed,{t0,u0[[iObserver]]}];
u0 = RK3Step[  u0,dt0];
t0+=dt0;
,{iStep,1,Round[tFin/dt0]}]

Jaké ψ vidí pozorovatel na fixním r

ω022 = 0.37367-0.08896 I;
ω122 = 0.34671-0.27391 I; (*podle https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5253841/  *)
Show[
ListLinePlot[Transpose@{observed[[;;,1]], observed[[;;,2,1]]},PlotTheme->"Scientific",PlotRange->{-0.21,0.21},ImageSize->Large],
Plot[ Re[ 0.41 Exp[-I ω022 (t-0.9)]],{t,15,80},PlotRange->Full]
]

Problém: Ukažte, že až na amplitudu a fázový posun je pro různá počáteční data (např. různý tvar, poloha nebo počet vlnek ) tzv. “ringdown”  stejný (jde o dominantní QNM)

Animace časové závislosti prostorového průběhu ψ

Table[
ListLinePlot[ Transpose@{xDomain,results[[k,2,;;,1]]} ,PlotRange->{-1,1},PlotTheme->"Scientific",PlotLabel->"t="<>ToString@results[[k,1]],ImageSize->666],
 {k,1,Length@results}
];
ListAnimate@%