Energy preserving model order reduction of the nonlinear Schr\"odinger equation

TL;DR

This paper develops an energy-preserving reduced order model for the 2D nonlinear Schrödinger equation, ensuring long-term stability and efficient computation using POD-Galerkin, DEIM, and DMD methods.

Contribution

It introduces a novel energy-preserving ROM for the 2D NLSE that combines POD-Galerkin with DEIM and DMD for efficient nonlinear computations.

Findings

The ROM preserves energy and mass, ensuring stability.

POD-DMD significantly speeds up computations compared to POD-DEIM.

Both methods accurately approximate the full order model.

Abstract

An energy preserving reduced order model is developed for two dimensional nonlinear Schr\"odinger equation (NLSE) with plane wave solutions and with an external potential. The NLSE is discretized in space by the symmetric interior penalty discontinuous Galerkin (SIPG) method. The resulting system of Hamiltonian ordinary differential equations are integrated in time by the energy preserving average vector field (AVF) method. The mass and energy preserving reduced order model (ROM) is constructed by proper orthogonal decomposition (POD) Galerkin projection. The nonlinearities are computed for the ROM efficiently by discrete empirical interpolation method (DEIM) and dynamic mode decomposition (DMD). Preservation of the semi-discrete energy and mass are shown for the full order model (FOM) and for the ROM which ensures the long term stability of the solutions. Numerical simulations illustrate the preservation of the energy and mass in the reduced order model for the two dimensional NLSE with and without the external potential. The POD-DMD makes a remarkable improvement in computational speed-up over the POD-DEIM. Both methods approximate accurately the FOM, whereas POD-DEIM is more accurate than the POD-DMD.

Tables

Table 1. Table 1: Example  5.1 : The computation time (in sec.) and speed-up factors for 10 POD and 15 DEIM/DMD modes

	CPU Time (sec.)	Speedup
FOM	8024.7	-
ROM with DEIM	763.9	10.5
ROM with DMD	8.2	974.3

Table 2. Table 2: Example  5.2 : The computation time (in sec.) and speed-up factors for 20 POD/DEIM/DMD modes

	CPU Time (sec.)	Speedup
FOM	309.7	-
ROM with DEIM	27.5	11.3
ROM with DMD	0.3	956.7

Equations

$$i\frac{\partial \psi (t,\mathbf{x})}{\partial t}$$

$$\psi (0,\mathbf{x})={\psi }_0(\mathbf{x}),\mathbf{x}\in \Omega ,$$

$$\psi (t,a,y)$$

$$\psi (t,x,c)$$

$$N(\psi (t,\mathbf{x})):=$$

$$E(\psi (t,\mathbf{x})):=$$

$$N(\psi (t,\mathbf{x}))=N({\psi }_0(\mathbf{x})),E(\psi (t,\mathbf{x}))=E({\psi }_0(\mathbf{x})),\forall t\in (0,T],$$

$$i⟨\frac{\partial \psi }{\partial t},\varphi ⟩$$

$$\displaystyle\langle\psi,\phi\big{)}$$

$$L^2(\Omega )=\{f:\Omega \mapsto \mathbb{R}|{\int }_{\Omega }[f(\mathbf{x}){]}^{2}d\mathbf{x}<\infty \},$$

$$⟨f(\mathbf{x}),g(\mathbf{x}){⟩}_{L^2(\Omega )}={\int }_{\Omega }f(\mathbf{x})g(\mathbf{x})d\mathbf{x},\Vert f(\mathbf{x})\Vert :=\Vert f(\mathbf{x}){\Vert }_{L^2}=\sqrt{⟨f(\mathbf{x}),f(\mathbf{x}){⟩}_{L^2(\Omega )}}.$$

$$W$$

$$H^1(\Omega )$$

$$⟨\frac{\partial r}{\partial t},\eta ⟩$$

$$⟨\frac{\partial s}{\partial t},\vartheta ⟩$$

$$⟨r,\eta ⟩=⟨r_0,\eta ⟩,⟨s,\vartheta ⟩=⟨s_0,\vartheta ⟩,$$

$$\mbox{$\left[\!\left[w\right]\!\right]$}=w|_{E}\bm{n}_{K}+w^{e}|_{E}\bm{n}_{K^{e}},\quad\mbox{$\left\{\!\!\left\{w\right\}\!\!\right\}$}=\frac{1}{2}\big{(}w|_{E}+w^{e}|_{E}\big{)},$$

$$\mbox{$\left[\!\left[\nabla w\right]\!\right]$}=\nabla w|_{E}\cdot\bm{n}_{K}+\nabla w^{e}|_{E}\cdot\bm{n}_{K^{e}},\quad\mbox{$\left\{\!\!\left\{\nabla w\right\}\!\!\right\}$}=\frac{1}{2}\big{(}\nabla w|_{E}+\nabla w^{e}|_{E}\big{)}.$$

$$\mbox{$\left[\!\left[w\right]\!\right]$}_{\omega}=w|_{E_{l}}\bm{n}_{l}-w|_{E_{m}}\bm{n}_{l},\quad\mbox{$\left\{\!\!\left\{w\right\}\!\!\right\}$}_{\omega}=\frac{1}{2}\big{(}w|_{E_{l}}+w|_{E_{m}}\big{)},$$

$$W_h$$

$$⟨r_h(0,\mathbf{x}),\eta (\mathbf{x})⟩=⟨r_0(\mathbf{x}),\eta (\mathbf{x})⟩,⟨s_h(0,\mathbf{x}),\vartheta (\mathbf{x})⟩=⟨s_0(\mathbf{x}),\vartheta (\mathbf{x})⟩.$$

$$⟨\frac{\partial r_h}{\partial t},\eta ⟩$$

$$⟨\frac{\partial s_h}{\partial t},\vartheta ⟩$$

$$a_h(z_h,\varphi )=$$

$$\displaystyle+\sum\limits_{K\in\mathcal{T}_{h}}\int\limits_{K}V(\bm{x})\,z_{h}\,\phi\;d\bm{x}-\sum_{\omega\in\mathcal{E}^{p}_{h}}\int\limits_{\omega}\Big{(}\mbox{$\left\{\!\!\left\{\nabla z_{h}\right\}\!\!\right\}$}_{\omega}\cdot\mbox{$\left[\!\left[\phi\right]\!\right]$}_{\omega}+\mbox{$\left\{\!\!\left\{\nabla\phi\right\}\!\!\right\}$}_{\omega}\cdot\mbox{$\left[\!\left[z_{h}\right]\!\right]$}_{\omega}\Big{)}\;d\bm{s}$$

$$\displaystyle+\sum\limits_{E\in\mathcal{E}^{0}_{h}}\frac{\kappa\alpha}{h_{E}}\int\limits_{E}\mbox{$\left[\!\left[z_{h}\right]\!\right]$}\cdot\mbox{$\left[\!\left[\phi\right]\!\right]$}\;d\bm{s}+\sum_{\omega\in\mathcal{E}^{p}_{h}}\frac{\kappa\alpha}{h_{E}}\int\limits_{\omega}\mbox{$\left[\!\left[z_{h}\right]\!\right]$}_{\omega}\cdot\mbox{$\left[\!\left[\phi\right]\!\right]$}_{\omega}\;d\bm{s},$$

$$b_{h,z}(r_h,s_h;\varphi )=$$

