Vlnová rovnice  
Přepis do prvního řádu v čase

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

Numerický integrátor pro MOL

rhs = Compile[ {{x,_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,iψ]] =u[[i,iχ]];
f[[i,iχ]] = (u[[i-1,iψ]]-2*u[[i,iψ]]+u[[i+1,iψ]])/(dx^2);
];
f]
,CompilationOptions->{"InlineExternalDefinitions"->True}];
(*<<"CompiledFunctionTools`"
CompilePrint[rhs]*)

MOL == Method of lines
1. nejprve zavedu konečnědim. prostor funkcí proměnné x
2. dostanu velkou sadu ODR
3. tu řeším spolehlivou num. metodou  

Zde numerický integrátor RK3

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

g0[x_Real]:=If[x>0.5&&x<0.75,{Sin[4Pi x],0},{0,0}]

nDomain = 300;
xDomain = Range[ 0.,1,1/nDomain];
uInit = g0/@ xDomain;

tFin = 1;
dtPrint = 0.05;

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

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

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@%

Exp[I s]
ExpS = Series[ %,{s,0,3}]//Normal

ParametricPlot[ {{Re[%],Im[%]},{Re[%%],Im[%%]}},{s,-Pi,Pi},PlotRange->({#,#}&)@{-1.1,1.1},PlotTheme->"Detailed",ImageSize->Large]

Simplify[Re[ExpS]^2+Im[ExpS]^2,Assumptions->s∈ Reals]
Solve[%==1,s]

FindRoot[ Abs[ExpS]==1,{s,2}]

Problém: 
Ukažte, že uvažovaná RK3 metoda dá pro U  jež je vlastní hodnotou operátoru pravé strany přesně 
U(t+dt)  = (1+i s-((s^2)/2)-((i (s^3))/6)) U(t)

{0,1,2,3,9,10}

With[ {n=10},Table[Cos[k x]/.k->k1/n Pi/Δx/.Δx->1/n, {k1,%}]]//N;
Plot[ %,{x,0,1},PlotTheme->"Scientific"]
ListLinePlot[ Transpose@Table[%%,{x,0,1,1/10}],PlotTheme->"Scientific"]

Eigenvalues[ {{0,1},{2(Cos[k Δx]-1)/Δx^2,0}} ]
FullSimplify[%,Assumptions->Δx>0]

With[ {n=10},
(Table[ %/.k->k1/n Pi/Δx/.Δx->1/n, {k1,0,10}]//N)
]

Problém: 
1. Nalezněte realizaci (∂_xx) čtvrtého řádu 
2. Ve všech bodech x, kde to jde ji použijte ve funkci rhs[...]
3. Ukažte v čem jsou výsledky lepší a v čem horší 
( využijte, že  pro malé časy znáte d’Alembetovo řešení vlnové rovnice )