Mass- and energy-conserved numerical schemes for nonlinear Schr\"odinger equations

TL;DR

This paper introduces new mass- and energy-conserving numerical schemes for nonlinear Schrödinger equations, with analyses and efficient solvers, demonstrating their accuracy and ability to capture blow-up phenomena.

Contribution

It presents a family of novel time-stepping schemes that conserve mass and energy, along with error analysis and efficient solvers, for nonlinear Schrödinger equations.

Findings

Schemes conserve mass and energy accurately.

Numerical experiments confirm convergence.

Methods effectively capture blow-up phenomena.

Abstract

In this paper, we propose a family of time-stepping schemes for approximating general nonlinear Schr\"odinger equations. The proposed schemes all satisfy both mass conservation and energy conservation. Truncation and dispersion error analyses are provided for each proposed scheme. Efficient fixed-point iterative solvers are also constructed to solve the resulting nonlinear discrete problems. As a byproduct, an efficient one-step implementation of the BDF schemes is obtained as well. Extensive numerical experiments are presented to demonstrate the convergence and the capability of capturing the blow-up phenomenon of the proposed schemes.

Tables

Table 1. Table 1: Parameters β j subscript 𝛽 𝑗 \beta_{j} for some proposed schemes.

Schemes	${\beta }_0$	${\beta }_1$	${\beta }_2$	${\beta }_3$	${\beta }_4$	${\beta }_5$
Crank-Nicolson	1
Leapfrog	$\frac{1}{2}$	$\frac{1}{2}$
M-BDF2	$\frac{3}{2}$	$-\frac{1}{2}$
M-BDF3	$\frac{11}{6}$	$-\frac{7}{6}$	$\frac{1}{3}$
M-BDF4	$\frac{25}{12}$	$-\frac{23}{12}$	$\frac{13}{12}$	$-\frac{1}{4}$
M-BDF5	$\frac{137}{60}$	$-\frac{163}{60}$	$\frac{137}{60}$	$-\frac{21}{20}$	$\frac{1}{5}$
M-BDF6	$\frac{147}{60}$	$-\frac{213}{60}$	$\frac{237}{60}$	$-\frac{163}{60}$	$\frac{31}{30}$	$-\frac{1}{6}$
4-Step Symmetric	$-\frac{1}{12}$	$\frac{7}{12}$	$\frac{7}{12}$	$-\frac{1}{12}$

Table 2. Table 2: Parameters α s , j subscript 𝛼 𝑠 𝑗 \alpha_{s,j} and β s , j subscript 𝛽 𝑠 𝑗 \beta_{s,j} .

s	${\alpha }_{s,0}$	${\alpha }_{s,1}$	${\alpha }_{s,2}$	${\alpha }_{s,3}$	${\alpha }_{s,4}$	${\alpha }_{s,5}$	${\alpha }_{s,6}$	${\beta }_{s,0}$	${\beta }_{s,1}$	${\beta }_{s,2}$	${\beta }_{s,3}$	${\beta }_{s,4}$	${\beta }_{s,5}$
2	$\frac{3}{2}$	$-2$	$\frac{1}{2}$					$\frac{3}{2}$	$-\frac{1}{2}$
3	$\frac{11}{6}$	$-3$	$\frac{3}{2}$	$-\frac{1}{3}$				$\frac{11}{6}$	$-\frac{7}{6}$	$\frac{1}{3}$
4	$\frac{25}{12}$	$-4$	3	$-\frac{4}{3}$	$\frac{1}{4}$			$\frac{25}{12}$	$-\frac{23}{12}$	$\frac{13}{12}$	$-\frac{1}{4}$
5	$\frac{137}{60}$	$-5$	5	$-\frac{10}{3}$	$\frac{5}{4}$	$-\frac{1}{5}$		$\frac{137}{60}$	$-\frac{163}{60}$	$\frac{137}{60}$	$-\frac{21}{20}$	$\frac{1}{5}$
6	$\frac{147}{60}$	$-6$	$\frac{15}{2}$	$-\frac{20}{3}$	$\frac{15}{4}$	$-\frac{6}{5}$	$\frac{1}{6}$	$\frac{147}{60}$	$-\frac{213}{60}$	$\frac{237}{60}$	$-\frac{163}{60}$	$\frac{31}{30}$	$-\frac{1}{6}$

Table 3. Table 3: Dispersion error rates of Crank-Nickson scheme.

$\lambda$	$k$	$\tau$	Dispersion errors	Error rates
2	1	1E-01	0.003324	–
		1E-02	3.323244E-05	1.9989
		1E-03	3.333332E-07	2.0000
		1E-04	3.333333E-09	2.0000

Table 4. Table 4: Dispersion error rates of the Leapfrog scheme.

$\lambda$	$k$	$\tau$	Dispersion errors	Error rates
2	1	1E-01	0.021354	–
		1E-02	2.253333E-04	1.9767
		1E-03	2.266667E-06	1.9974
		1E-04	2.250000E-08	2.0032

Table 5. Table 5 : H ​ ( ω ~ , τ ) 𝐻 ~ 𝜔 𝜏 H(\tilde{\omega},\tau) of the numerical dispersion relation equations of the modified BDF schemes.

M–BDFs	$H(\tilde{\omega },\tau )$
M–BDF2	$\frac{-1}{2}\left[3\cos \left(\tilde{\omega }\tau \right)-5\right]$
M–BDF3	$\frac{-1}{36}\left[47\cos \left(2\tilde{\omega }\tau \right)-83\right]$
M–BDF4	$\frac{-1}{72}\left[913\cos \left(\tilde{\omega }\tau \right)-394\cos \left(2\tilde{\omega }\tau \right)+75\cos \left(3\tilde{\omega }\tau \right)-666\right]$
M–BDF5	$\frac{-1}{1800}\left[54049\cos \left(\tilde{\omega }\tau \right)-30682\cos \left(2\tilde{\omega }\tau \right)+10587\cos \left(3\tilde{\omega }\tau \right)-1644\cos \left(4\tilde{\omega }\tau \right)+34110\right]$
M–BDF6	$\frac{-1}{1800}[131149\cos \left(\tilde{\omega }\tau \right)-85882\cos \left(2\tilde{\omega }\tau \right)+39537\cos \left(3\tilde{\omega }\tau \right)-11244\cos \left(4\tilde{\omega }\tau \right)$
	$+1470\cos \left(5\tilde{\omega }\tau \right)-76830]$

Table 6. Table 6: Dispersion error rates of M-BDF schemes.