$$r_h(t,\mathbf{x})=\sum _{i=1}^{n_K}\sum _{j=1}^{n_q}{r}_j^i(t){\phi }_j^i(\mathbf{x}),s_h(t,\mathbf{x})=\sum _{i=1}^{n_K}\sum _{j=1}^{n_q}{s}_j^i(t){\phi }_j^i(\mathbf{x}),$$

$$\mathbf{M}{\mathbf{z}}_t=\mathbf{J}(\mathbf{A}\mathbf{z}+\mathbf{b}(\mathbf{z})),$$

$$\mathbf{r}$$

$$\mathbf{s}$$

$$\bm{M}=\left[\begin{array}[]{cc}M&0\\ 0&M\\ \end{array}\right],\;\bm{A}=\left[\begin{array}[]{cc}A&0\\ 0&A\\ \end{array}\right],\;\bm{b}(\bm{z})=\left[\begin{array}[]{cc}b({\mathbf{r}},{\mathbf{s}};r_{h})\\ b({\mathbf{r}},{\mathbf{s}};s_{h})\end{array}\right],\;\bm{J}=\left[\begin{array}[]{cc}0&\text{Id}\\ -\text{Id}&0\\ \end{array}\right],$$

$$\dot{\mathbf{y}}=f(\mathbf{y}),$$

$$\frac{{\mathbf{y}}_{n+1}-{\mathbf{y}}_n}{\tau }={\int }_0^{1}f(\xi {\mathbf{y}}_{n+1}+(1-\xi ){\mathbf{y}}_n)d\xi ,$$

$$\mathbf{M}({\mathbf{z}}_{n+1}-{\mathbf{z}}_n)=\mathbf{J}\left(\frac{\tau }{2}\mathbf{A}({\mathbf{z}}_{n+1}+{\mathbf{z}}_n)+\tau {\int }_0^1\mathbf{b}(\xi {\mathbf{z}}_{n+1}+(1-\xi ){\mathbf{z}}_n)d\xi \right).$$

$$N_h({\psi }_h(t,\mathbf{x}))=\Vert {\psi }_h(t,\mathbf{x}){\Vert }^2={\int }_{\Omega }|{\psi }_h(t,\mathbf{x}){|}^{2}d\mathbf{x}={\int }_{\Omega }(r_h^2(t,\mathbf{x})+s_h^2(t,\mathbf{x}))d\mathbf{x}$$

Full text

††journal: myjournal

1 Introduction

The nonlinear Schrödinger equation (NLSE) is used frequently for modelling wave propagation phenomena in different areas of physics, chemistry and engineering. In quantum physics and chemistry, the NLSE is used in the study of Bose-Einstein condensation (BEC), where it is also called Gross-Pitaevskii equation (GPE). The numerical simulations are extremely important for predicting the long term behavior of NLSE/GPE, because experiments to create BECs are challenging and computationally expensive.

The structure preserving symplectic and multisymplectic Bridges06 ; Chen02 ; Islas01 integrators and energy preserving integrators Celledoni12ped using finite difference discretization in space, lead to long term stability of the solutions of Hamiltonian partial differential equations (PDEs) like the NLSE. In Li15 local energy-preserving multi-symplectic algorithms in time are constructed for coupled and two dimensional NLSE using pseudospectral methods or Gauss-Legendre collocation methods for the spatial discretization.

Due to the implicit nature of these integrators, at each time step, a coupled fully nonlinear system has to be solved by a fixed point or a Newton-Raphson method with high accuracy Hairer10 . The discrete energy is conserved when the nonlinear system is solved accurately up to the machine precision. Therefore solving NLSE might be extremely time consuming, especially for two and three dimensional (2D and 3D) problems.

The goal of this work is fast and accurate reduced order modeling for the 2D NLSE that preserves energy and mass. The proper orthogonal decomposition reduced order model (POD-ROM) has been widely used as a computationally efficient surrogate model in large-scale numerical simulations of complex systems. However, when it is applied to a Hamiltonian system like the NLSE, naive application of the POD method can destroy the Hamiltonian structure in the reduced order model. The stability of reduced models over long-time integration and the structure preserving properties has been recently investigated in the context of Lagrangian systems Lall03 ; Carlberg15 , and for Hamiltonian systems Beattie11 ; Chaturantabu16 . For linear and nonlinear Hamiltonian systems, a symplectic model reduction technique, proper symplectic decomposition (PSD) with symplectic Galerkin projection is constructed in Peng16 which captures the symplectic structure of Hamiltonian systems to ensure long term stability of the reduced model. The PSD is computed using three algorithms; complex SVD, cotangent lift and nonlinear programming. The PSD-DEIM method is extended with a greedy algorithm to parametric Hamiltonian systems in Hesthaven16 . For nonlinear Hamiltonian systems, the PSD is combined with the discrete empirical interpolation method (DEIM) to reduce the computational cost Peng16 ; Hesthaven16 . Recently the energy preserving average vector field (AVF) method is used as a time integrator to construct reduced order models for Hamiltonian systems Gong17 , for dissipative gradient systems like the Allen-Cahn equation Karasozen17 . In general, there exist no numerical integration method that preserves both the symplectic structure and the energy of a Hamiltonian systems. The snapshots are enriched using the gradients of the Hamiltonian at time instances in Gong17 . Besides the symplectic Hairer10 and multi-symplectic structure preserving geometric integrators Bridges06 , there are energy or Hamiltonian preserving integrators; the average vector field (AVF) method Quispel08 ; Celledoni12ped . The AVF is second order convergent in time, and preserves the time reversal symmetry of the Hamiltonian systems like the implicit mid-point rule. Energy preservation of various PDEs by the AVF method is studied in Celledoni12ped : Hamiltonian and dissipative PDEs, like the NLSE, Korteweg-de Vries equation, Sine-Gordon equation, Allen-Chan equation, Ginzburg Landau equation and Maxwell equation.

Currently the symplectic and energy preserving reduced order models are applied to one dimensional linear and non-linear Hamiltonian systems using finite differences in space; linear wave equation Peng16 ; Hesthaven16 ; Gong17 , NLSE Peng16 ; Hesthaven16 , Sine-Gordon equation Peng16 , Korteweg-de Vries equation Gong17 . The long term computation of the two dimensional NLSE with and without external potential is computationally challenging and more expensive than the one dimensional NLSE. In this paper we develop an energy stable reduced order model (ROM) for two dimensional NLSE, which inherits the energy preservation property of the full order model (FOM). For the space discretization we apply symmetric interior penalty discontinuous Galerkin (SIPG) method Arnold02 ; riviere08dgm , and for time discretization the energy preserving AVF method. The SIPG is more efficient by evaluations of the nonlinear terms than the continuous FEM Karasozen15 . The AVF method preserves the energy dissipation of gradient systems like the Allen-Cahn equation Karasozen17a and the skew-gradient structure like the diffusive FitzHugh-Nagumo equation Karasozen15 . Reduced order models, which preserve the energy dissipation of gradient and skew-gradient systems are developed in Karasozen15 ; Karasozen17 . The POD-Galerkin approach involves an offline-online splitting methodology. In the offline stage, the high fidelity or truth solutions are generated by numerical simulations of the discretized high dimensional FOM. The POD is then applied to compute an optimal subspace to fit the high fidelity data. A reduced system is constructed by projecting the high-dimensional system to this subspace. In the online stage, the reduced system is solved in the low-dimensional subspace. The primary challenge in producing the low dimensional models of the high dimensional discretized PDEs is the efficient evaluation of the nonlinearities (inner products) on the POD basis. Different POD Galerkin methods are developed to reduce the complexity of evaluating the nonlinear terms; gappy POD Everson95 ; Carlberg13 , missing point estimation Astrid08 ; Zimmermann16 , empirical and discrete empirical interpolation method (EIM/DEIM) Barrault04 ; chaturantabut10nmr . Recently the dynamic mode decomposition (DMD) has appeared as an alternative to produce low rank approximation of the nonlinear terms Alla16 . Like the DEIM, the DMD uses singular value decomposition (SVD) of the nonlinearities, which also uses sparse sampling through the gappy POD and DEIM/EIM. In Kutz13 a hybrid combination of the POD with DMD is applied to the one dimensional NLSE and to the Ginzburg-Landau equation. Recently in Charnyi17 , for the incompressible Navier-Stokes equation (NSE), energy, mass and angular momentum conserving Galerkin discretization methods are developed in convective, skew-symmetric, conservative, and rotational formulations. The conservation properties of ROMs such as energy balance and mass conservation for the NLSE are investigated theoretically and numerically Mohebujjaman17 .

