CITE THIS NOTEBOOK: Koopman analysis of the periodic Korteweg–de Vries equation by Jeremy Parker. Wolfram Community April 6 2023.
ORIGINAL ARTICLE: Jeremy Parker, Claire Valva, Koopman analysis of the periodic Korteweg–de Vries equation, Chaos 33, 043102 (2023); https://doi.org/10.1063/5.0137088.
ORIGINAL ARTICLE: Jeremy Parker, Claire Valva, Koopman analysis of the periodic Korteweg–de Vries equation, Chaos 33, 043102 (2023); https://doi.org/10.1063/5.0137088.
Article Abstract: The eigenspectrum of the Koopman operator enables the decomposition of nonlinear dynamics into a sum of nonlinear functions of the state space with purely exponential and sinusoidal time dependence. For a limited number of dynamical systems, it is possible to find these Koopman eigenfunctions exactly and analytically. Here, this is done for the Korteweg–de Vries equation on a periodic interval using the periodic inverse scattering transform and some concepts of algebraic geometry. To the authors’ knowledge, this is the first complete Koopman analysis of a partial differential equation, which does not have a trivial global attractor. The results are shown to match the frequencies computed by the data-driven method of dynamic mode decomposition (DMD). We demonstrate that in general, DMD gives a large number of eigenvalues near the imaginary axis and show how these should be interpreted in this setting. Dynamic mode decomposition (DMD) is a widely used computational method for analyzing spatiotemporal data from experiments, observations, and numerics. The connection to the mathematical Koopman operator means that we can understand the behavior of DMD by analytically applying the Koopman operator to integrable partial differential equations. One non-trivial example is the Korteweg–de Vries equation on a periodic domain, which admits both wavelike and soliton solutions, and can be solved analytically via the inverse scattering method.
The Koopman operator, also known as the composition operator, is a linear map on the space of observables for a given dynamical system. Koopman analysis, the eigendecomposition of this operator, enables the interpretation of complicated nonlinear dynamics in a linear framework. Koopman analysis is known for its close relationship with the machine learning method of dynamic mode decomposition (DMD). Through applying Koopman analysis to simple model systems, we can learn about the behaviour of DMD on real experimental and observational data of physical systems. In this work we apply Koopman analysis to the periodic Korteweg-de Vries (KdV) equation, exploiting the special properties of this famous integrable PDE, and this allows us to interpret the DMD results for weakly nonlinear wave equations. The notebook contains all the steps for a simple example of this, with sinusoidal initial conditions, making a crude approximation of the Koopman operator.
We consider periodic-in-space solutions to the KdV equation. In this case, we use a simple sinusoidal initial condition. Solving numerically:
In[]:=
ufun=NDSolveValue[{4D[u[t,x],t]==6u[t,x]D[u[t,x],x]+D[u[t,x],{x,3}],u[0,x]==Sin[x],u[t,0]==u[t,2Pi]},u,{t,0,5},{x,0,2Pi}]
Out[]=
InterpolatingFunction
We see that despite the initial condition, the solution seems to break down to a periodic pattern with two solitons, which nonlinearly interact.
In[]:=
Plot[ufun[1,x],{x,0,4Pi}]
Out[]=
The inverse scattering method tells us the solutions to the KdV equation (and other similar PDEs) can be represented as ‘scattering data’. For periodic solutions, the theory is rather complicated, and we refer readers to the paper and the textbooks referenced therein. Very briefly, the solution is associated with a Riemann surface defined by branch points which arise from solving the Schrödinger equation, using the initial condition for the KdV equation as the potential for the Schrödinger equation. Given that our solution seems to have two solitons, we might imagine that the solution is well approximated by a Riemann surface of genus 2, which means we have to find 5 branch points.
In[]:=
Now, we can calculate this scattering data for the sinusoidal initial condition, which has period 2π.
In[]:=
{U,W,B,Dp,E1,E2,E3,E4,E5}=ScatteringData[Sin,2Pi]
Out[]=
{{0.-0.999958,0.-1.99985},{0.-0.035859,0.+1.9302},{{-3.17073-5.55112×,-1.92915+1.11022×},{-1.92915-1.11022×,-7.24709-8.32667×}},{1.94911×+4.71533,-2.84661×+3.15202},-0.3784892211426357789021602661665435445451639135359955812276666653,-0.3476691233941287516171238422667866952987680067594475267560452207,0.5947999368669622075808130542114134716325662965985543451097328241,0.9180581545629546803538997269207609914854115527844906760982172243,1.293166264228380726592636941754073523836024066403538065944824668}
-17
10
-16
10
-16
10
-17
10
-12
10
-12
10
Let’s plot this Riemann surface of two sheets, as well as the contours used for the integrals. The contours used act as a choice of basis for the inverse scattering method.
In[]:=
PrettyPlot[E1_,E2_,E3_,E4_,E5_]:=(sheet1[z_]:=Sqrt[Abs[(z-E1)(z-E2)(z-E3)(z-E4)(z-E5)]]Exp[I(Arg[z-E1]+Arg[z-E2]+Arg[z-E3]+Arg[z-E4]+Arg[z-E5])/2];sheet2[z_]:=-sheet1[z];ps=ContourPlot3D[z^2==Re[Sqrt[(x+Iy-E1)(x+Iy-E2)(x+Iy-E3)(x+Iy-E4)(x+Iy-E5)]]^2,{x,-1.5,3},{y,-1.5,1.5},{z,-10,10},Mesh->None,ContourStyle->Opacity[0.8]];a1[ϕ_]:=(E2+E3)/2+Exp[Iϕ](E3-E2)/1.99;a2[ϕ_]:=(E4+E5)/2+Exp[Iϕ](E5-E4)/1.9;b1[ϕ_]:=(-1+(E2+E3)/2)/2+Exp[Iϕ](1+(E2+E3)/2)/2;b2[ϕ_]:=(-1+(E4+E5)/2)/2+Exp[Iϕ](1+(E4+E5)/2)/2;p1=ParametricPlot3D[{Re[a1[ϕ]],Im[a1[ϕ]],Re[sheet1[a1[ϕ]]]},{ϕ,2Pi,0},PlotStyle->Green];p2=ParametricPlot3D[{Re[a2[ϕ]],Im[a2[ϕ]],Re[sheet1[a2[ϕ]]]},{ϕ,2Pi,0},PlotStyle->Yellow];p3=ParametricPlot3D[{Re[b1[ϕ]],Im[b1[ϕ]],Re[sheet1[b1[ϕ]]]},{ϕ,0,Pi},PlotStyle->Red];p4=ParametricPlot3D[{Re[b1[ϕ]],Im[b1[ϕ]],Re[sheet2[b1[ϕ]]]},{ϕ,-Pi,0},PlotStyle->{Red,Dashed}];p5=ParametricPlot3D[{Re[b2[ϕ]],Im[b2[ϕ]],Re[sheet1[b2[ϕ]]]},{ϕ,0,Pi},PlotStyle->Blue];p6=ParametricPlot3D[{Re[b2[ϕ]],Im[b2[ϕ]],Re[sheet2[b2[ϕ]]]},{ϕ,-Pi,0},PlotStyle->{Blue,Dashed}];Style[Show[ps,p1,p2,p3,p4,p5,p6,PlotRange->All,BoxRatios->{1,1,1}]])
In[]:=
PrettyPlot[E1,E2,E3,E4,E5]
Out[]=
We can now compute the Koopman eigenfunctions and Koopman modes, as per the paper. The Koopman modes are sums with rather complicated coefficients, which are most easily calculated using a recursive function. We use Mathematica’s FunctionCompile to make this efficient.
In[]:=
Coeff=FunctionCompile[Function[{Typed[n,"MachineInteger"],Typed[m1,"MachineInteger"],Typed[m2,"MachineInteger"],Typed[B,TypeSpecifier["PackedArray"]["ComplexReal64",2]],Typed[mlim,"MachineInteger"]},Module[{Loop},Loop=Function[{q,Typed[m1s,TypeSpecifier["PackedArray"]["MachineInteger",1]],Typed[m2s,TypeSpecifier["PackedArray"]["MachineInteger",1]],m1tot,m2tot},If[q==n,If[m1==m1tot&&m2==m2tot,Complex[0.,0.],Exp[Sum[B[[1,1]]m1s[[k]]^2+2B[[1,2]]m1s[[k]]m2s[[k]]+B[[2,2]]m2s[[k]]^2,{k,1,n-1}]/2+(B[[1,1]](m1-m1tot)^2+2B[[1,2]](m1-m1tot)(m2-m2tot)+B[[2,2]](m2-m2tot)^2)/2]],(*Strictlyspeaking,thelasttermcanexceedmlim,thoughIexpectthesecontributionstobevanishinglysmall*)Sum[If[n1==0&&n2==0,Complex[0.,0.],Loop[q+1,Append[m1s,n1],Append[m2s,n2],m1tot+n1,m2tot+n2]],{n1,-mlim,mlim},{n2,-mlim,mlim}]]];If[n==1,If[m1==0&&m2==0,Complex[0.,0.],Exp[(B[[1,1]]m1^2+2B[[1,2]]m1m2+B[[2,2]]m2^2)/2]],Loop[1,CreateTypeInstance["PackedArray"::["MachineInteger",1],CreateTypeInstance["CArray"::["MachineInteger"],0],0],CreateTypeInstance["PackedArray"::["MachineInteger",1],CreateTypeInstance["CArray"::["MachineInteger"],0],0],0,0]]]]]
Out[]=
CompiledCodeFunction
In[]:=
KoopmanMode[x_,m1_,m2_]:=-2(Dot[U,{m1,m2}])^2Exp[Dot[U,{m1,m2}]x]Sum[(-1)^n/nCoeff[n,m1,m2,B,mlim],{n,1,nlim}]
In[]:=
KoopmanEigenfunction[m1_,m2_]:=Exp[Dot[Dp,{m1,m2}]]
In[]:=
Reconstruction[x_,t_]:=Re[Sum[Exp[Dot[Wt,{m1,m2}]]KoopmanMode[x,m1,m2]KoopmanEigenfunction[m1,m2],{m1,-mlim,mlim},{m2,-mlim,mlim}]]
Finally, let’s plot the reconstruction against the numerical solution of the PDE at three different times.
In[]:=
mlim=2;nlim=2;Plot[{Reconstruction[x,0],ufun[0,x],Reconstruction[x,0.5],ufun[0.5,x],Reconstruction[x,1],ufun[1,x]},{x,0,4Pi},PlotStyle->{{Black,Dashed},{Black},{Blue,Dashed},{Blue},{Red,Dashed},{Red}}]
Out[]=
And again with more sums included in the series:
In[]:=
mlim=3;nlim=4;Plot[{Reconstruction[x,0],ufun[0,x],Reconstruction[x,0.5],ufun[0.5,x],Reconstruction[x,1],ufun[1,x]},{x,0,4Pi},PlotStyle->{{Black,Dashed},{Black},{Blue,Dashed},{Blue},{Red,Dashed},{Red}}]
Out[]=
The supplementary materials included with the paper on the Chaos journal website have additional functions, including the (much slower) genus 3 representation.