Modified BDFs	$\lambda$	$k$	$\tau$	Dispersion errors	Error rates
M–BDF2	2	1	1E-02	3.752127E-04	–
			1E-03	3.750001E-06	2.0002
			1E-04	3.750000E-08	2.0000
M–BDF3			1E-02	2.467750E-05	–
			1E-03	2.499677E-05	1.9944
			1E-04	2.499996E-07	2.0000
M–BDF4			1E-02	2.521704E-05	–
			1E-03	2.500217E-07	2.0037
			1E-04	2.500002E-09	2.0000
M–BDF5			1E-02	2.500172E-05	–
			1E-03	2.500001E-07	2.0000
			1E-04	2.500000E-09	2.0000
M–BDF6			1E-02	2.500124E-05	–
			1E-03	2.500001E-07	2.0000
			1E-04	2.500005E-09	2.0000

Table 7. Table 7: Accuracy test of the Crank-Nicolson scheme for NLS equation ( 102 ) with the exact solution ( 104 ). ( h = 40 / N ) ℎ 40 𝑁 (h=40/N) and t = 2 𝑡 2 t=2 .

Crank-Nicolson		Real Part				Imaginary Part
$N$	$\tau$	$L^2$ error	Order	$L^{\infty }$ error	Order	$L^2$ error	Order	$L^{\infty }$ error	Order
4000	1/8	0.5194	–	0.5136	–	0.5290	–	0.4351	–
4000	1/16	0.1337	1.96	0.1073	2.26	0.1299	2.03	0.1057	2.04
4000	1/32	0.0319	2.07	0.0250	2.10	0.0311	2.06	0.0249	2.09

Table 8. Table 8: Accuracy test of the Leapfrog scheme for NLS equation ( 102 ) with the exact solution ( 104 ). ( h = 40 / N ) ℎ 40 𝑁 (h=40/N) and t = 2 𝑡 2 t=2 .

Leapfrog			Real Part				Imaginary Part
Initial values	$N$	$\tau$	$L^2$ error	Order	$L^{\infty }$ error	Order	$L^2$ error	Order	$L^{\infty }$ error	Order
Exact	4000	1/8	0.5194	–	0.5136	–	0.5290	–	0.4351	–
solutions	4000	1/16	0.1337	1.96	0.1073	2.26	0.1299	2.03	0.1057	2.04
	4000	1/32	0.0319	2.07	0.0250	2.10	0.0311	2.06	0.0249	2.09
Numerical	4000	1/8	0.4606	–	0.3583	–	0.5208	–	0.5005	–
solution	4000	1/16	0.1163	1.99	0.0971	1.88	0.1165	2.16	0.0969	2.37
(C-N)	4000	1/32	0.1594	0.74	0.1228	0.75	0.1596	0.74	0.1192	0.71

Table 9. Table 9: Accuracy test of modified BDF2 scheme for NLS equation ( 102 ) with the exact solution ( 104 ). ( h = 40 / N ) ℎ 40 𝑁 (h=40/N) and t = 2 𝑡 2 t=2 .

M–BDFs				Real Part				Imaginary Part
	Initial values	$N$	$\tau$	$L^2$ error	Order	$L^{\infty }$ error	Order	$L^2$ error	Order	$L^{\infty }$ error	Order
M–BDF2	Exact	4000	1/16	0.2557	–	0.2318	–	0.2627	–	0.2536	–
	solutions		1/32	0.0624	2.04	0.0558	2.06	0.0640	2.04	0.0618	2.04
			1/64	0.0149	2.07	0.0133	2.07	0.0152	2.07	0.0148	2.07
	Numerical	4000	1/4	0.7278	–	0.5425	–	0.7598	–	0.6963	–
	solution		1/8	0.3144	1.21	0.2751	0.98	0.3132	1.28	0.2664	1.39
	(C-N)		1/16	0.0535	2.55	0.0345	3.00	0.0506	2.63	0.0397	2.75
M–BDF3	Exact	4000	1/16	0.3848	–	0.3833	–	0.3833	–	0.3910	–
	solutions		1/32	0.0799	2.27	0.0770	2.32	0.0770	2.32	0.0799	2.29
			1/64	0.0168	2.25	0.0153	2.33	0.0153	2.33	0.0161	2.31
	Numerical	4000	1/4	0.7491	–	0.5313	–	0.7680	–	0.6889	–
	solution		1/8	0.3456	1.12	0.2821	0.91	0.3470	1.15	0.2892	1.25
	(C-N)		1/16	0.0555	2.64	0.0373	2.92	0.0505	2.78	0.0407	2.83
M–BDF4	Exact	4000	1/16	0.5192	–	0.5376	–	0.5443	–	0.4896	–
	solutions		1/32	0.1073	2.27	0.1043	2.37	0.1204	2.18	0.0999	2.29
			1/64	0.0267	2.01	0.0234	2.16	0.0310	1.96	0.0251	1.99
	Numerical	4000	1/4	0.7736	–	0.5937	–	0.7849	–	0.6807	–
	solution		1/8	0.3651	1.08	0.3245	0.87	0.3637	1.12	0.2969	1.20
	(C-N)		1/16	0.0584	2.64	0.0398	3.03	0.0534	2.77	0.0401	2.89
M–BDF5	Exact	4000	1/16	0.6498	1.47	0.6805	1.54	0.6567	1.57	0.5517	1.89
	solutions		1/32	0.1357	2.26	0.1321	2.37	0.1511	2.12	0.1199	2.20
			1/64	0.0361	1.91	0.0313	2.08	0.0413	1.87	0.0336	1.83
	Numerical	4000	1/4	0.7491	–	0.5313	–	0.7680	–	0.6889	–
	solution		1/8	0.3456	1.12	0.2821	0.91	0.3470	1.15	0.2891	1.25
	(C-N)		1/16	0.0555	2.64	0.0373	2.92	0.0505	2.78	0.0407	2.83
M–BDF6	Exact	4000	1/16	0.7683	–	0.8078	–	0.7492	–	0.5869	–
	solutions		1/32	0.1630	2.24	0.1589	2.35	0.1794	2.06	0.1389	2.08
			1/64	0.0444	1.88	0.0384	2.05	0.0505	1.83	0.0409	1.76
	Numerical	4000	1/4	0.8018	–	0.6350	–	0.8089	–	0.6751	–
	solution		1/8	0.3808	1.07	0.3515	0.85	0.3748	1.11	0.3017	1.16
	(C-N)		1/16	0.0579	2.72	0.0404	3.12	0.0529	2.83	0.0399	2.92