Efficient computation of the reduced order basis functions is an important aspect in reduced order modelling. The reduced order basis functions in POD-DEIM and POD-DMD are computed in the offline stage using singular value decomposition (SVD), which can be computationally demanding for large snapshot matrices. Here we use the randomized singular value decomposition (rSVD) Halko11a ; Mahoney11 ; Martinsson16 as a fast and accurate alternative to the deterministic SVD to reduce the computational cost in the offline stage Alla16a ; Bistrian17 ; Erichson16 . We compare the DEIM and DMD with respect to the accuracy and efficiency. We show that ROMs preserve well the energy and mass of the FOMs through numerical simulations for two dimensional NLSE with and without an external potential. It turns out that POD-DMD is several orders faster than the POD-DEIM, but less accurate than the POD-DEIM.

The remainder of the paper is organized as follows. The NLSE is reviewed briefly in Section 2. Space-time discretization of NLSE by the SIPG and AVF method is given in Section 3. In Section 4, the energy preserving ROM, and also the POD-DEIM and POD-DMD are presented. The preservation of the invariants like the energy and mass of the NLSE, efficiency and accuracy of the ROMs are illustrated through the numerical simulations in Section 5. The paper ends with some conclusions.

2 Nonlinear Schrödinger equation

The general form of the NLSE equation with a cubic nonlinearity and an external potential $V$ is given by

[TABLE]

with an initial condition

[TABLE]

where $\psi(t,\bm{x})$ is the complex valued wave function, $i$ is the complex unity ($i=\sqrt{-1}$), $\Delta=\nabla^{2}$ is the Laplace operator, $t$ is the time variable, $T>0$ is a terminal time, $\alpha$ is a positive constant, and $\Omega\subset\mathbb{R}^{d}\;(d=1,2,3)$ is a bounded polygonal computational domain. In this paper we consider two dimensional ($d=2$) NLSE, and we set the spatial variable $\mathbf{x}=(x,y)$. The parameter $\beta$ is a dimensionless constant with negative value for focusing (or attractive) and positive value for defocusing (or repulsive) nonlinearity. The external potential function $V(\mathbf{x})$ is a given real-valued function and its specific form depends on different applications Pitaevskii03 ; Sulem99 . For BEC, it is usually chosen as a harmonic confining potential, i.e., quadratic polynomial Antoine13 ; Bao13a . The boundary conditions are usually of homogeneous Dirichlet or periodic type. We consider here periodic boundary conditions in the space time domain $(t,x,y)\in(0,T]\times[a,b]\times[c,d]$

[TABLE]

For the NLSE (1), we relate the mass (or density) functional $N(\psi(t,\bm{x}))$ and the energy functional $E(\psi(t,\bm{x}))$, given as usual by

[TABLE]

The NLSE (1) is a Hamiltonian PDE which is well-known for Hamiltonian PDEs Sulem99 ; Bao13a that it conserves the mass (3a) and the energy (3b) as the time progresses, in other words, the mass and energy at all times is equal to the mass and energy at the initial time $t=0$:

[TABLE]

with the initial condition $\psi(0,{\mathbf{x}})=\psi_{0}(\bm{x})$.

There are other important dynamical properties of the solution $\psi$ to the NLSE (1); time reversibility or symmetry, time transverse or gauge invariance, dispersion relation for NLSE without external potential. The Gross-Pitaevskii equation (GPE) is NLSE for the macroscopic wave functions, governs the properties of BECs at temperatures far below the critical condensation temperature. The GPE includes the external potential term.

There exists a huge literature for the numerical solution of the NLSE, especially for 1D problems. Numerical methods for 2D NLSE include the alternating direction implicit (ADI) method Xu12 , implicit-explicit multistep method gao16a , high order compact finite difference schemes Wang13 , and local discontinuous Galerkin method Xu05 , the implicit Crank-Nicolson finite difference (CNFD), semi-implicit finite difference (SIFD), time splitting spectral and pseudo spectral methods Antoinegpelab1 ; Antoinegpelab2 .

3 Full order model

In this section we introduce the space-time discretization of the NLSE (1) leading to the so-called full order model (FOM), and we discuss the preservation of invariants for FOM under symmetric interior penalty Galerkin (SIPG) space discretization together with average vector field (AVF) time integrator.

Applying Green’s formula to the NLSE (1), the continuous weak formulation of the problem (1) can be written as follows: for a.e. $t\in(0,T]$, find $\psi:=\psi(t,\bm{x})\in W$ such that for any $\phi:=\phi(\bm{x})\in W$

[TABLE]

where $\langle\cdot,\cdot\rangle:=\langle\cdot,\cdot\rangle_{L^{2}(\Omega)}$ denotes the $L^{2}$-inner product on the space $L^{2}(\Omega)$ of square integrable functionals over a domain $\Omega\subset\mathbb{R}^{2}$, and together with the corresponding norm $\|\cdot\|_{L^{2}(\Omega)}$:

[TABLE]

[TABLE]

The solution space $W$ is set as the space of functions in the Sobolev space $H^{1}(\Omega)$ satisfying the periodicity on the boundary $\partial\Omega$ of the domain $\Omega$:

[TABLE]

We decompose the complex valued wave function $\psi(t,\bm{x})$ into its real and imaginary parts as $\psi(t,\bm{x})=r(t,\bm{x})+is(t,\bm{x})$, then (5) can be re-written as a pair of real-valued system:

[TABLE]

where $\psi_{0}(\bm{x})=r_{0}(\bm{x})+is_{0}(\bm{x})$. In the following, we introduce the SIPG method Arnold02 ; riviere08dgm for the spatial discretization of NLSE. We remark that among the three common interior penalty Galerkin discretizations, nonsymmetric interior penalty Galerkin (NIPG), and incomplete interior penalty Galerkin (IIPG), only SIPG leads to a Hamiltonian system of ODEs Karasozen15nls .

3.1 Semi-discretization in space

