- Email: [email protected]
Model order reduction for nonlinear Schrödinger equation
Applied Mathematics and Computation 258 (2015) 509–519
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Model order reduction for nonlinear Schrödinger equation Bülent Karasözen a,⇑, Canan Akkoyunlu b, Murat Uzunca c a b c
Department of Mathematics and Institute of Applied Mathematics, Middle East Technical University, 06800 Ankara, Turkey Department of Mathematics and Computer Sciences, Istanbul Kültür University, 34156 Istanbul, Turkey Department of Mathematics, Middle East Technical University, 06800 Ankara, Turkey
a r t i c l e
i n f o
Keywords: Nonlinear Schrödinger equation Proper orthogonal decomposition Model order reduction Error analysis
a b s t r a c t We apply the proper orthogonal decomposition (POD) to the nonlinear Schrödinger (NLS) equation to derive a reduced order model. The NLS equation is discretized in space by finite differences and is solved in time by structure preserving symplectic mid-point rule. A priori error estimates are derived for the POD reduced dynamical system. Numerical results for one and two dimensional NLS equations, coupled NLS equation with soliton solutions show that the low-dimensional approximations obtained by POD reproduce very well the characteristic dynamics of the system, such as preservation of energy and the solutions. Ó 2015 Elsevier Inc. All rights reserved.
1. Introduction The nonlinear Schrödinger (NLS) equation arises as the model equation with second order dispersion and cubic nonlinearity describing the dynamics of slowly varying wave packets in nonlinear fiber optics, in water waves and in Bose–Einstein condensate theory. We consider the NLS equation
iWt þ Wxx þ cjWj2 W ¼ 0;
ð1Þ
with the periodic boundary conditions Wðx þ L; tÞ ¼ Wðx; tÞ. Here W ¼ Wðx; tÞ is a complex valued function, c is a parameter pffiffiffiffiffiffiffi and i ¼ 1. The NLS equation is called ‘‘focusing’’ if c > 0 and ‘‘defocusing’’ if c < 0; for c ¼ 0, it reduces to the linear Schrödinger equation. In last two decades, various numerical methods were applied for solving NLS equation, among them are the well-known symplectic and multisymplectic integrators and discontinuous Galerkin methods. There is a strong need for model order reduction techniques to reduce the computational costs and storage requirements in large scale simulations, yielding low-dimensional approximations for the full high-dimensional dynamical system, which reproduce the characteristic dynamics of the system. Among the model order reduction techniques, the proper orthogonal decomposition (POD) is one of the most widely used method. It was first introduced for analyzing coherent structures and turbulent flow in numerical simulation of fluid dynamics equations [5]. It has been successfully used in different fields including signal processing, fluid dynamics, parameter estimation, control theory and optimal control of partial differential equations. In this paper, we apply the POD to the NLS equation. To the best of our knowledge, there is only one paper where POD is applied to NLS equation [4], where only one and two modes approximations of the NLS equation are used in the Fourier domain in connection with mode-locking ultra short laser applications. In this paper, the NLS equation being a semi-linear partial differential equation (PDE) is discretized in space and time by preserving the symplectic structure and
⇑ Corresponding author. E-mail addresses: [email protected] (B. Karasözen), [email protected] (C. Akkoyunlu), [email protected] (M. Uzunca). http://dx.doi.org/10.1016/j.amc.2015.02.001 0096-3003/Ó 2015 Elsevier Inc. All rights reserved.
510
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
the energy (Hamiltonian). Then, from the snapshots of the fully discretized dynamical system, the POD basis functions are computed using the singular value decomposition (SVD). The reduced model consists of Hamiltonian ordinary differential equations (ODEs), which indicates that the geometric structure of the original system is preserved for the reduced model. The semi-disretized NLS equations and the reduced equations are solved in time using Strang splitting and mid-point rule. A priori error estimates are derived for POD reduced model, which is solved by mid-point rule. It turns out that most of the energy of the system can be accurately approximated by using few POD modes. Numerical results for a NLS equation with soliton solutions confirm that the energy of the system is well preserved by POD approximation and the solution of the reduced model are close to the solution of the fully discretized system. The paper is organized as follows. In Section 2, the POD method and its application to semi-linear dynamical systems are reviewed. In Section 3, a priori error estimators are derived for the mid-point time-discretization of semi-linear PDEs. Numerical solution of the semi-discrete NLS equation and the POD reduced form are described in Section 4. In the last section, Section 5, the numerical results for the reduced order models of NLS equations are presented. 2. The POD approximation for semi-linear PDEs In the following, we briefly describe the important features of the POD reduced order modeling (ROM); more details can be found in [11]. In the first step of the POD based model order reduction, the set of snapshots, the discrete solutions of the nonlinear PDE, are collected. The snapshots are usually equally spaced in time corresponding to the solution of PDE obtained by finite difference or finite element method. The snapshots are then used to determine the POD bases which are much smaller than the snapshot set. In the last step, the POD reduced order model is constructed to obtain approximate solutions of the PDE. We mention that the choice of the snapshots representing the dynamics of the underlying PDE is crucial for the effectiveness of POD based reduced model. Let X be a real Hilbert space endowed with inner product h; iX and norm kkX . For y1 ; . . . ; yn 2 X, we set n V ¼ spanfy1 ; ; yn g, as the ensemble consisting of the snapshots yj j¼1 . In the finite difference context, the snapshots can be viewed as discrete solutions yj 2 Rm at time instances tj ; j ¼ 1; . . . ; n, and ½y1 ; . . . ; yn 2 Rmn denotes the snapshot matrix. Let fwk gdk¼1 denote an orthonormal basis of V of dimension d. Then, any yj 2 V can be expressed as
yj ¼
d X hyj ; wk iX wk ;
j ¼ 1; . . . ; n:
ð2Þ
k¼1
The POD is constructed by choosing the orthonormal basis such that for every l 2 f1; . . . ; dg, the mean square error between the elements yj ; 1 6 j 6 n, and the corresponding l th partial sum of (2) is minimized on average:
2 n l X X ~ ~ a y hy ; i u u j k X k ; j j ~ 1 ;...;u ~ l 2X u j¼1 k¼1
~i ; u ~ j iX ¼ dij ; hu
min
1 6 i; j 6 l;
ð3Þ
X
where aj ’s are non-negative weights. Throughout this paper, we take the space X ¼ Rm endowed with the weighted inner product hu; v iW ¼ uT W v with the diagonal elements of the diagonal matrix W, and also aj ’s are the trapezoidal weights so that we obtain all the computations in L2 -sense. Under these choices, the solution of the above minimization problem is given by the following theorem: Theorem 1 [11, Theorem 1.3.2]. Let Y ¼ ½y1 ; . . . ; yn 2 Rmn be a given snapshot matrix with rank d 6 minfm; ng. For a b ¼ W 1=2 YD1=2 2 Rmn . Further, symmetric, positive definite matrix W 2 Rmm and a diagonal matrix D ¼ diagfa1 ; . . . ; an g, set Y b¼U b RV b T be the SVD of Y b , where U b ¼ ½u b ¼ ½v ^1 ; . . . ; u ^ m 2 Rmm ; V ^1; . . . ; v ^ n 2 Rnn are orthogonal matrices and the matrix let Y mn b . Then, for any R2R is all zero but first d diagonal elements are the nonzero singular values, r1 P r2 P . . . P rd , of Y l 2 f1; . . . ; dg, the solution to
2 n l X X ~ k iW u ~k min m aj yj hyj ; u ; ~ 1 ;...;u ~ l 2R u j¼1 k¼1
~i ; u ~ j iW ¼ dij ; hu
1 6 i; j 6 l;
ð4Þ
W
l
^ i gi¼1 . Moreover, we have r2i ¼ ki is given by the orthonormal, with respect to the weighted inner product, vectors fwi gli¼1 ¼ fW 1=2 u bT Y b. where ki ’s are the eigenvalues of the correlation matrix Y We consider the following initial value problem for POD approximation
_ yðtÞ ¼ AyðtÞ þ f ðt; yðtÞÞ; m
m
t 2 ½0; T;
yð0Þ ¼ y0 ;
ð5Þ
where f : ½0; T R ! R is continuous in both arguments and locally Lipschitz-continuous with respect to the second argument. The semi-discrete form of NLS Eq. (1) is a semi-linear equation as (5) where the cubic nonlinear part is locally Lipschitz
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
511
l continuous. Suppose that we have determined a POD basis wj j¼1 of rank l 2 f1; . . . ; dg in Rm , then we make the ansatz
yl ðtÞ ¼
l X hyl ðtÞ; wj iW wj ; |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} j¼1
t 2 ½0; T:
ð6Þ
¼:ylj ðtÞ
Substituting (6) in (5), we obtain the reduced model l l X X ylj ðtÞAwj þ f ðt; yl ðtÞÞ; y_ lj ðtÞwj ¼ j¼1
t 2 ½0; T;
l X ylj ð0Þwj ¼ y0 :
j¼1
ð7Þ
j¼1
The POD approximation (7) holds after projection on the l dimensional subspace V l ¼ spanfw1 ; . . . ; wl g. From (7) and hwj ; wi iW ¼ dij , we get
y_ li ðtÞ ¼
l X ylj ðtÞhAwj ; wi iW þ hf ðt; yl ðtÞÞ; wi iW
ð8Þ
j¼1
for 1 6 i 6 l and t 2 ð0; T. Let us introduce the matrix
B ¼ fbij g 2 Rll ;
bij ¼ hAwj ; wi iW ;
the non-linearity F ¼ ðF 1 ; ; F l ÞT : ½0; T Rl ! Rl by
* F i ðt; yÞ ¼
f t;
l X
! yj wj ; wi
j¼1 l
and the vector y ¼
+
T ðyl1 ; . . . ; yll Þ
y_ l ðtÞ ¼ Byl ðtÞ þ Fðt; yl ðtÞÞ;
;
t 2 ½0; T;
y ¼ ðy1 ; ; yl Þ 2 Rl ;
W
: ½0; T ! Rl . Then, (8) can be expressed as
t 2 ð0; T:
ð9Þ
The initial condition of the reduced system is given by yl ð0Þ ¼ y0 with
y0 ¼ ðhy0 ; w1 iW ; . . . ; hy0 ; wl iW ÞT 2 Rl : The system (9) is called the POD-Galerkin projection for (5). The ROM is constructed with POD basis vectors fwi gli¼1 of rank l. In case of l d, the l-dimensional reduced system (9) is a low-dimensional approximation for (5). The POD basis can also be computed using eigenvalues and eigenvectors. We prefer singular value decomposition, because it is more accurate than the computation of the eigenvalues. The singular values decay up to machine precision, where the eigenvalues stagnate several orders above due the fact ki ¼ r2i [3]. The choice of l is based on heuristic considb which is erations combined with observing the ratio of the modeled energy to the total energy contained in the system Y expressed by the relative information content (RIC)
EðlÞ ¼
l X
r2i 100;
i¼1
b are normalized, i.e., Pd r2 ¼ 1 [8]. The total energy of the when the squared singular values of the snapshot matrix Y i¼1 i system is contained in a small number of POD modes. In practice, l is chosen by guaranteeing that EðlÞ capturing at least % 99 of total energy of the system. 3. POD error analysis for the mid-point rule In this section we derive a priori error estimates for the mid-point rule. A priori error estimates for POD method are obtained for linear and semi-linear parabolic equations using backward Euler and Crank–Nicholson (trapezoidal rule) time discretization in [7]. When implicit midpoint rule is applied for solving the reduced model (9), we obtain the discrete system for the sequence fY j gnj¼1 in V ln ¼ spanfwn1 ; . . . ; wnl g (l 6 d)
Y j Y j1 n 1 Y j þ Y j1 ; wni ; wi ¼ AðY j þ Y j1 Þ þ f t; ; 2 Dt 2 W W hY 1 ; wni iW ¼ hy0 ; wni iW ;
ð10Þ
i ¼ 1; . . . ; l;
where Y j denotes an approximation for yl at the time t j . We are interested in estimating use the following decomposition for the error
ð11Þ 2 m j¼1 aj yðt j Þ Y j W . For u 2 R . We
Pn
512
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
yðtj Þ Y j ¼ yðtj Þ Pln yðt j Þ þ Pln yðtj Þ Y j ¼ .lj þ #lj ; where
.lj ¼ yðtj Þ Pln yðtj Þ; #lj ¼ Pln yðtj Þ Y j , and Pln : Rm ! V ln is the projection given by Pln u ¼
l P n
l X hu; wni iW wni ;
W
i¼1
¼ 1:
l
Using that fwni gi¼1 is the POD basis of rank l, we have the estimate for the terms involving
.lj
2 n l n 2 n 2 d X X X X X n n aj yðtj Þ < yðtj Þ; wi > wi ¼ aj yðtj Þ P ln yðtj Þ ¼ aj .lj ¼ r2i : W W j¼1 i¼1 j¼1 j¼1 i¼lþ1
ð12Þ
W
l ¼ ð#l #l Þ=Dt, we obtain Next, we estimate the terms involving #lj . Using the notation @# j j j1
yðtj Þ yðtj1 Þ Y j Y j1 n Pln ; wi Dt Dt W
Dt Y j þ Y j1 Dt Y j þ Y j1 þ f y tj f ; wni ¼ A y tj þ hwlj þ zlj ; wni iW ; 2 2 2 2 W
l ; wn i ¼ h@# i W j
ð13Þ
where
zlj ¼ Pln
yðt j Þ yðt j1 Þ yðt j Þ yðt j1 Þ ; Dt Dt
wlj ¼
yðt j Þ yðt j1 Þ Dt : y_ tj 2 Dt
Choosing wni ¼ #lj þ #lj1 , and using Lipschitz-continuity of f and the Cauchy–Schwarz inequality in (13) and noting that
2
2 l ; #l þ #l i ¼ 1 h@# #lj #lj1 : j j j1 W W W Dt we arrive at
2 l #j
W
2 Y þY 6 #lj1 þ DtðkAkW þ Lf Þyðt j D2tÞ j 2 j1 #lj þ #lj1 W W W l l l l þ Dt wj þ zj #j þ #j1 : W
W
ð14Þ
W
Using the Taylor series approximation yðt j D2tÞ ðyðt j Þ þ yðt j1 ÞÞ=2, applying Cauchy–Schwarz and Young’s inequalities, and collecting the terms yields
2
2 þ kzlj k2W þ kwlj k2W ; ð1 c1 DtÞk#lj k2W 6 ð1 þ c1 DtÞk#lj1 k2W þ Dt c2 .lj þ .lj1 W
W
ð15Þ
where the positive constants c1 and c2 depend on kAkW , Lipschitz constant Lf and a constant C Tay related to the Taylor series approximation. For 0 < Dt 6 2c11 , using the fact that Dtj 6 T, we have
1 6 1 þ 2c1 Dt 1 c 1 Dt
and ð1 þ 2c1 DtÞ j 6 e2c1 T ; ð1 þ c1 DtÞ j 6 ec1 T :
ð16Þ
Summation on j in (15) with using (16) and Cauchy–Schwarz inequality yields, j 2 X l 2 l . þ .l 2 þ kzl k2 þ kwl k2 ; #j 6 C Dt2 k W k W k W k1 W W
ð17Þ
k¼1
where the positive constant C depend on c1 ; c2 ; T and j. Next, we estimate the term involving wlk :
Dt 2
2 Z T j j X X l 2 e 2 Dt 4 w ¼ Dt2 yðt k Þ yðt k1 Þ y_ t k Dt 6 C kyttt ðtÞk2 dt k W 2 D t 0 W k¼1 k¼1
e depending on y, but independent of n. for a constant C
ð18Þ
513
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
_ k Þ: Now, we estimate the term involving zlk . Adding and subtracting Pln yðt
2
l 2 yðt k Þ yðt k1 Þ l z ¼ Pl yðtk Þ yðtk1 Þ Pl yðt _ _ Þ þ P Þ yðt k k n n n k W Dt Dt W 2 yðt Þ yðt Þ k k1 2 þ 4kPl yðt _ kÞ _ k Þ yðt _ k Þk2W 6 ð2kPln kW þ 4Þ n yðt Dt W 2
D t yðt Þ yðt Þ k k1 þ 4kPl yðt _ k Þ yðt _ k Þk2W þ C Tay Dt2 6 6 n y_ tk 2 Dt W _ k Þ yðt _ k Þk2W þ C Tay Dt 2 ; ¼ 6kwlk k2W þ 4kPln yðt
ð19Þ
€ðnÞ for some n 2 ðtk D2t ; t k Þ. Note that ak ’s are the trapezoidal weights which satisfy Dt 6 ak , and combining of where C Tay ¼ y (18) and (19) leads
Dt 2
j n X X b Dt 4 ; _ k Þ yðt _ k Þk2W þ C kzlk k2W þ kwlk k2W 6 8 ak kPln yðt k¼1
ð20Þ
k¼1
b depends on C e ; ky k 2 where C ttt L ð0;T;Rm Þ and C Tay . Imposing the estimate (20) in (17), we obtain
2 n n X X l 2 b Dt 4 : _ k Þ yðt _ k Þk2W þ C C 8ak kPln yðt #j 6 4C ak k.lk kW þ C W
k¼1
ð21Þ
k¼1
In addition, we have that
Pn
k¼1
estimate to the term involving
#lj
n X
d X
j¼1
i¼lþ1
aj k#lj k2W 6 C Dt4 þ
Pd
_ j Þ yðt _ j Þk2W ¼ ak 6 T and kPln yðt
2 n _ i¼lþ1 jhyðt j Þ; wi iW j .
Using these identities, we arrive at the
as
r2i þ
!!
n X
_ j Þ; wni iW j2 aj jhyðt
ð22Þ
;
j¼1
where C is dependent on y; T, but independent of n and l. Finally, combining the estimates (12) and (22), we obtain the error estimate n X
n X
n X
n X
d X
j¼1
j¼1
j¼1
j¼1
i¼lþ1
aj kyðtj Þ Y j k2W ¼
aj k#lj þ .lj k2W 6 2
aj k#lj k2W þ 2
aj k.lj k2W 6 C E
2r2i þ
n X
!
!
_ j Þ; wni iW j2 þ Dt 4 ; aj jhyðt
j¼1
where C E is dependent on y; T, but independent of n and l. As for the backward Euler and Crank–Nicolson method [7], the error between the reduced and the unreduced solutions depend for the mid-point rule on the time discretization and on the number of not modelled POD snapshots. 4. Discretization of NLS equation One dimensional NLS Eq. (1) can be written by decomposing W ¼ p þ iq in real and imaginary components
pt ¼ qxx cðp2 þ q2 Þq;
qt ¼ pxx þ cðp2 þ q2 Þp;
ð23Þ
as an infinite dimensional Hamiltonian PDE in the phase space u ¼ ðp; qÞ
u_ ¼ D
dH ; du
H¼
Z
1 2 c 2 px þ q2x ðp2 þ q2 Þ dx; 2 2
D¼
0
1
1 0
T
:
After discretizing the Hamiltonian in space by finite differences,
H ¼ Dx
!
n X pjþ1 pj 2 qjþ1 qj 2 c 2 1 2 þ ðpj þ q2j Þ ; 2 Dx Dx 2 j¼1
ð24Þ
we obtain the semi-discretized Hamiltonian ode’s
pt ¼ Aq cqðp2 þ q2 Þ;
qt ¼ Ap þ cpðp2 þ q2 Þ;
where A is the circulant matrix
0 2 B 1 B B A¼B B @ 1
1 1 1 C 2 1 C .. .. .. C C: . . . C 1 2 1 A 1 2
ð25Þ
514
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
4.1. Reduced order model for NLS equation l l Suppose that we have determined POD bases wj j¼1 and /j j¼1 of rank l ¼ f1; . . . ; dg in Rm . Then we make the ansatz
pl ¼
l X
pj ðtÞwj ðxÞ;
ql ¼
j¼1
l X qj ðtÞ/j ðxÞ;
ð26Þ
j¼1
where pj ¼ hpl ; wj iW ; qj ¼ hql ; /j iW . Inserting (26) into (25), and using the orthogonality of the POD bases l /j j¼1 , we obtain for i ¼ 1; ; l the systems
p_ i
* ! !2 + * !3 + l l l l X X X X ¼ qj hA/j ; wi iW c qj /j pj wj ; wi c qj /j ; wi ; j¼1
q_ i
l wj j¼1 and
¼
j¼1
*
l X
pj hAwj ; /i iW þ c
j¼1
l X
j¼1
!
pj wj
j¼1
l X
!2
qj /j
j¼1
+
W
j¼1
*
þc
; /i W
l X
pj wj
j¼1
!3
+ ; /i
W
: W
After defining U ¼ ½/1 ; /2 ; ; /l 2 Rml ; W ¼ ½w1 ; w2 ; ; wl 2 Rml ; ðBÞij ¼ hA/j ; wi iW , we obtain
p_ ¼ Bq cWT ðUqÞ ðWpÞ2 cWT ðUqÞ3 ; q_ ¼ BT p þ cUT ðWpÞ ðUqÞ2 þ cUT ðWpÞ3
ð27Þ
with both ’’ operation and the powers are hold elementwise. The reduced NLS Eq. (27) is also Hamiltonian and is solved, as the unreduced semi-discretized NLS Eq. (1), with the symplectic midpoint method applying linear-nonlinear Strang splitting [6]: In order to solve (25) efficiently, we apply the second order linear, non-linear Strang splitting [6]
iut ¼ N u þ Lu;
Lu ¼ uxx ;
N u ¼ cjuj2 u:
The nonlinear parts of the equations are solved by Newton–Raphson method. In the numerical examples, the boundary conditions are periodic, so that the resulting discretized matrices are circulant. For solving the linear system of equations, we have used the Matlab toolbox smt [9], which is designed for solving linear systems with a structured coefficient matrix like the circulant and Toeplitz matrices. It reduces the number of floating point operations for matrix factorization to O ðn log nÞ. 5. Numerical results All weights in the POD approximation are taken equally as ai ¼ 1=n and W ¼ I. Then the average ROM error, difference between the numerical solutions of NLS equation and ROM is measured in the form of the error between the fully discrete NLS solution
ROM error ¼
!1=2 n 1X jjy ðt j Þ yl ðtj Þjj : n j¼1 h
The average Hamiltonian ROM error is given by n 1X ðHh ðt j Þ Hl ðtj ÞÞ2 n j¼1
!1=2 ;
where Hh ðtj Þ and Hl ðtj Þ refer to the discrete Hamiltonian errors at the time instance t j corresponding to the full-order and ROM solutions, respectively. The energy of the Hamiltonian PDEs is usually expressed by the Hamiltonian. It is well known that symplectic integrators like the midpoint rule can preserve the only quadratic Hamiltonians exactly. Higher order polynomials and nonlinear Hamiltonians are preserved by the symplectic integration approximately, i.e. the approximate Hamiltonians do not show any drift in long term integration. For large matrices, the SVD is very time consuming. Recently several randomized methods are developed [10], which are very efficient when the rank is very small, i.e, d minðm; nÞ. We compare the efficiency of MATLAB programs svd and fsvd (based on the algorithm in [10]) for computation of singular values for the NLS equations in this section, on a PC with AMD FX(tm)-8150 Eight-Core Processor and 32 Gb RAM. The accuracy of the SVD is measured by L2 norm, jjY U RV T jjW . The randomized version of SVD, the fast SVD fsvd, requires the rank of the matrix as input parameter, which can be determined by MATLAB’s rank routine. When the singular values decay rapidly and the size of the matrices is very large, then randomized methods [10] are more efficient than MATLAB’s svd. Computation of the rank with rank and singular values with fsvd requires much less time than the svd for one and two dimensional NLS equations (Table 1).
515
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519 Table 1 Comparison of svd and fsvd. Problem
1D NLS 2D NLS CNLS
Size of snapshot matrix
Column rank
Computing time rank routine
fsvd Computing time
Accuracy
svd Computing time
Accuracy
32 x 50001 6400 x 30001 128 x 2001
15 25 122
0.14 279.74 0.05
0.73 3.81 0.10
6.4e-14 2.0e-13 2.7e-14
194.48 1300.42 1.19
1.3e-16 2.5e-15 1.5e-16
5.1. One-dimensional NLS equation For the one-dimensional NLS Eq. (1), we have taken the example in [2] with c ¼ 2 and the periodic boundary conditions in pffiffiffi the interval ½L=2; L=2 with L ¼ 2 2p. The initial conditions are pðx; 0Þ ¼ 0:5ð1 þ 0:01 cosð2px=LÞÞ and qðx; 0Þ ¼ 0. As mesh sizes in space and time, Dx ¼ L=32 and Dt ¼ 0:01 are used, respectively. Time steps are bounded by the stability condition for 2
the splitting method [6]; Dt < 2DLx where L is the period of the problem. The discretized Hamiltonian is given by (24) with c ¼ 2. The singular values of the snapshot matrix are rapidly decaying (Fig. 1) so that only few POD modes would be sufficient to approximate the fully discretized NLS equation. In all figures we have plotted the normalized singular values of the real part of the solutions. The decay of the singular values of the imaginary part is similar. In Fig. 2, the relative errors are plotted. As expected with increasing number of POD basis functions l, the errors in the energy and the errors between the discrete solutions of the fully discretized NLS equation and the reduced order model decreases which confirm the error analysis given in Section 3. In Figs. 3 and 4, the evolution of the Hamiltonian error and the numerical solution at time t ¼ 500 are shown for the POD basis with l ¼ 4, where 99.99 % of the energy of the system is well preserved. These figures confirm that the reduced model well preserves the Hamiltonian, and the numerical solution is close to the fully discrete solution. 5.2. Two-dimensional NLS equation We consider the following two-dimensional NLS equation [12]
iWt þ Wxx þ Wyy þ jWj2 W ¼ 0
on ½0; 2p ½0; 2p
with the exact solution, Wðx; y; tÞ ¼ expðiðx þ y tÞÞ. The mesh size for spatial discretization and time step size are taken as Dx ¼ Dy ¼ 2p=80 and Dt ¼ 0:001, respectively. The discrete Hamiltonian is given by
! m 2 X piþ1;j pi;j 2 qiþ1;j qi;j 2 pi;jþ1 pi;j 2 qi;jþ1 qi;j 2 1 1 2 H ¼ DxDy þ þ þ pi;j þ q2i;j : 2 4 Dx Dx Dy Dy i;j¼1 Only 3 POD modes were sufficient to capture almost all of the energy of the system (Table 2). A comparison of the Hamiltonian errors in long term computation shows that the reduced order model with a few POD modes preserve the energy of the system very well (Fig. 6). The singular values of 2D NLS are decreasing not continuously as for 1D NLS equation (Fig. 5). 5.3. Coupled NLS equation We consider two coupled NLS equations (CNLS) with elliptic polarization with plane wave solutions [1] 0
10
−2
10
−4
10
−6
10
−8
10
−10
10
−12
10
0
5
10
Fig. 1. 1D NLS, Decay of the normalized singular values.
15
516
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519 −6
x 10 4 energy−error of POD model
data−error of POD model
0.2 0.15 0.1 0.05
3
2
1
0 0 −0.05 0
2
4
6
8
10
12
14
16
0
2
4
# of POD modes used
6
8
10
12
14
16
# of POD modes used
Fig. 2. 1D NLS, Decay of the ROM errors: solution (left), Hamiltonian (right).
−5
−5
x 10
2 0 Energy Error by Rom
Energy Error
0
−5
−10
−15
0
x 10
−2 −4 −6 −8 −10 −12
100
200
300
400
500
−14 0
100
200
time
300
400
time
Fig. 3. 1D NLS, Energy error: full-order model (left), ROM with 4 POD modes (right).
Fig. 4. 1D NLS, Envelope of the approximate solution jwj: full-order model (left), ROM with 4 POD modes (right).
Table 2 2D NLS, RIC and errors for the real part of the solution. # POD
Info (%)
ROM Hamiltonian error
ROM error
1 2
51.65 99.99
8.181e-002 6.116e-007
2.770e + 001 1.040e-003
500
517
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519 0
10
−2
10
−4
10
−6
10
−8
10
−10
10
−12
10
0
5
10
15
20
25
Fig. 5. 2D NLS, Decay of the normalized singular values.
−7
3.5
−4
x 10
2
3
1.5 Energy Error by Rom
Energy Error
2.5 2 1.5 1 0.5
1 0.5 0 −0.5 −1
0 −0.5 0
x 10
5
10
15 time
20
25
30
−1.5 0
5
10
15 time
20
25
30
Fig. 6. 2D NLS, Energy error: full-order model (left), ROM with 2 POD modes (right).
i
@ W1 @ 2 W1 þ þ ðjW1 j2 þ jW2 j2 ÞW1 ¼ 0; @t @x2
i
@ W2 @ 2 W2 þ þ ðjW2 j2 þ jW1 j2 ÞW2 ¼ 0; @t @x2
ð28Þ
using the initial conditions
W1 ðx; 0Þ ¼ ð0:5Þð1 0:1 cosð0:5xÞÞ;
W2 ðx; 0Þ ¼ ð0:5Þð1 0:1 cosð0:5xÞÞ:
The equations are solved over the space ½0; 8p and time interval [0, 100], respectively, with the mesh size and time steps dx ¼ 8p=128; Dt ¼ 0:05. The discrete Hamiltonian is given as [1]
0 ! !2 !2 !2 1 1 1 2 q1jþ1 q1j p2jþ1 p2j q2jþ1 q2j 1 @ pjþ1 pj A H ¼ Dx þ þ þ 2 Dx Dx Dx Dx j¼1
1 1 1 2 2 2 2 2 2 2 2 2 2 ððp1j Þ þ ðp2j Þ Þ þ ððq1j Þ þ ðq2j Þ Þ þ ððpj Þ þ ðp2j Þ Þððq1j Þ þ ðq2j Þ Þ ; þ 4 2 m X
where p1 ; q1 and p2 ; q2 denote the real and imaginary parts of w1 and w2 , respectively. Table 3 Coupled NLS, RIC and errors for the real part of W1 . #POD
RIC(%)
ROM Hamilton error
ROM error
2 3 4 5
99.58 99.98 99.99 99.99
1.879e-004 1.865e-004 1.213e-004 2.825e-005
5.060e-001 3.761e-001 6.491e-002 3.919e-003
518
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519 0
10
−5
10
−10
10
−15
10
0
20
40
60
80
100
120
140
Fig. 7. Coupled NLS, Decay of the normalized singular values for the real part of W1 .
0.035
0.03
0.03
0.025 Energy Error by Rom
Energy Error
0.025 0.02 0.015 0.01 0.005
0.015 0.01 0.005 0
0 −0.005 0
0.02
20
40
60 time
80
100
−0.005 0
20
40
60
80
100
time
Fig. 8. Coupled NLS, Hamiltonian error: full-order model (left) and ROM with 5 POD modes (right).
Fig. 9. Coupled NLS, Interaction of solitons jW1 j and jW2 j: full-order model (left) and ROM with 5 POD modes (right).
Figs. 8 and 9 and Table 3 show that only few POD modes are necessary to capture the dynamics of the CNLS equation. The singular values are decreasing not so rapidly (Fig. 7) as in case of single 1D and 2D NLS equations. 6. Conclusions A reduced model is derived for the NLS equation by preserving the Hamiltonian structure. A priori error estimates are obtained for the mid-point rule as time integrator for the reduced dynamical system. Numerical results show that the energy
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
519
and the phase space structure of the three different NLS equations are well preserved by using few POD modes. The number of the POD modes containing most of the energy depends on the decay of the singular values of the snapshot matrix, reflecting the dynamics of the underlying systems. In a future work, we will investigate the dependence of the ROM solutions on parameters for the CNLS equation by performing a sensitivity analysis. Acknowledgements The authors wish to thank the referees for their helpful suggestions and comments. This research was supported by the Middle East Technical University Scientific Research Fund (Project: BAP-07–05-2012–102). References [1] A. Aydn, B. Karasözen, Symplectic and multi-symplectic methods for coupled nonlinear Schrödinger equations with periodic solutions, Comput. Phys. Commun. 177 (2007) 66–583.. [2] A.L. Islas, D.A. Karpeev, C.M. Schober, Geometric integrators for the nonlinear Schrödinger equation, J. Comput. Phys. 173 (2001) 16–148. . [3] A. Studingeer, S. Volkwein, Numerical analysis of POD a-posteriori error estimation for optimal control, in Control and Optimization with PDE Constraints, eds: K. Bredies, C. Clason, K. Kunisch, G. von Winckel***, International Series of Numerical Mathematics 164 (2013) 137–158. [4] E. Schlizerman, E. Ding, O.M. Williams, J.N. Kutz, The Proper orthogonal decomposition for dimensionality reduction in mode-locked lasers and optical systems, Int. J. Optics. (2012) 31604. . [5] G. Berkooz, P. Holmes, J.L. Lumley, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, Cambridge Monographs on Mechanics, 1996. [6] J.A.C. Weideman, B.M. Herbst, Split-step methods for the solution of the nonlinear Schrödinger equation, SIAM J. Numer. Anal. 23 (1986) 485–507. . [7] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems, Numer. Math. 90 (2001) 17–148. . [8] O. Lass, S. Volkwein, Adaptive POD basis computation for parametrized nonlinear systems using optimal snapshot location, Comput. Optim. Appl. 58 (2014) 645–677. . [9] M. Redivo-Zaglia, G. Rodriguez, SMT: a Matlab structured matrices toolbox, Numer. Algorithms 59 (2012) 639–659. . [10] N. Halko, P.G. Martinsson, Y. Shkolnisky, M. Tygert, An algorithm for the principal component analysis of large data sets, SIAM J. Sci. Comput. 33 (2011) 2580–2594. . [11] S. Volkwein, Proper Orthogonal Decomposition: Theory and Reduced-Order Modelling, Lecture Notes, University of Konstanz, Department of Mathematics and Statistics, , 2013 [12] Ya-Ming Chen, Zhu Hua-Jun Zhu, He Song, Multi-symplectic splitting method for two-dimensional nonlinear Schrödinger equation, Commun. Theor. Phys. 56 (2011) 617–622. .
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Model order reduction for nonlinear Schrödinger equation Bülent Karasözen a,⇑, Canan Akkoyunlu b, Murat Uzunca c a b c
Department of Mathematics and Institute of Applied Mathematics, Middle East Technical University, 06800 Ankara, Turkey Department of Mathematics and Computer Sciences, Istanbul Kültür University, 34156 Istanbul, Turkey Department of Mathematics, Middle East Technical University, 06800 Ankara, Turkey
a r t i c l e
i n f o
Keywords: Nonlinear Schrödinger equation Proper orthogonal decomposition Model order reduction Error analysis
a b s t r a c t We apply the proper orthogonal decomposition (POD) to the nonlinear Schrödinger (NLS) equation to derive a reduced order model. The NLS equation is discretized in space by finite differences and is solved in time by structure preserving symplectic mid-point rule. A priori error estimates are derived for the POD reduced dynamical system. Numerical results for one and two dimensional NLS equations, coupled NLS equation with soliton solutions show that the low-dimensional approximations obtained by POD reproduce very well the characteristic dynamics of the system, such as preservation of energy and the solutions. Ó 2015 Elsevier Inc. All rights reserved.
1. Introduction The nonlinear Schrödinger (NLS) equation arises as the model equation with second order dispersion and cubic nonlinearity describing the dynamics of slowly varying wave packets in nonlinear fiber optics, in water waves and in Bose–Einstein condensate theory. We consider the NLS equation
iWt þ Wxx þ cjWj2 W ¼ 0;
ð1Þ
with the periodic boundary conditions Wðx þ L; tÞ ¼ Wðx; tÞ. Here W ¼ Wðx; tÞ is a complex valued function, c is a parameter pffiffiffiffiffiffiffi and i ¼ 1. The NLS equation is called ‘‘focusing’’ if c > 0 and ‘‘defocusing’’ if c < 0; for c ¼ 0, it reduces to the linear Schrödinger equation. In last two decades, various numerical methods were applied for solving NLS equation, among them are the well-known symplectic and multisymplectic integrators and discontinuous Galerkin methods. There is a strong need for model order reduction techniques to reduce the computational costs and storage requirements in large scale simulations, yielding low-dimensional approximations for the full high-dimensional dynamical system, which reproduce the characteristic dynamics of the system. Among the model order reduction techniques, the proper orthogonal decomposition (POD) is one of the most widely used method. It was first introduced for analyzing coherent structures and turbulent flow in numerical simulation of fluid dynamics equations [5]. It has been successfully used in different fields including signal processing, fluid dynamics, parameter estimation, control theory and optimal control of partial differential equations. In this paper, we apply the POD to the NLS equation. To the best of our knowledge, there is only one paper where POD is applied to NLS equation [4], where only one and two modes approximations of the NLS equation are used in the Fourier domain in connection with mode-locking ultra short laser applications. In this paper, the NLS equation being a semi-linear partial differential equation (PDE) is discretized in space and time by preserving the symplectic structure and
⇑ Corresponding author. E-mail addresses: [email protected] (B. Karasözen), [email protected] (C. Akkoyunlu), [email protected] (M. Uzunca). http://dx.doi.org/10.1016/j.amc.2015.02.001 0096-3003/Ó 2015 Elsevier Inc. All rights reserved.
510
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
the energy (Hamiltonian). Then, from the snapshots of the fully discretized dynamical system, the POD basis functions are computed using the singular value decomposition (SVD). The reduced model consists of Hamiltonian ordinary differential equations (ODEs), which indicates that the geometric structure of the original system is preserved for the reduced model. The semi-disretized NLS equations and the reduced equations are solved in time using Strang splitting and mid-point rule. A priori error estimates are derived for POD reduced model, which is solved by mid-point rule. It turns out that most of the energy of the system can be accurately approximated by using few POD modes. Numerical results for a NLS equation with soliton solutions confirm that the energy of the system is well preserved by POD approximation and the solution of the reduced model are close to the solution of the fully discretized system. The paper is organized as follows. In Section 2, the POD method and its application to semi-linear dynamical systems are reviewed. In Section 3, a priori error estimators are derived for the mid-point time-discretization of semi-linear PDEs. Numerical solution of the semi-discrete NLS equation and the POD reduced form are described in Section 4. In the last section, Section 5, the numerical results for the reduced order models of NLS equations are presented. 2. The POD approximation for semi-linear PDEs In the following, we briefly describe the important features of the POD reduced order modeling (ROM); more details can be found in [11]. In the first step of the POD based model order reduction, the set of snapshots, the discrete solutions of the nonlinear PDE, are collected. The snapshots are usually equally spaced in time corresponding to the solution of PDE obtained by finite difference or finite element method. The snapshots are then used to determine the POD bases which are much smaller than the snapshot set. In the last step, the POD reduced order model is constructed to obtain approximate solutions of the PDE. We mention that the choice of the snapshots representing the dynamics of the underlying PDE is crucial for the effectiveness of POD based reduced model. Let X be a real Hilbert space endowed with inner product h; iX and norm kkX . For y1 ; . . . ; yn 2 X, we set n V ¼ spanfy1 ; ; yn g, as the ensemble consisting of the snapshots yj j¼1 . In the finite difference context, the snapshots can be viewed as discrete solutions yj 2 Rm at time instances tj ; j ¼ 1; . . . ; n, and ½y1 ; . . . ; yn 2 Rmn denotes the snapshot matrix. Let fwk gdk¼1 denote an orthonormal basis of V of dimension d. Then, any yj 2 V can be expressed as
yj ¼
d X hyj ; wk iX wk ;
j ¼ 1; . . . ; n:
ð2Þ
k¼1
The POD is constructed by choosing the orthonormal basis such that for every l 2 f1; . . . ; dg, the mean square error between the elements yj ; 1 6 j 6 n, and the corresponding l th partial sum of (2) is minimized on average:
2 n l X X ~ ~ a y hy ; i u u j k X k ; j j ~ 1 ;...;u ~ l 2X u j¼1 k¼1
~i ; u ~ j iX ¼ dij ; hu
min
1 6 i; j 6 l;
ð3Þ
X
where aj ’s are non-negative weights. Throughout this paper, we take the space X ¼ Rm endowed with the weighted inner product hu; v iW ¼ uT W v with the diagonal elements of the diagonal matrix W, and also aj ’s are the trapezoidal weights so that we obtain all the computations in L2 -sense. Under these choices, the solution of the above minimization problem is given by the following theorem: Theorem 1 [11, Theorem 1.3.2]. Let Y ¼ ½y1 ; . . . ; yn 2 Rmn be a given snapshot matrix with rank d 6 minfm; ng. For a b ¼ W 1=2 YD1=2 2 Rmn . Further, symmetric, positive definite matrix W 2 Rmm and a diagonal matrix D ¼ diagfa1 ; . . . ; an g, set Y b¼U b RV b T be the SVD of Y b , where U b ¼ ½u b ¼ ½v ^1 ; . . . ; u ^ m 2 Rmm ; V ^1; . . . ; v ^ n 2 Rnn are orthogonal matrices and the matrix let Y mn b . Then, for any R2R is all zero but first d diagonal elements are the nonzero singular values, r1 P r2 P . . . P rd , of Y l 2 f1; . . . ; dg, the solution to
2 n l X X ~ k iW u ~k min m aj yj hyj ; u ; ~ 1 ;...;u ~ l 2R u j¼1 k¼1
~i ; u ~ j iW ¼ dij ; hu
1 6 i; j 6 l;
ð4Þ
W
l
^ i gi¼1 . Moreover, we have r2i ¼ ki is given by the orthonormal, with respect to the weighted inner product, vectors fwi gli¼1 ¼ fW 1=2 u bT Y b. where ki ’s are the eigenvalues of the correlation matrix Y We consider the following initial value problem for POD approximation
_ yðtÞ ¼ AyðtÞ þ f ðt; yðtÞÞ; m
m
t 2 ½0; T;
yð0Þ ¼ y0 ;
ð5Þ
where f : ½0; T R ! R is continuous in both arguments and locally Lipschitz-continuous with respect to the second argument. The semi-discrete form of NLS Eq. (1) is a semi-linear equation as (5) where the cubic nonlinear part is locally Lipschitz
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
511
l continuous. Suppose that we have determined a POD basis wj j¼1 of rank l 2 f1; . . . ; dg in Rm , then we make the ansatz
yl ðtÞ ¼
l X hyl ðtÞ; wj iW wj ; |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} j¼1
t 2 ½0; T:
ð6Þ
¼:ylj ðtÞ
Substituting (6) in (5), we obtain the reduced model l l X X ylj ðtÞAwj þ f ðt; yl ðtÞÞ; y_ lj ðtÞwj ¼ j¼1
t 2 ½0; T;
l X ylj ð0Þwj ¼ y0 :
j¼1
ð7Þ
j¼1
The POD approximation (7) holds after projection on the l dimensional subspace V l ¼ spanfw1 ; . . . ; wl g. From (7) and hwj ; wi iW ¼ dij , we get
y_ li ðtÞ ¼
l X ylj ðtÞhAwj ; wi iW þ hf ðt; yl ðtÞÞ; wi iW
ð8Þ
j¼1
for 1 6 i 6 l and t 2 ð0; T. Let us introduce the matrix
B ¼ fbij g 2 Rll ;
bij ¼ hAwj ; wi iW ;
the non-linearity F ¼ ðF 1 ; ; F l ÞT : ½0; T Rl ! Rl by
* F i ðt; yÞ ¼
f t;
l X
! yj wj ; wi
j¼1 l
and the vector y ¼
+
T ðyl1 ; . . . ; yll Þ
y_ l ðtÞ ¼ Byl ðtÞ þ Fðt; yl ðtÞÞ;
;
t 2 ½0; T;
y ¼ ðy1 ; ; yl Þ 2 Rl ;
W
: ½0; T ! Rl . Then, (8) can be expressed as
t 2 ð0; T:
ð9Þ
The initial condition of the reduced system is given by yl ð0Þ ¼ y0 with
y0 ¼ ðhy0 ; w1 iW ; . . . ; hy0 ; wl iW ÞT 2 Rl : The system (9) is called the POD-Galerkin projection for (5). The ROM is constructed with POD basis vectors fwi gli¼1 of rank l. In case of l d, the l-dimensional reduced system (9) is a low-dimensional approximation for (5). The POD basis can also be computed using eigenvalues and eigenvectors. We prefer singular value decomposition, because it is more accurate than the computation of the eigenvalues. The singular values decay up to machine precision, where the eigenvalues stagnate several orders above due the fact ki ¼ r2i [3]. The choice of l is based on heuristic considb which is erations combined with observing the ratio of the modeled energy to the total energy contained in the system Y expressed by the relative information content (RIC)
EðlÞ ¼
l X
r2i 100;
i¼1
b are normalized, i.e., Pd r2 ¼ 1 [8]. The total energy of the when the squared singular values of the snapshot matrix Y i¼1 i system is contained in a small number of POD modes. In practice, l is chosen by guaranteeing that EðlÞ capturing at least % 99 of total energy of the system. 3. POD error analysis for the mid-point rule In this section we derive a priori error estimates for the mid-point rule. A priori error estimates for POD method are obtained for linear and semi-linear parabolic equations using backward Euler and Crank–Nicholson (trapezoidal rule) time discretization in [7]. When implicit midpoint rule is applied for solving the reduced model (9), we obtain the discrete system for the sequence fY j gnj¼1 in V ln ¼ spanfwn1 ; . . . ; wnl g (l 6 d)
Y j Y j1 n 1 Y j þ Y j1 ; wni ; wi ¼ AðY j þ Y j1 Þ þ f t; ; 2 Dt 2 W W hY 1 ; wni iW ¼ hy0 ; wni iW ;
ð10Þ
i ¼ 1; . . . ; l;
where Y j denotes an approximation for yl at the time t j . We are interested in estimating use the following decomposition for the error
ð11Þ 2 m j¼1 aj yðt j Þ Y j W . For u 2 R . We
Pn
512
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
yðtj Þ Y j ¼ yðtj Þ Pln yðt j Þ þ Pln yðtj Þ Y j ¼ .lj þ #lj ; where
.lj ¼ yðtj Þ Pln yðtj Þ; #lj ¼ Pln yðtj Þ Y j , and Pln : Rm ! V ln is the projection given by Pln u ¼
l P n
l X hu; wni iW wni ;
W
i¼1
¼ 1:
l
Using that fwni gi¼1 is the POD basis of rank l, we have the estimate for the terms involving
.lj
2 n l n 2 n 2 d X X X X X n n aj yðtj Þ < yðtj Þ; wi > wi ¼ aj yðtj Þ P ln yðtj Þ ¼ aj .lj ¼ r2i : W W j¼1 i¼1 j¼1 j¼1 i¼lþ1
ð12Þ
W
l ¼ ð#l #l Þ=Dt, we obtain Next, we estimate the terms involving #lj . Using the notation @# j j j1
yðtj Þ yðtj1 Þ Y j Y j1 n Pln ; wi Dt Dt W
Dt Y j þ Y j1 Dt Y j þ Y j1 þ f y tj f ; wni ¼ A y tj þ hwlj þ zlj ; wni iW ; 2 2 2 2 W
l ; wn i ¼ h@# i W j
ð13Þ
where
zlj ¼ Pln
yðt j Þ yðt j1 Þ yðt j Þ yðt j1 Þ ; Dt Dt
wlj ¼
yðt j Þ yðt j1 Þ Dt : y_ tj 2 Dt
Choosing wni ¼ #lj þ #lj1 , and using Lipschitz-continuity of f and the Cauchy–Schwarz inequality in (13) and noting that
2
2 l ; #l þ #l i ¼ 1 h@# #lj #lj1 : j j j1 W W W Dt we arrive at
2 l #j
W
2 Y þY 6 #lj1 þ DtðkAkW þ Lf Þyðt j D2tÞ j 2 j1 #lj þ #lj1 W W W l l l l þ Dt wj þ zj #j þ #j1 : W
W
ð14Þ
W
Using the Taylor series approximation yðt j D2tÞ ðyðt j Þ þ yðt j1 ÞÞ=2, applying Cauchy–Schwarz and Young’s inequalities, and collecting the terms yields
2
2 þ kzlj k2W þ kwlj k2W ; ð1 c1 DtÞk#lj k2W 6 ð1 þ c1 DtÞk#lj1 k2W þ Dt c2 .lj þ .lj1 W
W
ð15Þ
where the positive constants c1 and c2 depend on kAkW , Lipschitz constant Lf and a constant C Tay related to the Taylor series approximation. For 0 < Dt 6 2c11 , using the fact that Dtj 6 T, we have
1 6 1 þ 2c1 Dt 1 c 1 Dt
and ð1 þ 2c1 DtÞ j 6 e2c1 T ; ð1 þ c1 DtÞ j 6 ec1 T :
ð16Þ
Summation on j in (15) with using (16) and Cauchy–Schwarz inequality yields, j 2 X l 2 l . þ .l 2 þ kzl k2 þ kwl k2 ; #j 6 C Dt2 k W k W k W k1 W W
ð17Þ
k¼1
where the positive constant C depend on c1 ; c2 ; T and j. Next, we estimate the term involving wlk :
Dt 2
2 Z T j j X X l 2 e 2 Dt 4 w ¼ Dt2 yðt k Þ yðt k1 Þ y_ t k Dt 6 C kyttt ðtÞk2 dt k W 2 D t 0 W k¼1 k¼1
e depending on y, but independent of n. for a constant C
ð18Þ
513
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
_ k Þ: Now, we estimate the term involving zlk . Adding and subtracting Pln yðt
2
l 2 yðt k Þ yðt k1 Þ l z ¼ Pl yðtk Þ yðtk1 Þ Pl yðt _ _ Þ þ P Þ yðt k k n n n k W Dt Dt W 2 yðt Þ yðt Þ k k1 2 þ 4kPl yðt _ kÞ _ k Þ yðt _ k Þk2W 6 ð2kPln kW þ 4Þ n yðt Dt W 2
D t yðt Þ yðt Þ k k1 þ 4kPl yðt _ k Þ yðt _ k Þk2W þ C Tay Dt2 6 6 n y_ tk 2 Dt W _ k Þ yðt _ k Þk2W þ C Tay Dt 2 ; ¼ 6kwlk k2W þ 4kPln yðt
ð19Þ
€ðnÞ for some n 2 ðtk D2t ; t k Þ. Note that ak ’s are the trapezoidal weights which satisfy Dt 6 ak , and combining of where C Tay ¼ y (18) and (19) leads
Dt 2
j n X X b Dt 4 ; _ k Þ yðt _ k Þk2W þ C kzlk k2W þ kwlk k2W 6 8 ak kPln yðt k¼1
ð20Þ
k¼1
b depends on C e ; ky k 2 where C ttt L ð0;T;Rm Þ and C Tay . Imposing the estimate (20) in (17), we obtain
2 n n X X l 2 b Dt 4 : _ k Þ yðt _ k Þk2W þ C C 8ak kPln yðt #j 6 4C ak k.lk kW þ C W
k¼1
ð21Þ
k¼1
In addition, we have that
Pn
k¼1
estimate to the term involving
#lj
n X
d X
j¼1
i¼lþ1
aj k#lj k2W 6 C Dt4 þ
Pd
_ j Þ yðt _ j Þk2W ¼ ak 6 T and kPln yðt
2 n _ i¼lþ1 jhyðt j Þ; wi iW j .
Using these identities, we arrive at the
as
r2i þ
!!
n X
_ j Þ; wni iW j2 aj jhyðt
ð22Þ
;
j¼1
where C is dependent on y; T, but independent of n and l. Finally, combining the estimates (12) and (22), we obtain the error estimate n X
n X
n X
n X
d X
j¼1
j¼1
j¼1
j¼1
i¼lþ1
aj kyðtj Þ Y j k2W ¼
aj k#lj þ .lj k2W 6 2
aj k#lj k2W þ 2
aj k.lj k2W 6 C E
2r2i þ
n X
!
!
_ j Þ; wni iW j2 þ Dt 4 ; aj jhyðt
j¼1
where C E is dependent on y; T, but independent of n and l. As for the backward Euler and Crank–Nicolson method [7], the error between the reduced and the unreduced solutions depend for the mid-point rule on the time discretization and on the number of not modelled POD snapshots. 4. Discretization of NLS equation One dimensional NLS Eq. (1) can be written by decomposing W ¼ p þ iq in real and imaginary components
pt ¼ qxx cðp2 þ q2 Þq;
qt ¼ pxx þ cðp2 þ q2 Þp;
ð23Þ
as an infinite dimensional Hamiltonian PDE in the phase space u ¼ ðp; qÞ
u_ ¼ D
dH ; du
H¼
Z
1 2 c 2 px þ q2x ðp2 þ q2 Þ dx; 2 2
D¼
0
1
1 0
T
:
After discretizing the Hamiltonian in space by finite differences,
H ¼ Dx
!
n X pjþ1 pj 2 qjþ1 qj 2 c 2 1 2 þ ðpj þ q2j Þ ; 2 Dx Dx 2 j¼1
ð24Þ
we obtain the semi-discretized Hamiltonian ode’s
pt ¼ Aq cqðp2 þ q2 Þ;
qt ¼ Ap þ cpðp2 þ q2 Þ;
where A is the circulant matrix
0 2 B 1 B B A¼B B @ 1
1 1 1 C 2 1 C .. .. .. C C: . . . C 1 2 1 A 1 2
ð25Þ
514
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
4.1. Reduced order model for NLS equation l l Suppose that we have determined POD bases wj j¼1 and /j j¼1 of rank l ¼ f1; . . . ; dg in Rm . Then we make the ansatz
pl ¼
l X
pj ðtÞwj ðxÞ;
ql ¼
j¼1
l X qj ðtÞ/j ðxÞ;
ð26Þ
j¼1
where pj ¼ hpl ; wj iW ; qj ¼ hql ; /j iW . Inserting (26) into (25), and using the orthogonality of the POD bases l /j j¼1 , we obtain for i ¼ 1; ; l the systems
p_ i
* ! !2 + * !3 + l l l l X X X X ¼ qj hA/j ; wi iW c qj /j pj wj ; wi c qj /j ; wi ; j¼1
q_ i
l wj j¼1 and
¼
j¼1
*
l X
pj hAwj ; /i iW þ c
j¼1
l X
j¼1
!
pj wj
j¼1
l X
!2
qj /j
j¼1
+
W
j¼1
*
þc
; /i W
l X
pj wj
j¼1
!3
+ ; /i
W
: W
After defining U ¼ ½/1 ; /2 ; ; /l 2 Rml ; W ¼ ½w1 ; w2 ; ; wl 2 Rml ; ðBÞij ¼ hA/j ; wi iW , we obtain
p_ ¼ Bq cWT ðUqÞ ðWpÞ2 cWT ðUqÞ3 ; q_ ¼ BT p þ cUT ðWpÞ ðUqÞ2 þ cUT ðWpÞ3
ð27Þ
with both ’’ operation and the powers are hold elementwise. The reduced NLS Eq. (27) is also Hamiltonian and is solved, as the unreduced semi-discretized NLS Eq. (1), with the symplectic midpoint method applying linear-nonlinear Strang splitting [6]: In order to solve (25) efficiently, we apply the second order linear, non-linear Strang splitting [6]
iut ¼ N u þ Lu;
Lu ¼ uxx ;
N u ¼ cjuj2 u:
The nonlinear parts of the equations are solved by Newton–Raphson method. In the numerical examples, the boundary conditions are periodic, so that the resulting discretized matrices are circulant. For solving the linear system of equations, we have used the Matlab toolbox smt [9], which is designed for solving linear systems with a structured coefficient matrix like the circulant and Toeplitz matrices. It reduces the number of floating point operations for matrix factorization to O ðn log nÞ. 5. Numerical results All weights in the POD approximation are taken equally as ai ¼ 1=n and W ¼ I. Then the average ROM error, difference between the numerical solutions of NLS equation and ROM is measured in the form of the error between the fully discrete NLS solution
ROM error ¼
!1=2 n 1X jjy ðt j Þ yl ðtj Þjj : n j¼1 h
The average Hamiltonian ROM error is given by n 1X ðHh ðt j Þ Hl ðtj ÞÞ2 n j¼1
!1=2 ;
where Hh ðtj Þ and Hl ðtj Þ refer to the discrete Hamiltonian errors at the time instance t j corresponding to the full-order and ROM solutions, respectively. The energy of the Hamiltonian PDEs is usually expressed by the Hamiltonian. It is well known that symplectic integrators like the midpoint rule can preserve the only quadratic Hamiltonians exactly. Higher order polynomials and nonlinear Hamiltonians are preserved by the symplectic integration approximately, i.e. the approximate Hamiltonians do not show any drift in long term integration. For large matrices, the SVD is very time consuming. Recently several randomized methods are developed [10], which are very efficient when the rank is very small, i.e, d minðm; nÞ. We compare the efficiency of MATLAB programs svd and fsvd (based on the algorithm in [10]) for computation of singular values for the NLS equations in this section, on a PC with AMD FX(tm)-8150 Eight-Core Processor and 32 Gb RAM. The accuracy of the SVD is measured by L2 norm, jjY U RV T jjW . The randomized version of SVD, the fast SVD fsvd, requires the rank of the matrix as input parameter, which can be determined by MATLAB’s rank routine. When the singular values decay rapidly and the size of the matrices is very large, then randomized methods [10] are more efficient than MATLAB’s svd. Computation of the rank with rank and singular values with fsvd requires much less time than the svd for one and two dimensional NLS equations (Table 1).
515
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519 Table 1 Comparison of svd and fsvd. Problem
1D NLS 2D NLS CNLS
Size of snapshot matrix
Column rank
Computing time rank routine
fsvd Computing time
Accuracy
svd Computing time
Accuracy
32 x 50001 6400 x 30001 128 x 2001
15 25 122
0.14 279.74 0.05
0.73 3.81 0.10
6.4e-14 2.0e-13 2.7e-14
194.48 1300.42 1.19
1.3e-16 2.5e-15 1.5e-16
5.1. One-dimensional NLS equation For the one-dimensional NLS Eq. (1), we have taken the example in [2] with c ¼ 2 and the periodic boundary conditions in pffiffiffi the interval ½L=2; L=2 with L ¼ 2 2p. The initial conditions are pðx; 0Þ ¼ 0:5ð1 þ 0:01 cosð2px=LÞÞ and qðx; 0Þ ¼ 0. As mesh sizes in space and time, Dx ¼ L=32 and Dt ¼ 0:01 are used, respectively. Time steps are bounded by the stability condition for 2
the splitting method [6]; Dt < 2DLx where L is the period of the problem. The discretized Hamiltonian is given by (24) with c ¼ 2. The singular values of the snapshot matrix are rapidly decaying (Fig. 1) so that only few POD modes would be sufficient to approximate the fully discretized NLS equation. In all figures we have plotted the normalized singular values of the real part of the solutions. The decay of the singular values of the imaginary part is similar. In Fig. 2, the relative errors are plotted. As expected with increasing number of POD basis functions l, the errors in the energy and the errors between the discrete solutions of the fully discretized NLS equation and the reduced order model decreases which confirm the error analysis given in Section 3. In Figs. 3 and 4, the evolution of the Hamiltonian error and the numerical solution at time t ¼ 500 are shown for the POD basis with l ¼ 4, where 99.99 % of the energy of the system is well preserved. These figures confirm that the reduced model well preserves the Hamiltonian, and the numerical solution is close to the fully discrete solution. 5.2. Two-dimensional NLS equation We consider the following two-dimensional NLS equation [12]
iWt þ Wxx þ Wyy þ jWj2 W ¼ 0
on ½0; 2p ½0; 2p
with the exact solution, Wðx; y; tÞ ¼ expðiðx þ y tÞÞ. The mesh size for spatial discretization and time step size are taken as Dx ¼ Dy ¼ 2p=80 and Dt ¼ 0:001, respectively. The discrete Hamiltonian is given by
! m 2 X piþ1;j pi;j 2 qiþ1;j qi;j 2 pi;jþ1 pi;j 2 qi;jþ1 qi;j 2 1 1 2 H ¼ DxDy þ þ þ pi;j þ q2i;j : 2 4 Dx Dx Dy Dy i;j¼1 Only 3 POD modes were sufficient to capture almost all of the energy of the system (Table 2). A comparison of the Hamiltonian errors in long term computation shows that the reduced order model with a few POD modes preserve the energy of the system very well (Fig. 6). The singular values of 2D NLS are decreasing not continuously as for 1D NLS equation (Fig. 5). 5.3. Coupled NLS equation We consider two coupled NLS equations (CNLS) with elliptic polarization with plane wave solutions [1] 0
10
−2
10
−4
10
−6
10
−8
10
−10
10
−12
10
0
5
10
Fig. 1. 1D NLS, Decay of the normalized singular values.
15
516
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519 −6
x 10 4 energy−error of POD model
data−error of POD model
0.2 0.15 0.1 0.05
3
2
1
0 0 −0.05 0
2
4
6
8
10
12
14
16
0
2
4
# of POD modes used
6
8
10
12
14
16
# of POD modes used
Fig. 2. 1D NLS, Decay of the ROM errors: solution (left), Hamiltonian (right).
−5
−5
x 10
2 0 Energy Error by Rom
Energy Error
0
−5
−10
−15
0
x 10
−2 −4 −6 −8 −10 −12
100
200
300
400
500
−14 0
100
200
time
300
400
time
Fig. 3. 1D NLS, Energy error: full-order model (left), ROM with 4 POD modes (right).
Fig. 4. 1D NLS, Envelope of the approximate solution jwj: full-order model (left), ROM with 4 POD modes (right).
Table 2 2D NLS, RIC and errors for the real part of the solution. # POD
Info (%)
ROM Hamiltonian error
ROM error
1 2
51.65 99.99
8.181e-002 6.116e-007
2.770e + 001 1.040e-003
500
517
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519 0
10
−2
10
−4
10
−6
10
−8
10
−10
10
−12
10
0
5
10
15
20
25
Fig. 5. 2D NLS, Decay of the normalized singular values.
−7
3.5
−4
x 10
2
3
1.5 Energy Error by Rom
Energy Error
2.5 2 1.5 1 0.5
1 0.5 0 −0.5 −1
0 −0.5 0
x 10
5
10
15 time
20
25
30
−1.5 0
5
10
15 time
20
25
30
Fig. 6. 2D NLS, Energy error: full-order model (left), ROM with 2 POD modes (right).
i
@ W1 @ 2 W1 þ þ ðjW1 j2 þ jW2 j2 ÞW1 ¼ 0; @t @x2
i
@ W2 @ 2 W2 þ þ ðjW2 j2 þ jW1 j2 ÞW2 ¼ 0; @t @x2
ð28Þ
using the initial conditions
W1 ðx; 0Þ ¼ ð0:5Þð1 0:1 cosð0:5xÞÞ;
W2 ðx; 0Þ ¼ ð0:5Þð1 0:1 cosð0:5xÞÞ:
The equations are solved over the space ½0; 8p and time interval [0, 100], respectively, with the mesh size and time steps dx ¼ 8p=128; Dt ¼ 0:05. The discrete Hamiltonian is given as [1]
0 ! !2 !2 !2 1 1 1 2 q1jþ1 q1j p2jþ1 p2j q2jþ1 q2j 1 @ pjþ1 pj A H ¼ Dx þ þ þ 2 Dx Dx Dx Dx j¼1
1 1 1 2 2 2 2 2 2 2 2 2 2 ððp1j Þ þ ðp2j Þ Þ þ ððq1j Þ þ ðq2j Þ Þ þ ððpj Þ þ ðp2j Þ Þððq1j Þ þ ðq2j Þ Þ ; þ 4 2 m X
where p1 ; q1 and p2 ; q2 denote the real and imaginary parts of w1 and w2 , respectively. Table 3 Coupled NLS, RIC and errors for the real part of W1 . #POD
RIC(%)
ROM Hamilton error
ROM error
2 3 4 5
99.58 99.98 99.99 99.99
1.879e-004 1.865e-004 1.213e-004 2.825e-005
5.060e-001 3.761e-001 6.491e-002 3.919e-003
518
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519 0
10
−5
10
−10
10
−15
10
0
20
40
60
80
100
120
140
Fig. 7. Coupled NLS, Decay of the normalized singular values for the real part of W1 .
0.035
0.03
0.03
0.025 Energy Error by Rom
Energy Error
0.025 0.02 0.015 0.01 0.005
0.015 0.01 0.005 0
0 −0.005 0
0.02
20
40
60 time
80
100
−0.005 0
20
40
60
80
100
time
Fig. 8. Coupled NLS, Hamiltonian error: full-order model (left) and ROM with 5 POD modes (right).
Fig. 9. Coupled NLS, Interaction of solitons jW1 j and jW2 j: full-order model (left) and ROM with 5 POD modes (right).
Figs. 8 and 9 and Table 3 show that only few POD modes are necessary to capture the dynamics of the CNLS equation. The singular values are decreasing not so rapidly (Fig. 7) as in case of single 1D and 2D NLS equations. 6. Conclusions A reduced model is derived for the NLS equation by preserving the Hamiltonian structure. A priori error estimates are obtained for the mid-point rule as time integrator for the reduced dynamical system. Numerical results show that the energy
B. Karasözen et al. / Applied Mathematics and Computation 258 (2015) 509–519
519
and the phase space structure of the three different NLS equations are well preserved by using few POD modes. The number of the POD modes containing most of the energy depends on the decay of the singular values of the snapshot matrix, reflecting the dynamics of the underlying systems. In a future work, we will investigate the dependence of the ROM solutions on parameters for the CNLS equation by performing a sensitivity analysis. Acknowledgements The authors wish to thank the referees for their helpful suggestions and comments. This research was supported by the Middle East Technical University Scientific Research Fund (Project: BAP-07–05-2012–102). References [1] A. Aydn, B. Karasözen, Symplectic and multi-symplectic methods for coupled nonlinear Schrödinger equations with periodic solutions, Comput. Phys. Commun. 177 (2007) 66–583.