Table 10. Table 10: Comparison of the blow-up times found by different methods (with N 1 = 2000 subscript 𝑁 1 2000 N_{1}=2000 and N 2 = 4000 subscript 𝑁 2 4000 N_{2}=4000 ).

$\tau$	Schemes	$t^{max}$		$t_1^R$		$t_2^R$		$u^{max}$ values
		$N_1$	$N_2$	$N_1$	$N_2$	$N_1$	$N_2$	$N_1$	$N_2$
0.02	C-N	0.62	0.64	0.62		0.66		24.6002	24.6537
	M-BDF2	0.60						10.6163	10.6164
	M-BDF3							10.3393	10.3550
	M-BDF4							10.5676	10.5936
	M-BDF5	0.68						14.3712	14.5215
	M-BDF6							21.1538	21.4550
0.01	C-N	0.65		0.65		0.67		35.2965	35.2392
	M-BDF2	0.64						15.3288	15.3288
	M-BDF3							14.2396	14.4062
	M-BDF4							14.5168	14.7713
	M-BDF5	0.68						20.7615	21.1336
	M-BDF6							29.9041	30.7954
0.005	C-N	0.67	0.675	0.67		0.68		48.2457	48.2666
	M-BDF2	0.665						20.4022	20.4023
	M-BDF3							20.6082	20.8447
	M-BDF4							21.5725	22.2511
	M-BDF5	0.685						29.6170	30.2689
	M-BDF6							43.8004	44.7018

Equations

$$i{u}_t$$

$$u(0)$$

$$\mathcal{M}(u)(t)$$

$$\mathcal{H}(u)(t)$$

$$i{d}_t{R}^{n+1}=$$

$$u^{n+1}=$$

$$d_t{R}^{n+1}=\frac{R^{n+1}-R^n}{\tau },R^{n+1/2}=\frac{R^{n+1}+R^n}{2}.$$

$$R^n=\sum _{j=0}^{k-1}{\beta }_j{u}^{n-j}.$$

$${\mathcal{M}}^n:=\Vert R^n{\Vert }_{L^2}^2,{\mathcal{H}}^n:=\frac{1}{2}\Vert \nabla R^n{\Vert }_{L^2}^2+\frac{\lambda }{4}\Vert R^n{\Vert }_{L^4}^4.$$

$$\displaystyle\frac{i}{2\tau}\int_{\mathcal{D}}\bigl{(}R^{n+1}-R^{n}\bigr{)}\bigl{(}\bar{R}^{n+1}+\bar{R}^{n}\bigr{)}\,dx$$

$$\displaystyle\qquad=\|\nabla R^{n+1/2}\|_{L^{2}}^{2}+\frac{\lambda}{2}\int_{\mathcal{D}}\bigl{(}|R^{n}|^{2}+|R^{n+1}|^{2}\bigr{)}\bigl{|}R^{n+1/2}\bigr{|}^{2}\,dx.$$

$$\displaystyle\frac{1}{2\tau}\Re\Bigl{[}\int_{\mathcal{D}}\bigl{(}R^{n+1}-R^{n}\bigr{)}\bigl{(}\bar{R}^{n+1}+\bar{R}^{n}\bigr{)}\,dx\Bigr{]}=0.$$

$${\mathcal{M}}^{n+1}={\mathcal{M}}^n,\forall n\ge 0.$$

$$i\Vert d_t{R}^{n+1}{\Vert }_{L^2}^2=$$

$$\displaystyle+\frac{\lambda}{4\tau}\int_{\mathcal{D}}\bigl{(}|R^{n}|^{2}+|R^{n+1}|^{2}\bigr{)}\bigl{(}R^{n+1}-R^{n}\bigr{)}\bigl{(}\bar{R}^{n+1}+\bar{R}^{n}\bigr{)}\,dx.$$

$${\mathcal{H}}^{n+1}={\mathcal{H}}^n,\forall n\ge 0.$$

$$\displaystyle\frac{i}{\tau}(u^{n+1}-u^{n})=-\frac{1}{2}\Delta\bigl{(}u^{n+1}+u^{n}\bigr{)}+\frac{\lambda}{4}\Bigl{(}\bigl{|}u^{n}\bigr{|}^{2}+\bigl{|}u^{n+1}\bigr{|}^{2}\Bigr{)}\bigl{(}u^{n+1}+u^{n}\bigr{)}.$$

$$T{E}^n=$$

$$\displaystyle-\frac{\lambda}{4}\Bigl{(}\bigl{|}u(t_{n})\bigr{|}^{2}+\bigl{|}u(t_{n+1})\bigr{|}^{2}\Bigr{)}\bigl{(}u(t_{n+1})+u(t_{n})\bigr{)}$$

$$=$$

$$i{u}_t(t)+\Delta u(t)-\lambda |u(t){|}^{2}u(t)=0.$$

$${\mathcal{A}}_1=$$

$$-\frac{i}{\tau }u(t_{n+1/2})+\frac{i}{2}{u}_t(t_{n+1/2})-\frac{i\tau }{8}{u}_{tt}(t_{n+1/2})+\frac{i{\tau }^2}{2^{3}3!}{u}^{(3)}(t_{n+1/2})+O({\tau }^3)$$

$$=$$

$${\mathcal{A}}_2=$$

$$+\frac{1}{2}\Delta u(t_{n+1/2})-\frac{\tau }{4}(\Delta u{)}_t(t_{n+1/2})+\frac{{\tau }^2}{16}(\Delta u{)}_{tt}(t_{n+1/2})+O({\tau }^3)$$

$$=$$

$$\displaystyle\bigl{|}u(t_{n})\bigr{|}^{2}=$$

$$=$$

$$\displaystyle\bigl{|}u(t_{n+1})\bigr{|}^{2}=$$

$$=$$

$${\mathcal{A}}_3=$$

$$=$$

$$T{E}^n=$$

$$\displaystyle\,\frac{i}{2\tau}(u^{n+1}-u^{n-1})=-\frac{1}{4}\bigl{(}\Delta u^{n+1}+2\Delta u^{n}+\Delta u^{n-1}\bigr{)}$$

$$\displaystyle\qquad+\frac{\lambda}{8}\Bigl{(}\Bigl{|}\frac{u^{n}+u^{n-1}}{2}\Bigr{|}^{2}+\Bigl{|}\frac{u^{n+1}+u^{n}}{2}\Bigr{|}^{2}\Bigr{)}\bigl{(}u^{n+1}+2u^{n}+u^{n-1}\bigr{)}.$$