Consider for $h>0$, $\{\mathcal{T}_{h}\}_{h}$ as a family of shape-regular simplicial triangulations of $\Omega$. Each mesh $\mathcal{T}_{h}$ consists of closed triangles such that $\overline{\Omega}=\bigcup_{K\in\mathcal{T}_{h}}\overline{K}$ holds. The diameter of an (triangular) element $K$ and the length of an edge $E$ are denoted by $h_{K}$ and $h_{E}$, respectively. We split the set of all edges $\mathcal{E}_{h}$ into the set $\mathcal{E}^{0}_{h}$ of interior edges and the set $\mathcal{E}^{p}_{h}$ of periodic boundary edge-pairs. An individual element of the set $\mathcal{E}^{p}_{h}$ is of the form $\omega=\{E_{l},E_{m}\}$ where $E_{l}\subset\partial K_{n_{l}}\cap\partial\Omega$, and $E_{m}\subset\partial K_{n_{m}}\cap\partial\Omega$ is the corresponding periodic edge-pair of $E_{l}$ with $l>m$. Let the edge $E$ be a common edge for two elements $K$ and $K^{e}$. For a piecewise continuous scalar function $w$, there are two traces of $w$ along $E$, denoted by $w|_{E}$ from inside $K$ and $w^{e}|_{E}$ from inside $K^{e}$. Then, the jump and average of $w$ across the edge $E$ are defined by:

[TABLE]

where $\bm{n}_{K}$ (resp. $\bm{n}_{K^{e}}$) denotes the unit outward normal to $\partial K$ (resp. $\partial K^{e}$). Similarly, for a piecewise continuous vector field $\nabla w$, the jump and average across the edge $E$ are given by

[TABLE]

For a boundary edge $E\in K\cap\partial\Omega$, we set \mbox{\left{!!\left{\nabla w\right}!!\right}}=\nabla w and \mbox{\left[!\left[w\right]!\right]}=w\bm{n}, where $\bm{n}$ is the unit outward normal vector to $\partial\Omega$. The treatment of the periodic boundary edges is the same as the interior edges, in other words, appropriate jump and average operators introduced in Vemaganti07 are used. Then, for each $\omega=\{E_{l},E_{m}\}\in\mathcal{E}^{p}_{h}$, we define the jump and average operators as follow

[TABLE]

where, for $l>m$, we use the fixed unit outward normal vector $\bm{n}_{l}$ to the boundary edge $E_{l}\subset\partial\Omega$.

For continuous finite elements methods (FEMs), the idea is to approximate the solutions $r$ and $s$ of (6) from a conforming, finite dimensional space $W_{h}\subset W$. On the other hand, we point out that in discontinuous Galerkin (dG) methods the space of solutions or test functions consist of piecewise discontinuous polynomials. That is, no continuity constraints are explicitly imposed on the state and test functions across the element interfaces. As a consequence, weak formulations must include jump terms across interfaces, and typically penalty terms are added to control the jump terms. We define the space of solution and test functions by

[TABLE]

where $\mathbb{P}^{q}(K)$ is the set of polynomials of degree at most $q$ in $K$. Note that the space of discrete solutions and the space of test functions are identical due to the weak treatment of boundary conditions in dG methods. Note also that the space $W_{h}$ is a non-conforming space, i.e., $W_{h}\not\subset W$.

Let $r_{h}(0,\bm{x}),s_{h}(0,\bm{x})\in W_{h}$ be the projections (orthogonal $L^{2}$-projection) of the initial conditions $r_{0}(\bm{x})$ and $s_{0}(\bm{x})$ onto $W_{h}$, i.e., $\forall\eta(\bm{x}),\vartheta(\bm{x})\in W_{h}$

[TABLE]

Then, the weak formulation of the system (6), discretized by the SIPG method, reads as: for a.e. $t\in(0,T]$, find $r_{h}:=r_{h}(t,\bm{x})$, $s_{h}:=s_{h}(t,\bm{x})\in W_{h}$ such that

[TABLE]

where the bi-linear form $a_{h}(\cdot,\cdot)$ and the non-linear forms $b_{h,z}(\cdot,\cdot;\cdot)$ for $z\in\{r,s\}$ are given for any $\phi\in W_{h}$ by

[TABLE]

where the parameter $\kappa$ in the above formulation is called penalty parameter which should be sufficiently large to ensure the stability of the SIPG discretization with a lower bound depending only on the polynomial degree $q$, see, e.g., riviere08dgm .

The SIPG discretized semi-discrete solutions of (9) are given as

[TABLE]

where $r_{j}^{i}(t),s_{j}^{i}(t)$, and $\varphi_{j}^{i}(\bm{x})$ are the unknown coefficients and the basis functions for $W_{h}$, respectively, for $j=1,2,\cdots,n_{q}$ and $i=1,2,\cdots,n_{K}$. The number $n_{K}$ denotes the number of (triangular) elements in $\mathcal{T}_{h}$, and $n_{q}$ is the local dimension on each element with the identity for the 2D problems $n_{q}=(q+1)(q+2)/2$, where $q$ is the degree of the polynomial order. Note that the degrees of freedom in dG methods are given by $N:=n_{K}\times n_{q}$, and throughout this paper we denote by $N$ the dimension of the high-fidelity model, i.e., the full order model (FOM). Inserting the expansions (11) into the system (9), we obtain the following semi-discrete Hamiltonian system

[TABLE]

where ${\mathbf{z}}:={\mathbf{z}}(t):=({\mathbf{r}}^{T},{\mathbf{s}}^{T})^{T}\in\mathbb{R}^{2N}$, and ${\mathbf{r}}:={\mathbf{r}}(t)\in\mathbb{R}^{N}$ and ${\mathbf{s}}:={\mathbf{s}}(t)\in\mathbb{R}^{N}$ are the unknown coefficient vectors for the solutions $r_{h}$ and $s_{h}$ with the ordered entries

[TABLE]

The other settings are given by

[TABLE]

where $M\in\mathbb{R}^{N\times N}$ is the usual mass matrix, $A\in\mathbb{R}^{N\times N}$ is the stiffness matrix corresponding to the bilinear form $a_{h}(\cdot,\cdot)$, and $b(\cdot,\cdot;z)\in\mathbb{R}^{N}$ are the vectors corresponding to the non-linear forms (inner products) $b_{h,z}(\cdot,\cdot;\cdot)$ for $z\in\{r,s\}$. The matrix ${\bm{J}}$ is the skew-symmetric matrix with Id is the $N$ dimensional identity matrix.

3.2 Time discretization

For the temporal discretization of the semi-discrete system (12), we use the energy preserving, implicit, second order convergent AVF method Celledoni12ped ; Quispel08 . For a system of ordinary differential equations (ODEs)

[TABLE]

the AVF method reads as

[TABLE]

where ${\bm{y}}_{n+1}$ is the unknown solution at time $t_{n+1}$, ${\bm{y}}_{n}$ is the known solution at time $t_{n}$, and $\tau=t_{n+1}-t_{n}$ is the uniform time-step size. The AVF method also preserves the energy of Hamiltonian systems and Poisson systems with non-constant skew-symmetric structure Hairer10 , like the KdV equation and bi-Hamiltonian systems Karasozen13 . For problems with polynomial Hamiltonians, the AVF method can be interpreted as the Runge–Kutta method Celledoni14 , where the integrals can be evaluated exactly, and the implementation is comparable to that of the implicit mid-point rule. Using the multi-symplectic formulation of Hamiltonian PDEs, the AVF method is applied as an energy-preserving wavelet collocation to the NLSE and Camassa-Holm equation Gong16 , and as local energy-preserving method to the NLSE and KdV equation Gong14 , to the 2D NLSE Li15 .

For time discretization, we divide the time interval $[0,T]$ into $N_{T}$ steps: $0=t_{0}<t_{1}<\cdots<t_{N_{T}}=T$ with the uniform time step size $\tau=t_{n}-t_{n-1}$, $n=1,2,\cdots,N_{T}$. We set $\bm{z}_{n}\approx\bm{z}(t_{n})$ as the approximate solution vector at the time instance $t=t_{n}$, $n=0,1,\ldots,N_{T}$, with $\bm{z}_{n}=[{\bm{r}}_{n}^{T},{\bm{s}}_{n}^{T}]^{T}$, and ${\bm{r}}_{n}\approx{\bm{r}}(t_{n})$ and ${\bm{s}}_{n}\approx{\bm{s}}(t_{n})$. For $t=0$, let $r_{h}(0,\bm{x}),s_{h}(0,\bm{x})\in W_{h}$ be the projections (orthogonal $L^{2}$-projections) of the initial conditions $r_{0}({\bm{x}}),s_{0}({\bm{x}})$ onto $W_{h}$, and let ${\bm{r}}_{0}$ and ${\bm{s}}_{0}$ be the corresponding initial coefficient vectors satisfying the expansions (11), and set $\bm{z}_{0}=[{\bm{r}}_{0}^{T},{\bm{s}}_{0}^{T}]^{T}$. Then, applying the AVF method to the semi-discrete system (12), the full discrete problem of the NLSE (1) reads as: for $n=0,1,\ldots,N_{T}-1$, find $\bm{z}_{n+1}$ satisfying

[TABLE]

3.3 Preservation of the invariants in FOM

The NLSE is time reversible or symmetric, i.e., the solutions are unchanged under the change of time as $t\rightarrow-t$ and taken conjugate in the equation (1). For the full discrete scheme, the time reversibility implies that the scheme remains unchanged under the operation: $(n,n+1)\longleftrightarrow(n+1,n)$ and $\psi^{n}\longleftrightarrow\psi^{n+1}$ which is also true for the midpoint rule (Crank-Nicolson finite difference scheme) Antoine13 and SIPG-AVF scheme. There are two other invariants. The first one is the time transverse or gauge invariance, i.e., the translation of the potential $V\rightarrow V+\epsilon$ changes the phase of the solution as $\psi(t,\bm{x})\rightarrow\psi(t,\bm{x})e^{-i\epsilon t}$. The density or modulus of the solutions, $\rho=|\psi|^{2}$, is unchanged. NLSE admits plane wave solutions, satisfying a dispersion relation. NLSE without external potential (when $V(\bm{x})=0$) generates plane wave solutions $\psi(t,\bm{x})=Ae^{i(\bm{(}k)\bm{x}\omega t)}$, where $A$ is the amplitude, $\omega$ is the time frequency, $\bm{k}$ is the spatial wave number. Like the Crank-Nicolson method, AVF method can not preserve the time translation invariant and the dispersion relation. For an overview of the preservation of invariants of NLSE at the discrete level and computational cost we refer to Antoine13 .

In the sequel, we show the conservation of discrete mass and discrete energy.

Theorem 3.1

The semi-discrete solution $\psi_{h}(t,\bm{x})$ through the FOM (12), discretized by the SIPG method, conserves the discrete mass

[TABLE]

exactly for all time.

Proof

To prove the discrete mass conservation, we utilize the weak formulation (9) of FOM. We choose the test functions $\eta=r_{h}$ in (9a) and $\vartheta=s_{h}$ in (9b). Then, summation of (9a) and (9b), and using the linearity of bilinear form $a_{h}$ in both arguments, yields

[TABLE]

by which we obtain

[TABLE]

From this equation we can see that the discrete mass $N_{h}(\psi_{h}(t,\bm{x}))$ is invariant for all time. ∎

For the discrete energy preservation, we first introduce the following mesh-dependent energy functional which can be regarded as SIPG counterpart of the continuous energy functional $E(\cdot)$ defined in (3b):

[TABLE]

where

[TABLE]

Note that for a continuous solution $\psi=r+is$, since $J_{h}(r)=J_{h}(s)=0$, we have $E_{h}(r,s)=E(\psi)$.

Theorem 3.2

The FOM (12) is a Hamiltonian ODE

[TABLE]

arising from the discrete energy functional $E_{h}(\bm{z})=E_{h}(r_{h},s_{h})$ defined in (15), where $\bm{J}$ is the skew-symmetric matrix and $\bm{z}=[\bm{r}^{T},\bm{s}^{T}]^{T}$ is the vector of unknown coefficients $\bm{r}$ of $r_{h}$ and $\bm{s}$ of $s_{h}$. Moreover, for the discrete energy $E_{h}$ in (15), the FOM (12) satisfies the invariance of the energy for all time.

Proof

Firstly, we claim that $\nabla_{\bm{z}}E_{h}(\bm{z})=\bm{A}\bm{z}+\bm{b}(\bm{z})$, for which then (16) is the FOM (12). For this, by the definition of the bilinear form $a_{h}$ in (10), note that the discrete energy functional $E_{h}(r_{h},s_{h})$ in (15) can be written as

[TABLE]

Then, Fréchet derivatives of $E_{h}(r_{h},s_{h})$ with respect to the components $r_{h}$ and $s_{h}$ in any directions $\eta,\vartheta\in W_{h}$ are:

[TABLE]

or in matrix-vector form

[TABLE]

Hence, with the definition of the matrices and vectors in (13), we obtain that

[TABLE]

Further, using the identity (18) and the FOM (12), the invariance of the energy $E_{h}(\bm{z})$ for all time follows as:

[TABLE]

since the matrix $\tilde{\bm{J}}:=\bm{M}^{-1}\bm{J}$ is again a skew-symmetric matrix due to the fact that $\bm{J}$ is a skew-symmetric matrix and $\bm{M}$ is symmetric.∎

We also note that since the system of ODEs (12) is Hamiltonian, the semi-discrete energy is the preserved by the AVF method Celledoni12ped .

4 Reduced order model

Because the computation of the FOM (12) is time consuming, in this section, we construct a small dimensional reduced order model (ROM) by utilizing the proper orthogonal decomposition (POD) method Kunisch01 . In addition to POD, the nonlinear vectors in the ROM are computed efficiently by discrete empirical interpolation method (DEIM) and dynamic mode decomposition (DMD). The low-rank approximation is computed in three steps: computation of the numerical solutions of the original high-dimensional system; dimensionality-reduction of the snapshot matrices by singular value decomposition (SVD); Galerkin projection of the dynamics on the low-rank subspace. The first two steps are known as the offline stage, and the last one is the online stage. Offline stage is usually expensive and online step should be fast to run in real time.

4.1 POD Galerkin projection

For the $2N$-dimensional FOM (12), the ROM of lower dimension $k\ll 2N$ is formed by the Galerkin projection of the system onto a $k$-dimensional reduced space

[TABLE]

resulting in the lower dimensional reduced solution $\psi_{r}(t,\bm{x})$ as:

[TABLE]

where, without lost of generality, we have assumed that the complex solution $\psi_{h}=r_{h}+is_{h}$ is the column vector of the real and imaginary solutions $r_{h}$ and $s_{h}$, respectively, i.e., $\psi_{h}(t,\bm{x}):=[r_{h}(t,\bm{x}),s_{h}(t,\bm{x})]^{T}$, by which the solution $\bm{z}(t)$ of the FOM (12) becomes the dG coefficient vector of $\psi_{h}(t,\bm{x})$ at time $t$. From the coefficients $z_{i}^{r}(t)$ in (19), we set the solution of the reduced system as $\bm{z}^{r}(t):=\left(z_{1}^{r}(t),\ldots,z_{k}^{r}(t)\right)^{T}$. The functions $\{u_{i}(\bm{x})\}_{i=1}^{k}$ in (19) are the orthogonal (in $L^{2}$-sense) reduced basis functions spanning the reduced space $W_{h}^{r}$. Belonging to the space $W_{h}^{r}\subset[W_{h}]^{2}$, the reduced basis functions are linear combination of the dG basis functions $\{\varphi_{j}\}$, given by