$$T{E}^n=$$

$$\displaystyle-\frac{\lambda}{8}\Bigl{[}\Bigl{|}\frac{u(t_{n})+u(t_{n-1})}{2}\Bigr{|}^{2}+\Bigl{|}\frac{u(t_{n+1})+u(t_{n})}{2}\Bigr{|}^{2}\Bigr{]}\bigl{(}u(t_{n+1})+2u(t_{n})+u(t_{n-1})\bigr{)}$$

$$=:$$

$${\mathcal{B}}_1=$$

Full text

1 Introduction.

In this paper we consider the following nonlinear Schrödinger (NLS) equation:

[TABLE]

where $\mathcal{D}\subset\mathbf{R}^{d}\,(d=1,2,3)$ is a bounded domain, $\lambda=\pm 1$, $T>0$ and $i=\sqrt{-1}$ stands for the imaginary unit. $u=u(x,t):\mathcal{D}_{T}\to\mathbb{C}$ is a complex-valued function. $f:\mathbf{R}_{+}\to\mathbf{R}_{+}$ is a given real-valued function, which could be different in different applications, e.g., see [15, 16] and the references therein. The best known $f$ is $f(s)=s$, which leads to the well-known nonlinear Schrödinger equation with cubic nonlinearity. To close the system, we also need to specify a boundary condition for $u$. In this paper we consider both homogeneous Dirichlet and periodic boundary condition (see Sections 7 and 8), although most of our derivations and proofs are independent of the boundary condition.

The above Schrödinger equation describes many physical phenomena in optics, mechanics, and plasma physics. Mathematically, the NLS equation is a prototypical dispersive wave equation, its solutions exhibit some intriguing properties such as energy conservation, soliton wave, and possible blow-ups [5, 20]. In particular, the equation preserves both the mass and the Hamiltonian energy, that is, the following quantities are constants in time:

[TABLE]

Here the dependence of $u$ on $x$ variable is suppressed for notational brevity. The case with positive $\lambda$ is called defocusing and with negative $\lambda$ is called focusing which allows for bright soliton solutions as well as breather solutions.

Dispersion and nonlinearity can interact to produce permanent and localized wave forms in nonlinear dispersive wave equations such as the Korteweg-de Vries (KdV) equation [8, 14, 23] and the cubic Schrödinger equation [6, 24]. A distinct feature of these equations is the infinite many conservation laws (conserved integrals as invariants), allowing for soliton solutions which emerge from collision unchanged over time. The quality of the numerical approximation hence hinges on how well the conserved integrals can be preserved at the discrete level. Numerical methods without this property may result in substantial phase and shape errors after long time integration. Indeed for some wave equations the invariant preserving high order numerical methods have been shown more accurate than non-conservative methods after long-time numerical integration (see [4, 11]).

For the nonlinear Schrödinger equation considered in this paper, a natural question is whether it is possible to design numerical schemes which conserve the mass and energy simultaneously. A lot of effort has been made to preserve the mass by high order spatial discretization such as spectral methods [1, 2, 3], and discontinuous Galerkin methods [12]. A modified numerical energy may be preserved approximately by the corresponding spatial discretization (see [10]). However, since those methods are based on a time-splitting technique [3], except the Crank-Nickson scheme, they are only mass-conserved. The objective of this work is to provide an attempt to address the above question. Specifically, in this paper we develop and analyze a family of mass- and energy-conserved time-stepping schemes for approximating the cubic and general nonlinear Schrödinger equations. It should be noted that the energy conservation is also achieved for a modified energy, instead of the original energy.

The rest of this paper is organized as follows. In Section 2 we present a general framework involving two sequences of time-stepping schemes, which is shown to preserve both mass and energy for arbitrary time-step sizes, for the cubic nonlinear Schrödinger equation. In Section 3 we present many specific examples of mass- and (modified) energy-conserved time stepping schemes which fit into the general framework, and derive the truncation errors for four of these schemes. In Section 4 we present an efficient iterative algorithm to solve the resulting nonlinear equations. In Section 5 we extend the framework and examples to the general Schrödinger equations with arbitrary nonlinearity. In Section 6 we present a dispersion error analysis and derive the convergence rates for the dispersion errors. In Section 7 we present numerical experiments to validate the theoretical results and to gauge the performance of the proposed schemes, especially the sharpness of the convergence rates. In Section 8 we present additional numerical experiments to demonstrate the capability of the proposed numerical schemes for resolving the blow-up phenomenon. The paper is completed with some concluding remarks and comments given in Section 9.

2 Semi-discretization in time: a general framework.

In this section we propose a family of energy-conserved time-stepping schemes for approximating the cubic nonlinear Schrödinger equation.

Let $\tau>0$ and $t=t_{n}=n\tau$ for $n=0,1,2,\cdots,N$ be a uniform mesh for $[0,T]$. Let $k$ be a positive integer. We propose the following general $k$-step time-stepping scheme for problem (1)–(2): Seeking $\{R^{n},u^{n}\}$ for $n=k,k+1,\cdots,N$ such that

[TABLE]

where we use notation

[TABLE]

Note that from (6) we see that $R^{n}$ is a linear combination of $u^{n},u^{n-1},\cdots,u^{n-k+1}$ given as follows:

[TABLE]

As expected, choosing the parameters $\{\beta_{j}\}$ is a delicate issue. The guideline we use (see the details in the next section) is to choose them such that $d_{t}R^{n+1}$ has a specific order of accuracy for approximating $u_{t}(t_{n+1})$. In Table 1 we list several sets of parameters $\{\beta_{j}\}$ to be used in scheme (5)–(6), which results in various specific schemes for the nonlinear Schrödinger equation (again, see the details in the subsequent section).

We also remark that scheme (5)–(6) produces two sequences, namely $\{R^{n}\}$ and $\{u^{n}\}$. The first sequence can be regarded as the auxiliary quantities which are generated by solving nonlinear equation (5), while the second sequence are obtained as linear combinations of the first one.

To prove a key mass- and energy-conservation property of scheme (5)–(6), we define the following discrete mass and energy

[TABLE]

We start with establishing the following conservation results for the solution of scheme (5)–(6).

Theorem 1**.**

The solution to scheme (5) and (6) satisfies $\mathcal{M}^{n}=\mathcal{M}^{0}$ and $\mathcal{H}^{n}=\mathcal{H}^{0}$ for all $n\geq 1$.

Proof.

We multiply equation (5) by $\bar{R}^{n+1/2}$, integrate it over $\mathcal{D}$ and use the integration by parts to get,