[TABLE]

Then, using the column vectors $U_{i}=\left(U_{1,i},\ldots,U_{2N,i}\right)^{T}$, which are the coefficients of the $i$-th reduced basis function, we construct the following matrix of POD modes

[TABLE]

To obtain the reduced basis functions $\{u_{i}(\bm{x})\}_{i=1}^{k}$, we need to solve the minimization problem Kunisch01

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. The above minimization problem is equivalent to the eigenvalue problem Kunisch01

[TABLE]

where $\mathcal{Z}:=[\bm{z}_{1},\ldots,\bm{z}_{N_{T}}]\in\mathbb{R}^{2N\times N_{T}}$ is the snapshot matrix whose $n$-th column vector $\bm{z}_{n}$ is the solution vector of the FOM at time $t_{n}$. Then, the matrix $\bm{U}$ of POD modes can be computed through the singular value decomposition (SVD) of the snapshot matrix $\mathcal{Z}$. In addition, between the solution vector $\bm{z}$ of FOM and the reduced solution vector $\bm{z}^{r}$ of ROM, we have the relation

[TABLE]

from where we can find the initial reduced vector $\bm{z}^{r}(0)$. For a more detail description, we refer to Karasozen15 . We finally obtain the following reduced system:

[TABLE]

with the reduced stiffness matrix and the reduced nonlinear vector

[TABLE]

Like the FOM (12), the reduced system (22) is solved in time by the AVF method.

4.2 Randomized singular value decomposition

The SVD is known to be computationally demanding for large snapshot matrices resulting from the space-time discretized high dimensional model of the PDEs. Recently randomized algorithms are used in reduced order modelling to accelerate the offline computations in DMD Alla16a ; Bistrian17 ; Erichson16 . Randomized methods for matrix computations provide an efficient computation of low-rank structures in data matrices, which are robust, reliable and computationally efficient and can be used to construct a smaller (compressed) matrix, which accurately approximates a high-dimensional data matrix. The randomized singular value decomposition (rSVD) is robust, reliable and computationally efficient, which approximates the high dimensional snapshot matrices by constructing a smaller (compressed) matrix Halko11a ; Mahoney11 ; Martinsson16 .

Given a snapshot matrix $Y\in{\mathbb{R}}^{m\times n}$ and the target rank $k\ll\min\{m,n\}$, a low rank approximation is constructed by creating a random matrix $\Omega\in{\mathbb{R}}^{n\times k}$, then the sampled matrix $X=Y\Omega\in{\mathbb{R}}^{m\times k}$. The entries of the random sampling matrix $\Omega$ are independent and identically distributed Gaussian random variables of zero mean and unit variance Alla16a , and they are created by the MatLab routine randn. Afterward the $QR$ decomposition of $X$ is computed to obtain the orthonormal matrix $Q\in{\mathbb{R}}^{m\times k}$, $X=QR$, so that

[TABLE]

is satisfied. In the last step, the SVD of the small (compressed) matrix $B=Q^{T}Y$ is computed. The approximation error of the rSVD can be decreased by introducing a small oversampling parameter $p$ (e.g. $p=2,5,10$). Instead of $k$ random vectors, $k+p$ are generated.

We use the rSVD for the SVD computations required in the POD reduced basis calculations, and the DEIM and DMD reduced basis calculations of the nonlinear terms.

4.3 Preservation of the invariants in ROM

In this section we show the preservation of discrete mass and discrete energy for the ROM (22). While the former is not difficult to show, we should state additional modifications to show the discrete energy preservation in ROM.

Theorem 4.1

The solution to the ROM (22) conserves the discrete mass

[TABLE]

exactly for all time.

Proof

The proof of the mass conservation in ROM follows very similar to the proof for the FOM in Theorem 3.1. But now, we consider the weak formulation not on the space $W_{h}$ of dimension $2N$ but on the reduced space $W_{h}^{r}\subset W_{h}$. ∎

Theorem 4.2

The ROM (22) with the discrete energy $E_{h}$ defined in (15), does not necessarily satisfy the invariance of the discrete energy for all time. But the updated ROM

[TABLE]

with the updated reduced matrix $\bm{A}^{r}=\bm{J}_{r}\bm{U}^{T}\bm{A}$ and the updated reduced vector $\bm{b}^{r}=\bm{J}_{r}\bm{U}^{T}\bm{b}(\bm{U}\bm{z}^{r})$, where the skew-symmetric matrix $\bm{J}_{r}$ is such that $\bm{U}^{T}\bm{J}=\bm{J}_{r}\bm{U}^{T}$, preserves the energy for all time.

Proof

Using the identity (18) for $\nabla_{\bm{z}}E_{h}(\bm{z})$ and the FOM (12), the reduced system (22) can be written as

[TABLE]

Then, for the discrete energy functional $E_{h}(\bm{U}\bm{z}^{r})$ of the reduced system, we have that

[TABLE]

which may be non-zero since the matrix $\bm{U}\bm{U}^{T}\bm{J}$ is not necessarily a skew-symmetric matrix. But, if we can find a skew-symmetric matrix $\bm{J}_{r}$ so that $\bm{U}^{T}\bm{J}=\bm{J}_{r}\bm{U}^{T}$, then we can obtain an equivalent equation of the updated ROM (23) as

[TABLE]

The equation (25), as a result the updated ROM (23), is now Hamiltonian since $\bm{J}_{r}$ is skew-symmetric. Then, we also get for the discrete energy functional $E_{h}$ that

[TABLE]

Hence the ROM (23) satisfies the invariance of the discrete energy for all time. Indeed, the skew-symmetric matrix $\bm{J}_{r}$ can be easily computed as

[TABLE]

where we used the $\bm{M}$-orthogonality of the POD modes $\bm{U}$, i.e., $\bm{U}^{T}\bm{M}\bm{U}=\bm{I}$. ∎

Remark 1

*The preservation of the reduced energy is valid for any reduced Hamiltonian system of ODEs (23). *

4.4 Approximation of the nonlinearities

Although the reduced model (23) is of small dimension, the computation of the nonlinear vector $\bm{b}^{r}(\bm{z}^{r})$ still depends on the dimension $2N$ of the FOM. In this section, we give two different approaches to reduce the computational complexity due to the nonlinear vector in the ROM (23): discrete empirical interpolation method (DEIM) chaturantabut10nmr and dynamic mode decomposition (DMD) Alla16a ; Alla16 .

4.4.1 Discrete empirical interpolation method (DEIM)

The DEIM aims to find an approximation $\hat{\bm{b}}(\bm{U}\bm{z}^{r})$ to the nonlinear vector $\bm{b}(\bm{U}\bm{z}^{r})$, full dimensional nonlinear part of the reduced nonlinear vector $\bm{b}^{r}(\bm{z}^{r})=\bm{J}_{r}\bm{U}^{T}\bm{b}(\bm{U}\bm{z}^{r})$, by projecting it onto a subspace of the space generated by the non-linear vectors and spanned by a basis $\{\bm{Q}_{i}\}_{i=1}^{m}$ of dimension $m\ll 2N$:

[TABLE]

where $\bm{Q}:=[\bm{Q}_{1}\;\ldots\;\bm{Q}_{m}]\in\mathbb{R}^{2N\times m}$ is the DEIM basis matrix, and $a(t)$ is the corresponding coefficient vector. Since the system (26) is overdetermined, a projection matrix $\bm{P}=[e_{\mathfrak{p}_{1}},\ldots,e_{\mathfrak{p}_{m}}]\in\mathbb{R}^{2N\times m}$ with $e_{\mathfrak{p}_{i}}$ is the $i$-th column of the identity matrix $\bm{I}\in\mathbb{R}^{2N\times 2N}$ is computed. Then, the reduced model (23) can be rewritten as:

[TABLE]

where the matrix $\bm{B}:=\bm{U}^{T}\bm{Q}(\bm{P}^{T}\bm{Q})^{-1}$ is computed once in the off-line stage (finite element discretization), and the reduced nonlinear vector $\bm{b}^{r}_{\text{deim}}(\bm{z}^{r}):=\bm{P}^{T}\bm{b}(\bm{U}\bm{z}^{r})$ requires only $m\ll 2N$ integral evaluations. Using the definition of reduced matrix $\bm{A}^{r}=\bm{J}_{r}\bm{U}^{T}\bm{A}$, we also have an equivalent formulation of (27) that

[TABLE]

where $\hat{\bm{b}}(\bm{U}\bm{z}^{r})=\bm{Q}(\bm{P}^{T}\bm{Q})^{-1}\bm{P}^{T}\bm{b}(\bm{U}\bm{z}^{r})$ is the DEIM approximation of the nonlinear part $\bm{b}(\bm{U}\bm{z}^{r})$. When DEIM approximation is not used, it requires $2N$ integral evaluations. On the other hand, the computation of the Jacobian of the nonlinear vector requires $2N\times n_{p}$ integral evaluations without DEIM, but it is only $m\times n_{p}$ with DEIM approximation.

For the details of the computation of the reduced non-linear vectors we refer to the greedy DEIM algorithm chaturantabut10nmr . For continuous finite element and finite volume discretizations, the number of flops for the computation of bilinear form and nonlinear term depends on the maximum number of neighbor cells Drohmann12 . In the case of dG discretization, due to its local nature, it depends only on the number of nodes in the local cells. For instance, in the case of SIPG with linear elements ($n_{q}=3$), it contributes only 3 nonzero integrals for each degree of freedom, hence 3 integrals have to be computed on a single triangular element Karasozen15 , whereas in the case of continuous finite elements, integral computations on $6$ neighbor cells are needed Heinkenschloss14 . Since the AVF method is an implicit time integrator, at each time step, a non-linear system of equations has to be solved by Newton’s method. The reduced Jacobian has a block diagonal structure for the SIPG discretization, which is easily invertible Karasozen15 , and requires only $O(n_{q}N)$ operations with DEIM.

Because the nonlinearity in the ROM (23) is approximated by DEIM, the discrete energy is also preserved approximately in (27), not exactly as for the ROM (23), since we have that

[TABLE]

But an upper bound for the preservation of the discrete energy by POD-DEIM can be derived as follows:

[TABLE]

where in the last row, we have added $-\bm{U}^{T}\nabla_{\bm{z}}E_{h}(\bm{U}\bm{z}^{r})$ which does not affect the identity since $[\bm{U}^{T}\nabla_{\bm{z}}E_{h}(\bm{U}\bm{z}^{r})]^{T}\bm{J}_{r}[\bm{U}^{T}\nabla_{\bm{z}}E_{h}(\bm{U}\bm{z}^{r})]=0$. Using that $\bm{U}^{T}\nabla_{\bm{z}}E_{h}(\bm{U}\bm{z}^{r})=\bm{U}^{T}\bm{A}\bm{z}^{r}+\bm{U}^{T}\bm{b}(\bm{U}\bm{z}^{r})$, we obtain

[TABLE]

[TABLE]

Finally, using the DEIM approximation error (chaturantabut10nmr, , Lemma 3.2)

[TABLE]

we obtain an upper bound for the change of the energy with respect to time:

[TABLE]

According to the bound (30), the discrete energy is preserved since $\|\frac{d}{dt}E_{h}\|\rightarrow 0$ as $\|(\bm{I}-\bm{Q}\bm{Q}^{T})\bm{b}(\bm{U}\bm{z}^{r})\|\rightarrow 0$, and that $\bm{Q}$ is orthonormal. We remark that similar error bounds for the reduced order energy are derived for the symplectic DEIM in Hesthaven16 ; Peng16 .

4.4.2 Dynamic mode decomposition (DMD)

The DMD extracts dynamically relevant spatio-temporal information content from a numerical or experimental data sets Kutz16 . It is a powerful equation-free, data-driven method to analyze complex systems. Without explicit knowledge of the dynamical system, the DMD algorithm determines eigenvalues, eigenmodes, and spatial structures for each mode. The order of orthogonal POD modes are decreasing by the energy level of POD singular values. The DMD modes are not orthogonal and each DMD mode is associated with a growth rate and each mode has a single frequency. Again the DMD modes are ordered in form of decreasing energy; unlike the POD modes, they give the energy fluctuations at different frequencies. After building the DMD basis functions of rank $m$, we approximate the full dimensional nonlinear part $\bm{b}(\bm{U}\bm{z}^{r})$ of the nonlinear term $\bm{b}^{r}=\bm{J}_{r}\bm{U}^{T}\bm{b}(\bm{U}\bm{z}^{r})$ in (23) following Alla16a ; Alla16 .

DMD is a special case of the Koopman operator Koopman31 approximating nonlinear systems via an associated infinite dimensional system. The connection between the DMD and Koopman operator was established in Mezic13 ; Rowley12 ; Schmid10dmd . The Koopman operator ${\mathcal{K}}$ acts on a set of scalar observable functions $g:{\mathcal{M}}\rightarrow{\mathbb{C}}$

[TABLE]

for the nonlinear dynamical system

[TABLE]

where ${\bm{y}}\in{\mathcal{M}}$, an n-dimensional manifold. The DMD determines the Koopman eigenvalues and modes directly from the data, when the observable is considered as state space, $g({\bm{y}})={\bm{y}}$. For construction of the DMD modes, we follow (Chapter 1,Kutz16 ). The nonlinear system (31) is approximated locally by the following linear system

[TABLE]

by constructing the discrete dynamical system

[TABLE]

for $J+1$ trajectories, and $\bm{A}=exp(\mathcal{A}\Delta t)$. Low rank eigen-decomposition of the matrix $\bm{A}$ is constructed by minimizing

[TABLE]

in the least squares sense for all trajectories $\bm{y}_{i},\;i=0,1,\dots,J$. For minimization of the approximation error (33), we consider the snapshot matrices ${\bm{G}}$ and ${\bm{G}}^{\prime}$ formed as:

[TABLE]

Then, we find the unknown matrix ${\bm{A}}_{G}$ satisfying ${\bm{G}}^{\prime}={\bm{A}}_{G}{\bm{G}}$, which is the solution of the minimization problem in the Frobenius norm

[TABLE]

so that ${\bm{A}}_{G}={\bm{G}}^{\prime}{\bm{G}}^{{\dagger}}$, where ${\dagger}$ denotes the Moore-Penrose pseudo inverse. The DMD modes are computed by the exact DMD algorithm Tu14 .

In our case, the exact DMD algorithm, Algorithm 2, is applied to the snapshot matrices ${\bm{G}}$ and ${\bm{G}}^{\prime}$ formed by the snapshots of the nonlinear terms ${\bm{b}}(\bm{z})$ in the FOM (12) at $N_{T}+1$ equally spaced time instances as:

[TABLE]