[TABLE]

Taking the imaginary part of the resulting equation to get

[TABLE]

It follows from the identity \Re\bigl{[}(a-b)(\bar{a}+\bar{b})\bigr{]}=|a|^{2}-|b|^{2} for $a,b\in\mathbb{C}$ that

[TABLE]

To show the second conservation property, we multiply (5) by $d_{t}\bar{R}^{n+1}$, integrate the equation over $\mathcal{D}$ and use the integration by parts to get

[TABLE]

Taking the real part on the equation and applying identity \Re\bigl{[}(|a|^{2}+|b|^{2})(a-b)(\bar{a}+\bar{b})\bigr{]}=|a|^{4}-|b|^{4} for $a,b\in\mathbb{C}$ yield

[TABLE]

The proof is completed. ∎

3 Specific schemes and their truncation error analysis.

In this section we propose a number of specific schemes by defining $R^{n}$ in terms of $u^{n},u^{n-1},\cdots,u^{n-k+1}$. In other words, we shall specify the choice of parameters $\{\beta_{j}\}_{j=0}^{k-1}$ for each scheme.

3.1 A modified Crank-Nicolson scheme ($k=1$).

By setting $R^{n}=u^{n}$, the modified Crank-Nicolson scheme is defined as

[TABLE]

The local truncation error (LTE) of the modified Crank-Nicolson scheme (16) is defined by

[TABLE]

where $u(t)$ is the true solution of the nonlinear Schrödinger equation as follows:

[TABLE]

Here and below the dependence of $u$ on $x$ variable is suppressed for notational brevity.

The first lemma establishes the local truncation error for the modified Crank-Nicolson scheme.

Lemma 2**.**

The local truncation error of the modified Crank-Nicolson scheme (16) ($k=1$) is $O(\tau^{2})$.

Proof.

By using a Taylor series expansion of $u(t_{n+1})$ and $u(t_{n})$ about $u(t_{n+1/2})$, the accurate order of derivative term $\mathscr{A}_{1}$ and Laplace term $\mathscr{A}_{2}$ are $O(\tau^{2})$.

[TABLE]

Note that

[TABLE]

which leads to

[TABLE]

Recall the definition of the true solution $u(t)$ and using (19)-(21) in (17) gives

[TABLE]

Thus the local truncation error of above Crank-Nicolson scheme is $O(\tau^{2})$. The proof is completed. ∎

3.2 A Leapfrog scheme ($k=2$).

For $k=2$ we set $\beta_{0}=\beta_{1}=\frac{1}{2}$ in (6) so that $R^{n}=\frac{1}{2}(u^{n}+u^{n-1})$, which inserted into (5)) leads to the following Leapfrog scheme:

[TABLE]

The local truncation error of the leapfrog scheme (23) is defined by

[TABLE]

The next lemma establishes the local truncation error for the above Leapfrog scheme.

Lemma 3**.**

The local truncation error of the Leapfrog scheme ($k=2$) (23) is $O(\tau^{2})$.

Proof.

We formally apply the Taylor series expansions of $u(t_{n+1})$ and $u(t_{n-1})$ about $u(t_{n})$ for $\mathscr{B}_{1}$ to get

[TABLE]

Applying the Taylor series expansions of $u(t_{n+1})$ and $u(t_{n-1})$ about $u(t_{n})$ again for $\mathscr{B}_{2}$, we have

[TABLE]

By using the following facts for $\mathscr{B}_{3}$

[TABLE]

the order of nonlinear term can be estimated as

[TABLE]

Combing (25)-(29) together in (3.2) leads to

[TABLE]

Thus the local truncation error of above Leapfrog scheme is $O(\tau^{2})$. The proof is completed. ∎

3.3 Modified BDF schemes ($k=s$).

Let $s>1$. We recall that the s-step BDF scheme approximates the time derivative $u_{t}(t_{n+1})$ as follows:

[TABLE]

where $\alpha_{s,j}$ are given in Table 2. Our idea is to rewrite the above BDF expression as a first order backward difference, that is,

[TABLE]

which turns out is possible. Somehow this simple reformulation has not been seen in the literature before.

In order to determine the coefficients $\beta_{s,j}$, we solve them using the following identity:

[TABLE]

This holds true for all $u^{j}$ as long as the following matrix equation is satisfied:

[TABLE]

It is easy to check that

[TABLE]

Thus our modified BDF schemes are defined as

[TABLE]

where

[TABLE]

and $\beta_{s,j}$ are given in the following table (note that since BDF methods with $s>6$ are not zero-stable, so we only present s-step BDF with $s\leq 6$ here).

Set $\widetilde{R}(t_{n}):=\sum\limits_{j=0}^{s-1}\beta_{s,j}u(t_{n-j})$ as a linear combination of exact solution values, then the local truncation error of the modified BDF schemes (56) is defined by

[TABLE]

Lemma 4**.**

The local truncation error of the modified BDF schemes (56) ($k=s$) is $O(\tau^{2})$.

Proof.

Recall the approximation of the derivative yields $\mathscr{T}_{1}=iu_{t}(t_{n+1})+O(\tau^{s})$. Using the Taylor series expansion of $u(t_{n+1-j})$ about $u(t_{n})$ we have

[TABLE]

which implies

[TABLE]

where we have used the facts that

[TABLE]

For the nonlinear term, we have the following estimates:

[TABLE]

[TABLE]

where we use $\sum\limits_{j=0}^{s-1}\beta_{s,j}=\sum\limits_{j=1}^{s}\beta_{s,{j-1}}=1$, $\sum\limits_{j=1}^{s-1}j\beta_{s,j}\neq 0$ and $\sum\limits_{j=1}^{s}j\beta_{s,{j-1}}\neq 0$. Since $\sum\limits_{j=1}^{s-1}j\beta_{s,j}+\sum\limits_{j=1}^{s}j\beta_{s,{j-1}}=0$, hence,

[TABLE]

Similar to the estimates in (59), we have $\widetilde{R}(t_{n+1})+\widetilde{R}(t_{n})=2u(t_{n+1})+O(\tau^{2})$ and

[TABLE]

Combining (58) and (63) in (57), we obtain

[TABLE]

Thus the local truncation error is $O(\tau^{2})$. ∎

Remark 1**.**

As a by-product, the above construction also gives a (one-step) backward Euler reformulation for BDF schemes. Recall that the $s$-stage BDF scheme for $u^{\prime}(t)=f(t,u(t))$ is defined as

[TABLE]

Since

[TABLE]