Then, we obtain the time dependent DMD approximation to the nonlinear part ${\bm{b}}(\bm{U}\bm{z}^{r}(t))$ in the ROM (23) as

[TABLE]

where $U^{\text{DMD}}=[U_{1},\ldots,U_{m}]$ are DMD basis functions of rank $m$ of the nonlinear vector ${\bm{b}}(\bm{z})$, $\alpha=[\alpha_{1},\ldots,\alpha_{m}]$ is the initial vector $\alpha=(U^{\text{DMD}})^{{\dagger}}{\bm{b}}(\bm{z}_{0})$ and $\omega_{j}=\log{(\lambda_{j})}/\Delta t$, $j=1,\ldots m$. After plugging this term into (23), we obtain the following linear ROM:

[TABLE]

The POD-DMD ROM (35) corresponds to the semi-discretized linear Schrödinger equation with an additional time dependent term coming from the linearization by POD-DMD approximation. The reduced energy of the POD-DMD contains only quadratic terms in (3b) like linear Schrödinger equation Faou09 . The AVF method is equivalent to the midpoint rule for quadratic potentials which is also energy preserving Celledoni12ped ; Cohen11 . Because the POD-DMD approximation of the cubic nonlinearity of the NLSE is reduced to a time dependent function appearing as additional term in the linear Schrödinger equation, the mass and energy are preserved with a lower accuracy than by POD-DEIM, which is shown by the numerical results in the next section. Although the dimension of the system (35) is the same as (27), the main advantage of the system (35) is that it is linear and we do not need to use the Newton’s method. Therefore the POD-DMD is significantly much faster than POD and POD-DEIM.

5 Numerical results

In this section we present numerical results for the NLSE (1) on a 2D localized rectangle $\Omega=[a,b]^{2}$ with periodic boundary conditions. In all simulations we use linear dG basis functions on a uniform $32\times 32$ rectangular grid with $2048$ triangular elements. Because the NLSE generate wave type solutions, it is not possible to capture the dynamics with few modes. For numerical simulations, the number of POD modes are fixed according to the relative information content or energy criterion $\varepsilon_{k}$ given by

[TABLE]

where $\sigma_{i}$ is the $i$-th singular value in the SVD of the snapshot matrix $\mathcal{Z}$ for computing the matrix $\bm{U}$ of POD modes, and $d_{z}$ is the rank of the snapshot matrix $\mathcal{Z}$. We set the number of POD modes as $\min_{k}\varepsilon_{k}>0.9999$ which sufficiently reflects the system characteristics. In the POD and DEIM/DMD modes computation procedure, we use rSVD algorithm, Algorithm 1, for the computation of SVD with the random matrix $\Omega$ and the oversampling parameter $p=2$ . In all numerical experiments we have used for the FOM and ROMs the same time-step sizes. All simulations are performed on a Windows 10 machine with Intel Core i7, 2.5 GHz and 8 GB using MATLAB R2014.

5.1 Defocusing NLSE with progressive wave solutions

We first consider defocusing NLSE (1) ($\beta=2$) with the progressive plane wave solution Wang13 )

[TABLE]

where $\omega=c_{1}^{2}+c_{2}^{2}-\beta|A|^{2}$ and $\alpha=2$. The initial data is evaluated by taking $t=0$ using the exact solution. Numerical solutions are obtained with the linear dG elements in the spatial domain $\Omega=[0,2\pi]^{2}$ and for the final time $T=5$, with the spatial and temporal mesh sizes $h=\pi/16$ and $\tau=0.001$, respectively. The parameter values are taken as $A=1$, $c_{1}=1$ and $c_{2}=1$.

In Fig. 1 the decay of singular values are shown for the snapshot matrices of the solutions and nonlinear terms. At the very beginning a steep decay of singular values is observed. With increasing number of the modes, the singular values decay rather slowly. We choose in the ROM the number of POD modes $k=10$ which satisfies the energy criterion (36) with $\varepsilon_{10}>0.9999$. In Fig. 2 the FOM solution and the errors between FOM and ROMs are shown at the final time. The number of DEIM and DMD modes are chosen to obtain the same level of accuracy for the error between the FOM and ROM solutions Alla16 . Both POD-DEIM and POD-DMD approximations with 15 DEIM/DMD modes for 10-rank truncation (POD modes) have almost the same accuracy.

In Fig. 3, the energy and mass errors are shown for the FOM and the ROM in $L^{\infty}$-norm. The ROM errors are larger than the FOM errors as expected, but both mass and energy are well preserved.

On the other hand, for fixed number of POD modes $k=10$, the relative errors of the solution, of the discrete energy and of the discrete mass in Fig. 4 decay monotonically with increasing number of the DEIM/DMD modes, and after they stagnate. Similar results were obtained for 1D NLSE in Alla16 .

It can be clearly seen from Table 1 that the POD-DMD is much faster than the POD-DEIM, whereas there is no significant difference in the relative errors in Fig. 4. The advantage of the DMD lies in the fact that the ROM becomes linear, whereas in the case of DEIM, the ROM is still semi-linear and the nonlinearity have to be evaluated.

In Fig. 5, we give the discrete energy and the discrete mass errors for the ROMs in $L^{\infty}$-norm, using a larger time-step size $\tau=0.01$ than for the FOM. The discrete energy and the discrete mass are still preserved but with a lower accuracy.

We show in Figure 6 that the discrete energy and the discrete mass are not conserved when a non-conservative time integrator like the backward Euler method is used.

5.2 NLSE with an external potential

We consider for $\beta=1$ the 2D defocusing NLSE (1) with the harmonic external trap potential Galati13

[TABLE]

The initial data is taken as the Gaussian

[TABLE]

in the space-time domain $\Omega=[-8,8]^{2}$ and $T=3$, with the uniform spatial and temporal mesh sizes $h=0.5$ and $\tau=0.01$, respectively. The diffusion constant is set as $\alpha=0.5$. The decay of the singular values in Fig. 7 is similar to the previous Example, i.e. NLSE without the external potential.

The POD-DEIM and POD-DMD solutions in Fig. 8 have again the same accuracy, but requiring more POD, DEIM and DMD modes than the NLSE without the external potential. We set the number of POD modes $k=20$ according to the energy criterion $\varepsilon_{20}>0.9999$, further increase would not contribute a valuable accuracy.

Again the discrete energy and the discrete mass are well preserved in Fig. 9 by the FOM and the ROMs.

The solution, discrete energy and discrete mass errors in Fig. 10 reach a plateau by increasing number of DEIM/DMD modes, and the POD-DEIM is much more accurate than the POD-DMD. The computational efficiency of the ROM with the DMD is clearly seen in Table 2. Similar results are obtained for the one dimensional Burger’s equation and NLSE in Alla16 .

6 Conclusions

We have performed a comparative study using the energy preserving ROM applied to 2D NLSE, showing the numerical efficiency and accuracy of the reduced order approximations and the complexity reduction of the nonlinear terms with DEIM and DMD. The discrete energy and the discrete mass preservation of POD-DEIM and POD-DMD reduced order models have been shown to be accurate and efficient for capturing the spatio-temporal dynamics of the 2D NLSE with substantial reduction in both dimension and computational time. This is clearly demonstrated in numerical simulations by the comparative computation times (speed-ups) and relative errors of the reduced order systems with respect to FOMs. The POD-DMD is always faster and the POD-DEIM is in general more accurate.

Acknowledgments: The authors would like to thank the reviewers for the comments and suggestions that helped to improve the manuscript.