then we can rewrite the BDF scheme as

[TABLE]

where

[TABLE]

Hence, each BDF scheme can be implemented as a one-step backward Euler scheme as (65) shows.

3.4 A four-step symmetric scheme ($k=4$).

To define this scheme, we set

[TABLE]

which fits (8) with $\beta=\frac{1}{12}(-1,7,7,-1)^{\top}$.

Lemma 5**.**

The local truncation error of the four-step symmetric scheme ($k=4$) is $O(\tau^{2})$.

Proof.

The derivative term is $O(\tau^{4})$ as follows

[TABLE]

For the Laplace term, we obtain

[TABLE]

For the nonlinear term, we have the following estimates:

[TABLE]

The proof is complete. ∎

4 An efficient fixed-point nonlinear solver.

To solve the nonlinear equation (5), we adapt the fixed-point iterative algorithm of [12] to the job. The proposed algorithm is defined below.

Algorithm 1

Step 1: Given $u^{l}$ for $l=0,1,...,n-k+1$, set

[TABLE]

Step 2: Update $R^{n+1}$ as follows: define $\{w^{l}\}_{l=0}^{L}$ iteratively by solving

[TABLE]

such that $\|w^{L}-w^{L-1}\|\leq\delta$ for some prescribed tolerance level $\delta$, then set

[TABLE]

Step 3: Update $u^{n+1}$ from $R^{n+1}$ by (6), that is,

[TABLE]

We note that any spatial discretization method, such as finite element, spectral and discontinuous Galerkin methods, can be employed in combination with the above algorithm to solve the nonlinear Schrödinegr equation (cf. [12]).

5 Extensions to Schrödinger equations with arbitrary nonlinearity.

We consider the following initial-boundary value problem for the general nonlinear Schrödinger equation :

[TABLE]

We extend scheme (5)–(6) as follows for problem (73)–(74): Seeking $\{R^{n},u^{n}\}$ for $n=k,k+1,\cdots,N$ such that

[TABLE]

where $G(a,b)$ is the following two variable function:

[TABLE]

For example, $G(a,b)=\sum_{j=0}^{p}a^{p-j}b^{j}$ if $f(s)=s^{p+1}$. Again, from (76) we have

[TABLE]

Define

[TABLE]

we have the following mass- and energy-conservation property of scheme (75)–(76).

Lemma 6**.**

The solution to scheme (75) and (76) satisfies $\mathcal{M}^{n}_{g}=\mathcal{M}^{0}_{g}$ and $\mathcal{H}^{n}_{g}=\mathcal{H}^{0}_{g}$ for all $n\geq 1$.

The nonlinear solver, Algorithm 1, now is replaced by the following modified algorithm.

Algorithm 2

Step 1: Given $u^{l}$ for $l=0,1,...,n-k+1$, set

[TABLE]

Step 2: Update $R^{n+1}$ as follows: define $\{w^{l}\}_{l=0}^{L}$ iteratively by solving

[TABLE]

such that $\|w^{L}-w^{L-1}\|\leq\delta$ for some prescribed tolerance level $\delta$, then set

[TABLE]

Step 3: Update $u^{n+1}$ from $R^{n+1}$ by (76), that is,

[TABLE]

Again, we remark that any spatial discretization method, such as finite element, spectral and discontinuous Galerkin methods, can be employed in combination with the above algorithm to solve the nonlinear Schrödinegr equation.

6 Dispersion error analysis.

In this section, we analyze the difference between the exact and numerical dispersion relations for the nonlinear Schrödinger equation and investigate ways to reduce the dispersive error generated by our mass- and energy-conserved time-stepping scheme (5)–(6). To minimize the numerical phase error while accurately solving the Schrödinger equation, the idea of preserving dispersion relation equation (DRE), which was proposed earlier in the area of computational aeroacoustics by Tam and Webb [19], is adopted. We refer the reader to [13] for a discussion of other structure-preserving algorithms for solving ordinary differential equations.

Consider the cubic nonlinear Schrödinger equation,

[TABLE]

Substituting the plane wave solution $u=\exp(i(kx-\omega t))$ into equation (82), the relation between the angular frequency $\omega$ and the wave number $k$ is given by [17]

[TABLE]

To derive the numerical dispersion relation equation for scheme (5)–(6), the discrete plane wave solution of the form $u^{n}=e^{i(kx-\tilde{\omega}n\tau)}$ is utilized, where $\tilde{\omega}$ is the numerical angular frequency.

Lemma 7**.**

The numerical dispersion relation of the Crank-Nickson scheme for the cubic nonlinear Schrödinger equation (82) is given by

[TABLE]

Proof.

Substituting $u^{n}=e^{i(kx-\tilde{\omega}n\tau)}$ into the discrete Crank-Nicolson scheme (16) we get

[TABLE]

Multiplying the both sides of (85) by $e^{i\tilde{\omega}\tau/2}$ to obtain

[TABLE]

By using the identities $e^{ix}-e^{-ix}=2i\sin x$ and $e^{ix}+e^{-ix}=2\cos x$, it follows from (86) that

[TABLE]

Hence, (84) holds. The proof is complete. ∎

To analyze the difference between the exact and numerical dispersions for nonlinear Schrödinger equations, we define the following dispersion error

[TABLE]

Table 3 shows the computed dispersion errors and the convergence order for the modified Crank-Nicolson scheme. The numerical results indicate that this scheme has a second order dispersion error.

Lemma 8**.**

The numerical dispersion relation of the Leapfrog scheme for the cubic nonlinear Schrödinger equation (82) is given by

[TABLE]

Proof.

Setting $u^{n}=e^{i(kx-\tilde{\omega}n\tau)}$ and using the identities $e^{ix}+e^{-ix}=2\cos x$ and $\cos 2\theta=2\cos^{2}\theta-1$ in \bigl{|}\frac{1}{2}\bigl{(}u^{n}+u^{n-1}\bigr{)}\bigr{|}^{2} and \bigl{|}\frac{1}{2}\bigl{(}u^{n}+u^{n+1}\bigr{)}\bigr{|}^{2} we get

[TABLE]

Substituting $u^{n}=e^{i(kx-\tilde{\omega}n\tau)}$ into the Leapfrog scheme (23) and using (89), we get the following equality:

[TABLE]

Using the identities $e^{ix}-e^{-ix}=2i\sin x$, $e^{ix}+e^{-ix}=2\cos x$, $\sin 2\theta=2\sin\theta\cos\theta$, and $\cos 2\theta=2\cos^{2}\theta-1$ in (90), we then obtain the following numerical dispersion relation equation of the Leapfrog scheme for the cubic nonlinear Schrödinger equation:

[TABLE]

Hence, (88) holds. The proof is complete. ∎

Table 4 shows the computed dispersion errors and convergence rates for the Leapfrog scheme. The numerical results indicate that the Leapfrog scheme also has a second order dispersion error.

Lemma 9**.**

The numerical dispersion relations of the modified BDF schemes for the cubic nonlinear Schrödinger equation (82) are given by

[TABLE]

where $H(\tilde{\omega},\tau)$ is given in Table 5 below.

Proof.

We only present the proofs of the numerical dispersion relation equations for the modified BDF2 and BDF3 schemes because the proofs for the remaining ones are similar.

We first consider the modified BDF2 scheme. Setting $u^{n}=e^{i(kx-\tilde{\omega}n\tau)}$ and using the identities $e^{ix}+e^{-ix}=2\cos x$ and $\cos 2\theta=2\cos^{2}\theta-1$ in \bigl{|}\frac{1}{2}\bigl{(}3u^{n}-u^{n-1}\bigr{)}\bigr{|}^{2}, we get

[TABLE]

A similar calculation applies to \bigl{|}\frac{1}{2}\bigl{(}3u^{n+1}-u^{n}\bigr{)}\bigr{|}^{2}. Substituting $u^{n}=e^{i(kx-\tilde{\omega}n\tau)}$ into the modified BDF2 scheme (56) and using (92) yield

[TABLE]

where H_{1}(\omega,\tau)=\frac{1}{2}\bigl{(}5-3\cos(\tilde{\omega}\tau)\bigr{)}.

It is easy to obtain the following equation from (92):

[TABLE]

Using the identities $e^{ix}-e^{-ix}=2i\sin x$ and $e^{ix}+e^{-ix}=2\cos x$ in (94), we get

[TABLE]

Since $3e^{-i\tilde{\omega}\tau/2}-e^{i\tilde{\omega}\tau/2}\neq 0$ in (95), the numerical dispersion relation of the modified BDF2 scheme can be written as

[TABLE]

The desired equation (91) holds by letting $H(\omega,\tau)=H_{1}(\omega,\tau)$.

Next, we consider the modified BDF3 scheme. Using the identities $e^{ix}+e^{-ix}=2\cos x$ and $\cos{2\theta}=2\cos^{2}{\theta}-1$ we get

[TABLE]

Substituting $u^{n}=e^{i(kx-\tilde{\omega}n\tau)}$ into the modified BDF3 scheme, we obtain

[TABLE]

where H_{2}(\omega,\tau)=\frac{1}{36}\bigl{(}83-47\cos(2\tilde{\omega}\tau)\bigr{)}.

Similar to the derivation of (95), using identities $e^{ix}-e^{-ix}=2i\sin x$, $e^{ix}+e^{-ix}=2\cos x$, $e^{ix}-e^{-ix}=2i\sin x$, and $e^{ix}+e^{-ix}=2\cos x$ in (98), we get

[TABLE]

where we have used the following identity in (100):

[TABLE]

Since $(11e^{-i\tilde{\omega}\tau/2}-9e^{i\tilde{\omega}\tau/2}+4\cos{(\tilde{\omega}\tau/2)}e^{i\tilde{\omega}\tau})\neq 0$ in (100), the numerical dispersion relation for the modified BDF3 can be written as

[TABLE]

which gives (91) after setting $H(\omega,\tau)=H_{2}(\omega,\tau)$. The proof is complete. ∎

Table 6 shows the computed dispersion errors and convergence rates for the modified BDF schemes. The numerical results indicate that these modified BDF schemes have a second order dispersion error.

7 Numerical experiments: validating the convergence rates.

In this section, we present several 1D numerical tests to illustrate our theoretical results, in particular, to verify the rates of convergence of the proposed time-stepping schemes. Our computations are done using the software package FEniCS and the linear finite element method is employed for the spatial discretization in all our numerical tests.

We consider the cubic nonlinear Schrödinger equation [18, 22, 2] (i.e., $f(s)=s$, $\lambda=-2$)

[TABLE]

where initial data $u_{0}={\rm sech}(x)\exp(2ix)$ is chosen so that the exact solution is given by [9]

[TABLE]

We solve problem (102)–(103) with the periodic boundary condition in $[-20,20]$ by a few selected schemes from the family of the time–stepping schemes proposed in Section 3. In Theorem 1, we proved the mass- and energy-conservation properties of scheme (5)–(6) without any restrictions on $h$ and $\tau$. To show the conservation properties of problem (102)–(103), uniform spatial and temporal meshes are used with $h=2^{-5}$ and $\tau=2^{-6}$ and the discrete mass and energy are defined in (9).

The evolution of the $L^{2}$-norm of $R^{n}$ and $u^{n}$ for one trajectory of the modified Crank-Nicolson scheme and Leapfrog scheme are shown in Figure 1. For the Crank-Nicolson scheme, we also show in Figure 1 that the evolution of the $L^{2}$-norm is exactly conserved although at each time step, the nonlinear equation is not exactly solved. Meanwhile, the energy is also exactly conserved as shown in Figure 1.

It should be noted that the Leapfrog scheme is a multi-step method, it requires two starting values, which are usually generated by a one-step method. In order to choose a suitable starting one-step scheme, we present computational results using different starting values generated by the exact solution and the Crank-Nicolson scheme. The computed evolution of $L^{2}$-norm and the energy of the Leapfrog scheme are shown in Figure 1 and Figure 1, from which we observe that the evolution has a large oscillation at the early stage, and gradually becomes stable.

The evolution of the $L^{2}$-norm for one trajectory of the modified BDF schemes are show in Figure 1,1 and Figure 2, 2,2. As expected, all of them gradually become stable and are exactly conserved. The differences of these oscillations (the zoom-in graphics at the beginning stage of the evolution) are obtained in Figure 2. We observe that all modified BDF schemes start with some oscillations and the amplitudes of the oscillations quickly diminish with time. In addition, similar phenomenon for the energy are also seen in Figure 1,1, Figure 2,2, 2 and Figure 2.

The accuracy of a proposed method is examined numerically by comparing the solution obtained on a sequence of coarse (time) meshes with the exact solution given in (104). The computed errors and rates of the Crank-Nicolson scheme are shown in Table 7. We observe that the $L^{2}$-norm errors decrease by a factor $2$ when the step-size $\tau$ is halved. Hence, a second order convergence rate is verified.

Recall that the Leapfrog scheme is a multi-step numerical method. In order to choose a suitable starting (one-step) scheme, we present some convergence results using difference starting values generated by the exact solution and by the Crank-Nickson scheme in Table 8. It is clear to see that the $L^{2}$-norm error rate is $2$, which confirms our theoretical result.

We run the same tests for the proposed modified BDF schemes, the computed results are shown in Table 9, we again observe that the $L^{2}$-norm error rate is $2$.

We conclude this section by presenting a convergence and performance comparison of the Leapfrog scheme and the modified BDF schemes in Figure 3.

8 Numerical experiments: capturing the blow-up time.

Our aim in this section is to present a numerical study of the blow-up phenomenon for the quintic nonlinear Schrödinger equation, which is known to be very delicate to simulate in order to have an accurate understanding of this behavior. In particular, we focus on predicating and capturing the blow-up time using several proposed schemes. Existing numerical results have shown that starting with an initial condition of a given amplitude, one can claim that the solution has a singularity as soon as its amplitude becomes three (or more) times bigger than the initial amplitude [7].

We consider the following quintic nonlinear Schrödinger problem:

[TABLE]

with periodic boundary conditions [12]. The initial condition is chosen as $u_{0}=1.6e^{-x^{2}}$. Since the initial energy is negative, it is known that the a blow-up in the solution must occur in finite time [7].

8.1 Comparison of $L^{\infty}$-norm profiles of the computed solutions by different schemes.

In this subsection, we want to test whether all or which of our proposed time-stepping schemes will be able to capture the blow-up phenomenon.

Figure 4 shows the simulation results of various schemes. We clearly see the formation of a singularity and that the used mesh size is small enough to capture the essential feature of the blow-up by all but the Leapfrog scheme in Figure 4, As a result, we conclude that the Leapfrog scheme is not, but all other proposed schemes are, capable of capturing the blow-up phenomenon of the quintic nonlinear Schrödinger equation (105).

The further comparison of $L^{\infty}$-norm profiles obtained by other schemes (excluding the Leapfrog scheme) are shown in Figure 4–4. The tests in Figure 4 indicate that the Crank-Nickson scheme and the modified BDF schemes are capable of capturing the blow-up phenomenon.

In order to study whether the blow-up phenomenon will affect the mass- and energy-conservation results of Theorem 1, we present the time evolution of the mass and energy of $R^{n}$ in Figure 5–5. As expected, the mass of $R^{n}$ is exactly conserved and the energy is also conserved before and after the blow-up time, in spite of a sharply increase in energy at the blow-up time. The similar behaviors of $u^{n}$ are observed from Figure 5–5.

8.2 Comparison of $L^{\infty}$–norm profiles obtained using different time-step sizes.

To better understand the sensitivity of the blow-up simulations to the time-step size, we analyze the capability of different schemes for capturing the blow-up phenomenon and provide three criteria for the blow-up time in this subsection.

The comparison of $L^{\infty}$–norm profiles obtained by different time-step sizes for the proposed time-stepping schemes are shown in Figure 6. We observe similar behavior for most schemes. Different simulation results for the Leapfrog scheme are shown in Figure 6 and for the linearized scheme of [21] in Figure 6, although small enough time-step sizes are used. As already mentioned earlier, our numerical tests show that the Crank-Nickson scheme (see Figure 6) and the modified BDF schemes (see Figure 6–6) are capable of capturing the blow-up phenomenon. These results also consist with our previous results.

In order to provide some criterion for identifying the blow-up time, we first propose three different criteria and then to present a comparison of them on different schemes in Table 10. The first criterion is to identify $t^{\max}$ corresponding to the time at which $\|u^{n}\|_{L^{\infty}}$ takes its maximum as shown in Figure 7, where $u^{\max}$ represents the maximum value. One difficulty with this criterion is that we may get different $t^{\max}$ with different schemes. To be specific, the modified lower order BDF schemes (i.e., M–BDF2, M–BDF3 and M–BDF4) identify the same earliest $t^{\max}$, while the modified higher order BDF schemes (i.e., M–BDF5 and M–BDF6) capture the same latest $t^{\max}$. In addition, the $t^{\max}$ found by the Crank-Nickson scheme is in the middle of the above two values, so it is inconclusive that which $t^{\max}$ is the most accurate.

However, it should be noted that the Crank-Nickson scheme is a preferable scheme for capturing the blow-up phenomenon because it finds the largest $\|u^{\max}\|_{L^{\infty}}$ as shown in Table 10.

The other two criteria identify $t_{1}^{n}$ and $t_{2}^{n}$ in Figure 7, where $t_{1}^{n}$ represents the time point at which the energy of $R^{n}$ is the smallest and $t_{2}^{n}$ denotes the time point at which the energy of $R^{n}$ has the maximum increase. From Table 10 we observe that all time-stepping schemes identify the same blow-up time using these two criteria, which shows the robustness of both criteria.

9 Conclusion.

In this paper we present a family of mass- and energy-conserved time-stepping schemes for general nonlinear Schrödinger equations. This includes the modified Crank-Nicolson scheme, Leapfrog scheme and modified BDF schemes as well as a four-step symmetric scheme. We have shown that the proposed schemes have second-order convergence while preserving both mass and energy in the discrete setting without any mesh constraint. We also derive the dispersion relation equation for each of the proposed schemes and numerically show the convergence orders for the numerical dispersions. Extensive numerical experiments have been presented to illustrate the performance of the proposed schemes and to validate the theoretical results of the paper. Additional numerical experiments have also be provided to test the capability of the proposed schemes for capturing the blow-up phenomenon of the quintic nonlinear Schrödinger equation. Various criteria are proposed for identifying the blow-up time and their effectiveness is also extensively examined. It is a bit disappointing that all proposed time-stepping schemes of this paper only have second order accuracy and second order truncation errors. A very interesting question is whether it is possible to improve these schemes into higher order schemes while still conserving both mass and (a modified) energy. Another challenging question is whether it is possible to construct mass- and energy-conserved linear schemes (that is, only a linear problem needs to be solved at each time step). These open questions are worthy of further investigation and will be addressed in a further work.

Acknowledgments. The work of the first author was partially supported by the NSF Grant DMS-1620168, and the work of the second author was partially supported by NSF Grant DMS-1812666. Part of the third author’s work was done during a recent visit to the University of Tennessee at Knoxville (UTK), the author would like to thank Department of Mathematics of UTK for the support and hospitality. The visit was financially supported by a scholarship from the author’s home institution, Northwestern Polytechnical University of China.