Stepsize-adaptive integrators for dissipative solitons in cubic-quintic complex Ginzburg-Landau equations

TL;DR

This survey evaluates adaptive exponential integrators for solving stiff cubic-quintic complex Ginzburg-Landau equations, focusing on accuracy near pulsating and exploding solitons with multiple time scales.

Contribution

It introduces and compares stepsize-adaptive exponential integrator schemes, including new formulas and a comoving frame approach for improved performance.

Findings

Adaptive integrators effectively handle multiple time scales.

The comoving frame improves phase rotation resolution.

Thorough numerical comparisons demonstrate method strengths.

Abstract

This paper is a survey on exponential integrators to solve cubic-quintic complex Ginzburg-Landau equations and related stiff problems. In particular, we are interested in accurate computation near the pulsating and exploding soliton solutions where different time scales exist. We explore stepsize-adaptive variations of three types of exponential integrators: integrating factor (IF) methods, exponential Runge-Kutta (ERK) methods and split-step (SS) methods, and their embedded versions for computation and comparison. We present the details, derive formulas for completeness, and consider seven different stepsize-adaptive integrating schemes to solve the cubic-quintic complex Ginzburg-Landau equation. Moreover, we propose using a comoving frame to resolve fast phase rotation for better performance. We present thorough comparisons and experiments in the one- and two-dimensional cubic-quintic complex Ginzburg-Landau equations.

Figures

Tables

Table 1. Table 1: Stiff order conditions for ERK up to order 3. Here, z = h ​ 𝐋 𝑧 ℎ 𝐋 z=h\mathbf{L} . J 𝐽 J denotes an arbitrary square matrix. ψ j , i ​ ( h ​ 𝐋 ) = ∑ k = 2 i − 1 a i ​ k ​ ( h ​ 𝐋 ) ​ c k j − 1 / ( j − 1 ) ! − c i j ​ φ j ​ ( c i ​ h ​ 𝐋 ) subscript 𝜓 𝑗 𝑖 ℎ 𝐋 superscript subscript 𝑘 2 𝑖 1 subscript 𝑎 𝑖 𝑘 ℎ 𝐋 superscript subscript 𝑐 𝑘 𝑗 1 𝑗 1 superscript subscript 𝑐 𝑖 𝑗 subscript 𝜑 𝑗 subscript 𝑐 𝑖 ℎ 𝐋 \psi_{j,i}(h\mathbf{L})=\sum_{k=2}^{i-1}a_{ik}(h\mathbf{L})c_{k}^{j-1}/(j-1)!-c_{i}^{j}\varphi_{j}(c_{i}h\mathbf{L}) . For the table up to order 5, see   [ 39 ] .

Order	Index	Stiff order condition
1	1	${\sum }_{i=1}^s{b}_i\left(z\right)={\phi }_1\left(z\right)$
2	2	${\sum }_{i=2}^s{b}_i\left(z\right)c_i={\phi }_2\left(z\right)$
	3	${\sum }_{j=1}^{i-1}{a}_{ij}\left(z\right)=c_i{\phi }_1\left(c_{i}z\right),i=2,\dots ,s$
3	4	${\sum }_{i=2}^s{b}_i\left(z\right)\frac{c_i^2}{2!}={\phi }_3\left(z\right)$
	5	${\sum }_{i=2}^s{b}_i\left(z\right)J{\psi }_{2,i}\left(z\right)=0$

Table 2. Table 2: Butcher table of IF4(3). Here, z = h ​ 𝐋 𝑧 ℎ 𝐋 z=h\mathbf{L} .

0
$\frac{1}{2}$	$\frac{1}{2}{e}^{z/2}$
$\frac{1}{2}$	0	$\frac{1}{2}$
1	0	0	$e^{z/2}$
1	$\frac{1}{6}{e}^z$	$\frac{1}{3}{e}^{z/2}$	$\frac{1}{3}{e}^{z/2}$	$\frac{1}{6}$
$b_i$	$\frac{1}{6}{e}^z$	$\frac{1}{3}{e}^{z/2}$	$\frac{1}{3}{e}^{z/2}$	$\frac{1}{6}$
${\overline{b}}_i$	$\frac{1}{6}{e}^z$	$\frac{1}{3}{e}^{z/2}$	$\frac{1}{3}{e}^{z/2}$	$\frac{1}{15}$	$\frac{1}{10}$

Table 3. Table 3: Butcher table of IF5(4). Here, z = h ​ 𝐋 𝑧 ℎ 𝐋 z=h\mathbf{L} .

$0$
$\frac{1}{5}$	$\frac{1}{5}{e}^{z/5}$
$\frac{3}{10}$	$\frac{3}{40}{e}^{3z/10}$	$\frac{9}{40}{e}^{z/10}$
$\frac{4}{5}$	$\frac{44}{45}{e}^{4z/5}$	$-\frac{56}{15}{e}^{3z/5}$	$\frac{32}{9}{e}^{z/2}$
$\frac{8}{9}$	$\frac{19372}{6561}{e}^{8z/9}$	$-\frac{25360}{2187}{e}^{31z/45}$	$\frac{64448}{6561}{e}^{53z/90}$	$-\frac{212}{729}{e}^{4z/45}$
$1$	$\frac{9017}{3168}{e}^z$	$-\frac{355}{33}{e}^{4z/5}$	$\frac{46732}{5247}{e}^{7z/10}$	$\frac{49}{176}{e}^{z/5}$	$-\frac{5103}{18656}{e}^{z/9}$
$1$	$\frac{35}{384}{e}^z$	$0$	$\frac{500}{1113}{e}^{7z/10}$	$\frac{125}{192}{e}^{z/5}$	$-\frac{2187}{6784}{e}^{z/9}$	$\frac{11}{84}$
$b_i$	$\frac{35}{384}{e}^z$	$0$	$\frac{500}{1113}{e}^{7z/10}$	$\frac{125}{192}{e}^{z/5}$	$-\frac{2187}{6784}{e}^{z/9}$	$\frac{11}{84}$
${\overline{b}}_i$	$\frac{5179}{57600}{e}^z$	$0$	$\frac{7571}{16695}{e}^{7z/10}$	$\frac{393}{640}{e}^{z/5}$	$-\frac{92097}{339200}{e}^{z/9}$	$\frac{187}{2100}$	$\frac{1}{40}$

Table 4. Table 4: Butcher table of ERK4(3)2(2). Here, φ i = φ i ​ ( h ​ 𝐋 ) subscript 𝜑 𝑖 subscript 𝜑 𝑖 ℎ 𝐋 \varphi_{i}=\varphi_{i}(h\mathbf{L}) , φ i , j = φ i ​ ( c j ​ h ​ 𝐋 ) subscript 𝜑 𝑖 𝑗 subscript 𝜑 𝑖 subscript 𝑐 𝑗 ℎ 𝐋 \varphi_{i,j}=\varphi_{i}(c_{j}h\mathbf{L}) .

0
$\frac{1}{2}$	$\frac{1}{2}{\phi }_{1,2}$
$\frac{1}{2}$	0	$\frac{1}{2}{\phi }_{1,3}$
1	$\frac{1}{2}{\phi }_{1,3}\left({\phi }_{0,3}-1\right)$	0	${\phi }_{1,3}$
1	${\phi }_1-3{\phi }_2+4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$4{\phi }_3-{\phi }_2$
$b_i$	${\phi }_1-3{\phi }_2+4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$4{\phi }_3-{\phi }_2$	0
${\overline{b}}_i$	${\phi }_1-3{\phi }_2+4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	0	$4{\phi }_3-{\phi }_2$

Table 5. Table 5: Butcher table of ERK4(3)3(3).

0
$\frac{1}{2}$	$\frac{1}{2}{\phi }_{1,2}$
$\frac{1}{2}$	$\frac{1}{2}{\phi }_{1,3}-{\phi }_{2,3}$	${\phi }_{2,3}$
1	${\phi }_{1,4}-2{\phi }_{2,4}$	0	$2{\phi }_{2,4}$
1	${\phi }_1-3{\phi }_2+4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$4{\phi }_3-{\phi }_2$
$b_i$	${\phi }_1-3{\phi }_2+4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$4{\phi }_3-{\phi }_2$	0
${\overline{b}}_i$	${\phi }_1-3{\phi }_2+4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	0	$4{\phi }_3-{\phi }_2$

Table 6. Table 6: Butcher table of ERK4(3)4(3).

$0$
$\frac{1}{2}$	$\frac{1}{2}{\phi }_{1,2}$
$\frac{1}{2}$	$\frac{1}{2}{\phi }_{1,3}-{\phi }_{2,3}$	${\phi }_{2,3}$
$1$	${\phi }_{1,4}-2{\phi }_{2,4}$	${\phi }_{2,4}$	${\phi }_{2,4}$
$\frac{1}{2}$	$\frac{1}{2}{\phi }_{1,5}-2a_{5,2}-a_{5,4}$	$a_{5,2}$	$a_{5,2}$	$\frac{1}{4}{\phi }_{2,5}-a_{5,2}$
$b_i$	${\phi }_1-3{\phi }_2+4{\phi }_3$	$0$	$0$	$-{\phi }_2+4{\phi }_3$	$4{\phi }_2-8{\phi }_3$
${\overline{b}}_i$	${\phi }_1-3{\phi }_2+4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$2{\phi }_2-4{\phi }_3$	$-{\phi }_2+4{\phi }_3$	$0$

Table 7. Table 7: Butcher table of ERK5(4)5(4).

$0$
$\frac{1}{2}$	$\frac{1}{2}{\phi }_{1,2}$
$\frac{1}{2}$	$\frac{1}{2}{\phi }_{1,3}-\frac{1}{2}{\phi }_{2,3}$	$\frac{1}{2}{\phi }_{2,3}$
$\frac{1}{4}$	$\frac{1}{4}{\phi }_{1,4}-\frac{1}{8}{\phi }_{2,4}$	$0$	$\frac{1}{8}{\phi }_{2,4}$
$\frac{1}{2}$	$\frac{1}{2}{\phi }_{1,5}-\frac{3}{2}{\phi }_{2,5}+2{\phi }_{3,5}$	$0$	$-\frac{1}{2}{\phi }_{2,5}+2{\phi }_{3,5}$	$2{\phi }_{2,5}-4{\phi }_{3,5}$
$\frac{1}{5}$	$\frac{1}{5}{\phi }_{1,6}-\frac{2}{25}{\phi }_{2,6}-\frac{1}{2}{a}_{64}$	$0$	$0$	$a_{6,4}$	$\frac{2}{25}{\phi }_{2,6}-\frac{1}{2}{a}_{6,4}$
$\frac{2}{3}$	$\frac{2}{3}{\phi }_{1,7}+\frac{125}{162}{a}_{64}-a_{75}-a_{76}$	$0$	$0$	$-\frac{125}{162}{a}_{64}$	$a_{75}$	$a_{76}$
$1$	${\phi }_{1,8}-a_{85}-a_{86}-a_{87}$	$0$	$0$	$0$	$a_{85}$	$a_{86}$	$a_{87}$
$1$	${\phi }_1-b_6-b_7-b_8$	$0$	$0$	$0$	$0$	$b_6$	$b_7$	$b_8$
$b_i$	${\phi }_1-b_6-b_7-b_8$	$0$	$0$	$0$	$0$	$b_6$	$b_7$	$b_8$	$0$
${\overline{b}}_i$	${\phi }_1-b_6-b_7-b_8$	$0$	$0$	$0$	$0$	$b_6$	$b_7$	$0$	$b_8$

Table 8. Table 8: Characteristics of the seven embedded schemes. ‘order’: the order of the local error of the embedded scheme. ‘ nn ’: the number of evaluations of 𝐍 ​ ( t , y ) 𝐍 𝑡 𝑦 \mathbf{N}(t,y) in a single step. ‘ nab ’: the number of distinct a i ​ j subscript 𝑎 𝑖 𝑗 a_{ij} and b i subscript 𝑏 𝑖 b_{i} entries in a Butcher table. ‘ E n + 1 subscript 𝐸 𝑛 1 E_{n+1} ’: the expression for the estimated local error.

scheme	order	nn	nab	$E_{n+1}$
IF4(3)	4	5	-	$\frac{1}{10}h\left[𝐍(t_n+h,Y_5)-𝐍(t_n+h,Y_4)\right]$
IF5(4)	5	7	-	$h\{\frac{-71}{57600}𝐍(t_n,Y_1)+\frac{71}{16695}𝐍(t_n+\frac{3h}{10},Y_3)$
				$-\frac{71}{1920}𝐍(t_n+\frac{4h}{5},Y_4)+\frac{17253}{339200}𝐍(t_n+\frac{8h}{9},Y_4)$
				$-\frac{22}{525}𝐍(t_n+h,Y_5)+\frac{1}{40}𝐍(t_n+h,Y_6)\}$
ERK4(3)2(2)	4	5	4	$b_{4}h\left[𝐍(t_n+h,Y_5)-𝐍(t_n+h,Y_4)\right]$
ERK4(3)3(3)	4	5	8	$b_{4}h\left[𝐍(t_n+h,Y_5)-𝐍(t_n+h,Y_4)\right]$
ERK4(3)4(3)	4	5	11	$b_{5}h\left[𝐍(t_n+h,Y_5)-\frac{1}{2}\left(𝐍(t_n+h/2,Y_2)+𝐍(t_n+h/2,Y_3)\right)\right]$
ERK5(4)5(4)	5	9	23	$b_{8}h\left[𝐍(t_n+h,Y_9)-𝐍(t_n+h,Y_8)\right]$
SS4(3)	4	24	-	$\left(y_1-y_2\right)/2$

Table 9. Table 9: Nonstiff 5th-order conditions.

	Tree	condition
1		$\sum {\beta }^{r,0}{\alpha }_r^{j,0}{c}_j^3=\frac{1}{20}$
2		$\sum {\beta }^{r,1}{c}_r^3=\frac{1}{20}$
3		$\sum {\beta }^{r,0}{\alpha }_r^{j,0}{\alpha }_j^{k,0}{c}_j{c}_k=\frac{1}{40}$
4		$\sum {\beta }^{r,1}{\alpha }_r^{j,0}{c}_r{c}_j=\frac{1}{40}$
5		$\sum {\beta }^{r,0}{\alpha }_r^{j,0}{\alpha }_j^{k,1}{c}_j=\frac{1}{40}$
6		$\sum {\beta }^{r,1}{\alpha }_r^{j,1}{c}_r=\frac{1}{40}$
7		$\sum {\beta }^{r,0}{\alpha }_r^{j,0}{\alpha }_j^{k,0}{c}_k^2=\frac{1}{60}$
8		$\sum {\beta }^{r,1}{\alpha }_r^{j,0}{c}_j^2=\frac{1}{60}$
9		$\sum {\beta }^{r,0}{\alpha }_r^{j,1}{c}_j^2=\frac{1}{60}$
10		$\sum {\beta }^{r,2}{c}_r^2=\frac{1}{60}$
11		$\sum {\beta }^{r,0}{\alpha }_r^{j,0}{\alpha }_j^{k,0}{\alpha }_k^{p,0}{c}_p=\frac{1}{120}$
12		$\sum {\beta }^{r,1}{\alpha }_r^{j,0}{\alpha }_j^{k,0}{c}_k=\frac{1}{120}$
13		$\sum {\beta }^{r,0}{\alpha }_r^{j,0}{\alpha }_j^{k,0}{\alpha }_k^{p,1}=\frac{1}{120}$
14		$\sum {\beta }^{r,1}{\alpha }_r^{j,0}{\alpha }_j^{k,1}=\frac{1}{120}$

Equations

$$A_t=\mu A+(D_r+i{D}_i)\Delta A+({\beta }_r+i{\beta }_i)|A{|}^{2}A+({\gamma }_r+i{\gamma }_i)|A{|}^{4}A$$

$${\dot{a}}_k=[\mu -(D_r+i{D}_i)q^2{k}^2]a_k+({\beta }_r+i{\beta }_i){\mathcal{F}}_k(|A{|}^{2}A)+({\gamma }_r+i{\gamma }_i){\mathcal{F}}_k(|A{|}^{4}A)$$

$${\dot{a}}_{mn}=[\mu -(D_r+i{D}_i)(q_x^2{m}^2+q_y^2{n}^2)]a_{mn}+({\beta }_r+i{\beta }_i){\mathcal{F}}_{mn}(|A{|}^{2}A)+({\gamma }_r+i{\gamma }_i){\mathcal{F}}_{mn}(|A{|}^{4}A)$$

$$L=L_x=L_y=50$$

$$y^{\prime }(t)=f(t,y)=\mathbf{L}y+\mathbf{N}(t,y),$$

$$z^{\prime }(t)=e^{-t\mathbf{L}}\mathbf{N}(t,e^{t\mathbf{L}}z).$$

$$e^{-t\mathbf{L}}\cdot {\frac{\partial N(t,y)}{\partial y}|}_{y=e^{t\mathbf{L}}z}\cdot e^{t\mathbf{L}}$$

$$Y_i=e^{h{c}_i\mathbf{L}}{y}_n+h\sum _{j=1}^{i-1}{a}_{ij}{e}^{h{\alpha }_{ij}\mathbf{L}}\mathbf{N}(t_n+c_{j}h,Y_j),y_{n+1}=e^{h\mathbf{L}}{y}_n+h\sum _{i=1}^s{b}_i{e}^{h{\beta }_i\mathbf{L}}\mathbf{N}(t_n+c_{i}h,Y_i),$$

$$y(t_{n+1})=e^{h\mathbf{L}}y(t_n)+h{\int }_0^1{e}^{(1-\theta )h\mathbf{L}}\mathbf{N}(t_n+\theta h,y(t_n+\theta h))d\theta .$$

$$\mathbf{N}(t_n+\theta h,y(t_n+\theta h))=\sum _{i=1}^s\prod _{j\ne i}\frac{\theta -c_j}{c_i-c_j}\mathbf{N}(t_n+c_{i}h,y(t_n+c_{i}h)).$$

$$y(t_{n+1})=e^{h\mathbf{L}}y(t_n)+h\sum _{i=1}^s{b}_i(h\mathbf{L})\mathbf{N}(t_n+c_{i}h,y(t_n+c_{i}h)),b_i(h\mathbf{L})={\int }_0^1{e}^{(1-\theta )h\mathbf{L}}\prod _{j\ne i}\frac{\theta -c_j}{c_i-c_j}d\theta .$$

$$Y_i=e^{h{c}_i\mathbf{L}}{y}_n+h\sum _{j=1}^{i-1}{a}_{ij}(h\mathbf{L})\mathbf{N}(t_n+c_{j}h,Y_j),y_{n+1}=e^{h\mathbf{L}}{y}_n+h\sum _{i=1}^s{b}_i(h\mathbf{L})\mathbf{N}(t_n+c_{i}h,Y_i).$$

$${\phi }_j(z)={\int }_0^1{e}^{(1-\theta )z}\frac{{\theta }^{j-1}}{(j-1)!}d\theta ,j\ge 1.$$

$$y(t_{n+1})=e^{h\mathbf{L}}y(t_n)+h\sum _{r=0}^{\infty }{h}^r{\mathbf{N}}^{(r)}{\int }_0^1\frac{e^{(1-\theta )h\mathbf{L}}{\theta }^r}{r!}d\theta =e^{h\mathbf{L}}y(t_n)+h\sum _{r=0}^{\infty }{h}^r{\mathbf{N}}^{(r)}{\phi }_{r+1}(h\mathbf{L}).$$

$${\phi }_{j+1}(z)=\frac{{\phi }_j(z)-1/j!}{z},{\phi }_j(0)=\underset{z\to 0}{\lim }{\phi }_j(z)=\frac{1}{j!},j\ge 0$$

$${\phi }_1(z)=\frac{e^z-1}{z},{\phi }_2(z)=\frac{e^z-z-1}{z^2},{\phi }_3(z)=\frac{e^z-z^2/2-z-1}{z^3}.$$

$$y_{n+1}={\varphi }_{\mathbf{N}}(b_{s}h)\circ e^{a_{s}h\mathbf{L}}\circ \cdots \circ {\varphi }_{\mathbf{N}}(b_{1}h)\circ e^{a_{1}h\mathbf{L}}\circ y_n$$

$$y^{\prime }(t)=\mathbf{N}(t,y).$$

$$E_{n+1}={\bar{y}}_{n+1}-y_{n+1}=h\sum _{i=1}^{s+1}({\bar{b}}_i-b_i)e^{h{\beta }_i\mathbf{L}}\mathbf{N}(t_n+c_{j}h,Y_j).$$

$${\bar{b}}_2=x{b}_5,{\bar{b}}_3=y{b}_5,{\bar{b}}_5=0$$

$$a_{5,2}=\frac{1}{2}{\phi }_{2,5}-{\phi }_{3,4}+\frac{1}{4}{\phi }_{2,4}-\frac{1}{2}{\phi }_{3,5}$$

$$a_{64}=\frac{8}{25}{\phi }_{2,6}-\frac{32}{125}{\phi }_{3,6},a_{75}=\frac{125}{1944}{a}_{64}-\frac{16}{27}{\phi }_{2,7}+\frac{320}{81}{\phi }_{3,7},a_{76}=\frac{3125}{3888}{a}_{64}+\frac{100}{27}{\phi }_{2,7}-\frac{800}{81}{\phi }_{3,7},$$

$$\varphi =\frac{5}{32}{a}_{64}-\frac{1}{28}{\phi }_{2,6}+\frac{36}{175}{\phi }_{2,7}-\frac{48}{25}{\phi }_{3,7}+\frac{6}{175}{\phi }_{4,6}+\frac{192}{35}{\phi }_{4,7}+6{\phi }_{4,8},a_{85}=\frac{208}{3}{\phi }_{3,8}-\frac{16}{3}{\phi }_{2,8}-40\varphi ,$$

$$a_{86}=-\frac{250}{3}{\phi }_{3,8}+\frac{250}{21}{\phi }_{2,8}+\frac{250}{7}\varphi ,a_{87}=-27{\phi }_{3,8}+\frac{27}{14}{\phi }_{2,8}+\frac{135}{7}\varphi ,b_6=\frac{125}{14}{\phi }_2-\frac{625}{14}{\phi }_3+\frac{1125}{14}{\phi }_4,$$

$$b_7=-\frac{27}{14}{\phi }_2+\frac{162}{7}{\phi }_3-\frac{405}{7}{\phi }_4,b_8=\frac{1}{2}{\phi }_2-\frac{13}{2}{\phi }_3+\frac{45}{2}{\phi }_4$$

$$(a_1,b_1,\dots ,a_s,b_s)=(b_s,a_s,\dots ,b_1,a_1),$$

$$y_1$$

$$y_2$$

$$y_{n+1}=\frac{1}{2}(y_1+y_2),{\bar{y}}_{n+1}=\frac{1}{2}(y_1-y_2).$$

$$(a_1,a_2,a_3)=(0.268330095781759925,-0.187991618799159782,0.919661523017399857).$$

$$||E_{n+1}|{|}_{\infty }<\text{rtol}\cdot ||y_{n+1}|{|}_{\infty }$$

$$h_{attempt}=s\cdot h,s=v{\left(\frac{\text{rtol}\cdot ||y_{n+1}|{|}_{\infty }}{||E_{n+1}|{|}_{\infty }}\right)}^{1/p}.$$

Full text

††journal: Computer Physics Communications

1 Introduction

The complex Ginzburg-Landau equation [1, 2] is one of the most frequently studied nonlinear equations in physics and applied mathematics community. Derived as a general amplitude equation near the onset of instability in the context of pattern formation [1], it has applications in various fields of physics ranging from nonlinear optics [3] to superconductivity [4]. In order to preserve the invariance of the equation under a phase change $A\to Ae^{i\phi}$, complex Ginzburg-Landau equations only have odd-order nonlinear terms. The form with a cubic term is frequently used to study spatiotemporal chaos and intermittent traveling waves. Recently, complex Ginzburg-Landau equation with an additional quintic term

[TABLE]

has attracted attention for its peculiar dissipative soliton solutions in one- and two-dimensional cases [5, 3, 6, 7, 8, 9]. Here, $A(x,t)$ or $A(x,y,t)$ is a complex field. $\Delta$ stands for the Laplace operator $\partial_{xx}$ or $\partial_{xx}+\partial_{yy}$. A dissipative soliton is a self-localized structure that arises in spatially extended dissipative systems. It maintains its shape temporally similar to a constantly propagating solitary wave packet in an optical medium. Dissipative solitons appear as a result of the balance between dispersion and nonlinearity and the balance between the gain and loss of energy. Dispersion spreads the shape while nonlinearity focuses it. A nontrivial internal energy flow and the dependence on an external energy source differentiate dissipative solitons from solitons in integrable systems. Apart from a phase shift, solitons in integrable systems remain unchanged during soliton-soliton interactions; while dissipative solitons are free from energy or momentum conservation during scattering and annihilation. For a more comprehensive discussion on dissipative solitons, see ref. [10]. While dissipative solitons lack most of the very special properties of soliton solutions of integrable Hamiltonian systems, they are generic and physically important, both because such structures are observed for wide ranges of equation parameters [11, 12], and because many types of (pulsating) soliton solutions are observed in laser experiments [13, 14, 15].

By exploring the parameter space of (1), Akhmediev et al.  [3, 6] identified several types of dissipative soliton solutions, such as plain soliton, pulsating soliton, zig-zag soliton, creeping soliton, composite soliton, exploding soliton, chaotic soliton and so on, some of which are shown in figure 1. These discoveries stimulate further investigation in soliton solutions in both one- and two-dimensional cubic-quintic complex Ginzburg-Landau equation, and the study of their bifurcations. For example, Soto-Crespo, Akhmediev and Chiang [7] found two coexisting solitons with a high and a low energy respectively. Chang, Soto-Crespo, Vouzas and Akhmediev [16, 17] found a new class of pulsating solitons with large ratios of maximal to minimal energies as shown in figure 1(b). The energy may change more than two orders of magnitude in each period. Tsoy and Akhmediev [18] studied bifurcations from stationary to pulsating solitons based on reducing the infinite-dimensional dynamics (1) to a five-dimensional model. Meanwhile, Mancas and Choudhury [19] obtained a three-dimensional model of (1) by a variational method in the study of pulsating and snake solitons. Among all the discoveries, solitons which undergo intermittent explosions stand out. As shown in figure 1(c), an exploding soliton moves uniformly for the most time with occasional substantial changes, after which it quickly restores back to the pre-explosion profile. Both symmetric and asymmetric explosions [20] are recorded and the center of the soliton shifts after asymmetric explosions. For the two-dimensional cubic-quintic complex Ginzburg-Landau equation, the center of an asymmetric exploding soliton exhibits a random walk with an anomalous diffusion rate  [8, 9], which is unexpected for a deterministic system.

All these novel phenomena expand our knowledge about dissipative solitons in spatially extended chaotic systems, and also propose a numerical challenge, i.e., how to efficiently integrate dissipative solitons. There is little concern about solitons that propagate uniformly such as plain or composite solitons. But for pulsating or exploding solitons whose profiles change periodically or intermittently with interleaved fast and slow dynamics, the fast-changing parts such as the high-amplitude parts in figure 1(b) and the exploding parts in figure 1(c) require smaller integration time steps. Thus stepsize-adaptive integration schemes are preferred over constant time-stepping schemes. The purpose of this paper is to introduce several stepsize-adaptive schemes to integrate dissipative solitons in cubic-quintic complex Ginzburg-Landau equation.

For complex Ginzburg-Landau equations, finite difference method [21], Fourier pseudo-spectral method [22] and Chebyshev-Tau spectral method [23] are applied as spatial discretization methods. In this paper, we explore the Fourier pseudo-spectral approach, for its spectral accuracy when the geometry of the solution is smooth and regular. Periodic boundary condition $x\in[0,L]$ or $(x,y)\in[0,L_{x}]\times[0,L_{y}]$ is enforced for the one- and two-dimensional cases respectively. In the Fourier space, equation (1) takes form of

[TABLE]

for the one-dimensional case, and

[TABLE]

for the two-dimensional case. Here, $q=2\pi/L$, $q_{x}=2\pi/L_{x}$ and $q_{y}=2\pi/L_{y}$. $a_{k}$ and $a_{mn}$ are the $k$th Fourier mode of $A(x,t)$ and $mn$th Fourier mode of $A(x,y,t)$ respectively. $\mathcal{F}_{k}(\cdot)$ and $\mathcal{F}_{mn}(\cdot)$ denote the $k$th and $mn$th component of the one- and two-dimensional discrete Fourier transform. The linear part is represented by a diagonal matrix in the frequency domain. The nonlinear part is first transformed back to the physical domain by the inverse fast Fourier transform (FFT), integrated in the physical domain and then transformed again to the frequency domain. This is why this strategy takes name “pseudo-spectral”. In soliton simulations, the size of the domain is set to

[TABLE]

which is large enough to hold a single soliton in the domain. $1024$ and $1024^{2}$ Fourier modes are used respectively in (2) and (3). The number of Fourier modes has been doubled in numerical experiments but without noticeable improvement in accuracy.

For the time derivative, the linear parts of (2) and (3) have quadratic structures and thus are stiff. Also, to model pulsating and exploding soliton solution accurately, popular single-step integrators such as Runge-Kutta methods requires an extremely small step size. Therefore, we explore exponential integrators to integrate the linear part explicitly, and use stepsize-adaptation for exponential integrators.

In this paper, we investigate the performance of three different types of stepsize-adaptive exponential integrators. Our motivation is to tackle cubic-quintic complex Ginzburg-Landau equation, and to the best of our knowledge, this is the first work to apply these methods to dissipative soliton integration in cubic-quintic complex Ginzburg-Landau equations. The main contributions of this paper are as follows. First, We not only convert two popular integrating factor methods and four popular exponential Rung-Kutta methods to their stepsize-adaptive versions, but also consider one stepsize-adaptive split-step method with symmetric coefficients. We note that the Runge-Kutta methods in the interaction picture used in the quantum mechanics community are equivalent to the integrating factor methods that will be explained in sect. 2.1. Second, we formulate an embedded lower-order method in a 4th-order exponential Runge-Kutta method in table 6 without an additional internal stage and are the first to study the embedded 5th-order exponential Runge-Kutta method which will be introduced in sect. 3. Third, we utilize invariant dynamical structures in the system to accelerate the integration process. We compute traveling waves in cubic-quintic complex Ginzburg-Landau equations and integrate the system in a comoving frame with respect to the traveling waves. This change of frame promotes the performance of exponential Runge-Kutta methods, which will be covered in sect. 6. Finally, we experiment and compare these methods, to numerically integrate dissipative solitons in cubic-quintic complex Ginzburg-Landau equations effectively.

This paper is organized as follows. In sect. 2, we review three different types of exponential integrators: integrating factor (IF) methods, exponential Rung-Kutta (ERK) methods and split-step (SS) methods. In sect. 3, we discuss how to embed a lower-order scheme in an exponential integrator and introduce seven representative schemes. The strategy to update step size and the metrics to evaluate the performance of an integrator is discussed in sect. 4. Sect. 5 is devoted to the numerical experiments on the one-dimensional cubic-quintic complex Ginzburg-Landau equation. We compare the performance of constant time-stepping schemes with that of stepsize-adaptive schemes. In sect. 6, we introduce the idea of using a comoving frame to alleviate the fast-phase rotation problem. Sect. 7 is devoted to the discussion of the performance of stepsize-adaptive schemes in the two-dimensional cubic-quintic complex Ginzburg-Landau equation. We summarize our discoveries in sect. 8.

2 Review of exponential integrators

Exponential integrators are formulated to solve a system of ODEs of semilinear type

[TABLE]

which is usually spatially discretized from a parabolic PDE such as the cubic-quintic complex Ginzburg-Landau equation (1) considered in this paper. Here, matrix $\mathbf{L}$ has a large norm or an unbounded spectrum which makes the system stiff, and thus an ordinary explicit method is forced to use small time steps to ensure a stable integration process. To cope with the stiffness, exponential integrators treat the linear part exactly and approximate the nonlinear part by expansion series. There are basically three different types of single-step exponential integrators: integrating factor (IF) methods, exponential Runge-Kutta (ERK) methods, and split-step (SS) methods. Both IF and ERK are Runge-Kutta based methods, which fill out Butcher tables by exponential functions of $\mathbf{L}$ and derive order conditions similar to ordinary Runge-Kutta methods. ERK methods are very popular in the applied mathematics community. See ref. [24] for a detailed survey of IF and ERK methods. SS methods split a single integration step into several substeps and integrate the system forward by considering only the linear and nonlinear effect interchangeably. SS methods are frequently used in physics community especially in the field of nonlinear optics. In this section, we give an overview of IF, ERK and SS methods applied to solving semilinear problem (4).

2.1 IF methods

The integrating factor (IF) method is also called Lawson method [25] which dates back to 1967. It alleviates the stiffness of the linear part in (4) by a change of variables, i.e., $z(t)=e^{-t\mathbf{L}}y(t)$, resulting in a nonlinear system of the new variable

[TABLE]

Here, $e^{-t\mathbf{L}}$ is the integrating factor. This transformation effectively stabilizes this system since the linear part is integrated explicitly and the Jacobian of the new system

[TABLE]

has the same spectrum as that of the original one ${\partial N(t,y)}/{\partial y}$, which is assumed to be nonstiff. Then popular time-stepping schemes are free to be applied to solve the new system (5). Lawson [25] first used the 4th-order Runge-Kutta method to integrate the new system, after whom, various single-/multiple-step schemes have been implemented for (5). (See [26] for all these variants.) In this paper, we apply a general $s$-stage explicit Runge-Kutta scheme $Y_{i}=y_{n}+h\sum_{j=1}^{i-1}a_{ij}f(t_{n}+c_{j}h,Y_{j})$ with $y_{n+1}=y_{n}+h\sum_{i=1}^{s}b_{i}f(t_{n}+c_{i}h,Y_{i})$ to equation (5), and obtain

[TABLE]

where $\beta_{i}=1-c_{i}$ and $\alpha_{ij}=c_{i}-c_{j}$. The local error is estimated by embedding a lower-order Runge-Kutta method in (6). In Sect. 3, we consider two representative stepsize-adaptive IF schemes.

In the quantum mechanics community, IF is referred to as Runge-Kutta methods in the interaction picture and was first used to integrate the time-dependent Gross-Pitaevskii equation by Caradoc-Davies [27]. The transformation used in such case was $z(t)=e^{-(t-t^{\prime})\mathbf{L}}y(t)$ with $t^{\prime}=t_{n}+\frac{1}{2}h_{n}$ during period $t_{n}$ to $t_{n+1}$. Here $h_{n}=t_{n-1}-t_{n}$. The name originates form the exponential transformation from the Schrödinger picture to the interaction picture and implies that the latter will simplify the calculation. Later, this transformation is applied to the (generalized) nonlinear Schrödinger equation by Johan Hult [28]. He compared it with SS methods and made it well recognized in the optical fibers community.

We note here that IF methods have a disadvantage of producing large error coefficients when the linear term $\mathbf{L}$ has a large norm, and they do not preserve fixed points of the original system (4). (For an improvement to generalize IF methods, see [29].) However, even with these drawbacks, IF methods have a merit of an easy implementation, and we consider IF methods for cubic-quintic complex Ginzburg-Landau equations.

2.2 ERK methods

Originally used in computational electrodynamics [30], exponential Rung-Kutta (ERK) methods are widely used in other fields of physics [31, 32, 33] and often referred to as exponential time differencing methods [34, 35]. An ERK method integrates the linear part exactly as IF methods. However, instead of using a change of variables, it resorts to an exact integration of (4) from $t_{n}$ to $t_{n+1}$ by the variation-of-constants formula

[TABLE]

Here, $h=t_{n+1}-t_{n}$ is the time step. In order to approximate the integral part above, a polynomial interpolation at $s$ non-confluent quadrature nodes $c_{1},\cdots,c_{s}$ is applied to $\mathbf{N}(t_{n}+\theta h,y(t_{n}+\theta h))$, i.e.,

[TABLE]

So we get

[TABLE]

Similar to an explicit Runge-Kutta method, $y(t_{n}+c_{i}h)$ can be approximated as a combination of $y(t_{n}+c_{j}h)$ for $j<i$. Then an ERK method becomes

[TABLE]

This can also be obtained by replacing $a_{ij}e^{h\alpha_{ij}\mathbf{L}}$ and $b_{i}e^{h\beta_{i}\mathbf{L}}$ in (6) with functions $a_{ij}(h\mathbf{L})$ and $b_{i}(h\mathbf{L})$. Here, coefficient functions $a_{ij}(h\mathbf{L})$ and $b_{i}(h\mathbf{L})$ are chosen as linear combinations of functions $\varphi_{j}(h\mathbf{L})$,

[TABLE]

Here $b_{i}(h\mathbf{L})$ in (8) is indeed a linear combination of $\varphi_{j}(h\mathbf{L})$. Furthermore, to gain a quick intuition for the choice of $\varphi_{j}(h\mathbf{L})$, we can expand the nonlinear part in (7) with respect to $\theta$ and get

[TABLE]

Here, $\mathbf{N}^{(r)}$ is the $r$th derivative of $\mathbf{N}$ with respect to $\theta$. If we define $\varphi_{0}(z)=e^{z}$, then $\varphi_{j}(z)$ has a recursion relation

[TABLE]

with the first few terms

[TABLE]

Order conditions of ERK methods were developed in two directions of nonstiff order conditions and stiff order conditions. Friedli [36] first derived nonstiff order conditions up to order 5 by matching Taylor series expansions of the exact and the numerical solutions, which were later extended to an arbitrary order by using B-series [37]. By expanding $a_{ij}(h\mathbf{L})$ and $b_{i}(h\mathbf{L})$ in power series of $h\mathbf{L}$ and truncating them to a certain order, nonstiff order conditions are obtained by matching the coefficients of $h$ with the same order on both sides of (9).

However, this process is limiting when $h\mathbf{L}$ has a large norm, which implies nonstiff order conditions are blind to the stiffness of the problem. To account for the stiffness, stiff order conditions up to order 4 were first given in [38]. Luan and Ostermann [39] gave the 5th-order conditions by a perturbation analysis after reformulating the scheme as a perturbation of the exponential Euler method. Later, the authors generalized the stiff order conditions to an arbitrary order by using exponential B-series [40]. One of the benefits of stiff order conditions is that they preserve the equilibria of the original system [38]. Stiff order conditions up to order 3 are listed in table 1.

Given a scheme with a certain nonstiff order, order reduction  [38] may appear if this scheme has a lower stiff order, However, stiff order conditions are rather restrictive, and under favorable conditions [38], schemes can show a higher order of convergence than the order predicted by the general stiff order conditions. For quite a few physically interesting PDEs [35], order reduction does not appear. The order behavior of an ERK method applied to a specific system is subtle. Thus, in this paper, we sort ERK methods by their nonstiff orders, and formulate stepsize-adaptive schemes whose stiff orders match their nonstiff orders.

2.3 SS methods

The main idea of split-step (SS) methods is that if the velocity field of a physical system can be decomposed as a sum of several separable sub-processes, then integration in one step can be approximated by several consecutive substeps. In each of the substeps only one sub-process takes effect. SS methods were first proposed in the 1950s by Bagrinovskii and Godunov [41]. It was also formulated by Strang [42] as an alternating-direction difference scheme, which has been widely used in integrating Hamiltonian systems [43] and PDEs of semilinear type [35].

To solve equation (4), we split the velocity field into one linear and another nonlinear part. For an $s$-stage SS method, one step of integration is decomposed into $2s$ substeps as follows,

[TABLE]

here, $\circ$ is a composition operator, which means that the integration result of one sub-process is the input to the next sub-process. $\phi_{\mathbf{N}}$ denotes the integration operator induced only by the nonlinear part:

[TABLE]

The local error in each step can be expressed in terms of commutator $[\mathbf{L},\mathbf{N}]=\mathbf{L}\mathbf{N}-\mathbf{N}\mathbf{L}$. For the order condition theory of SS methods, see [43, 44, 45, 46].

3 Embedded exponential integrators

In this paper, we explore time adaptive versions of numerical schemes considered in the previous sections for cubic-quintic complex Ginzburg-Landau equations. Numerical schemes such as SS [47, 8], Adams-Bashforth [48], ERK [49], Runge-Kutta in interaction picture [28] have been applied to integrate cubic-quintic complex Ginzburg-Landau equations. However, constant time-stepping schemes are not efficient to integrate pulsating or exploding soliton solutions which have different time scales as indicated in figure 1. In literature, step doubling and embedded methods are the two main approaches in stepsize control for ordinary Runge-Kutta methods. For exponential integrators, performance is a stringent concern and thus embedded methods are more preferable with its less induced overhead. Guided by this idea, Whalen, Brio and Moloney [50] incorporated time-step adaptation into several IF and ERK schemes by embedding lower-order schemes. W. Auzinger and his authors [51, 45, 52] made time-step adaptation possible in an SS method by embedding a lower- or same-order method.

In this section, we explore and introduce seven representative embedded schemes. IF4(3) and IF5(4) are IF based methods, where the two numbers $a(b)$ indicate that the scheme is $a$th-order accurate with an embedded $b$th-order scheme. ERK4(3)2(2), ERK4(3)3(3), ERK4(3)4(3) and ERK5(4)5(4) are ERK based methods, where four numbers $a(b)c(d)$ meant that the scheme has nonstiff order $a$ and stiff order $c$, and the embedded scheme has nonstiff order $b$ and stiff order $d$. For these six IF or ERK based schemes, we follow the first-same-as-last (FSAL) rule to embed lower-order schemes. That is, the last stage is evaluated at the same point as the first stage of the next step. The last scheme is SS4(3) which is based on an SS method which is 4th-order accurate and there is an embedded 3rd-order scheme.

3.1 IF4(3) and IF5(4)

To estimate the local integration error in a general $s$-stage IF method described by formula (6), we embed a lower-order scheme of form $\bar{y}_{n+1}=e^{h\mathbf{L}}y_{n}+h\sum_{i=1}^{s+1}\bar{b}_{i}e^{h\beta_{i}\mathbf{L}}\mathbf{N}(t_{n}+c_{i}h,Y_{i})$. The estimated local error can be expressed as

[TABLE]

Here, the embedded scheme has one more stage, that the summation is from $1$ to $s+1$. We consider two classical embedded Runge-Kutta schemes proposed by Dormand and Princes [53, 54] whose Butcher tables are tuned to minimize the truncation error coefficients. One is 4th-order accurate and the other is 5th-order accurate, both have one-order-lower embedded schemes. Their corresponding IF schemes are IF4(3) in table 2 and IF5(4) in table 3. We note that IF4(3) was used by Balac and Mahé [55] in the interaction picture context. Following FSAL rule, in their Butcher tables, the last intermediate stage is the same as the stage of evaluating $y_{n+1}$. The expressions for their local error estimation are listed in table 8.

As mentioned in [50], caution should be taken when ordinary embedded Runge-Kutta methods are imported into IF methods. If $\alpha_{ij}<0$ in (6), backward propagation in the intermediate stage $e^{\alpha_{ij}\mathbf{L}}$ is troublesome for some systems with unbounded negative linear parts.

3.2 ERK4(3)2(2), ERK4(3)3(3), ERK4(3)4(3) and ERK5(4)5(4)

We consider four different embedded ERK methods. The first one is ERK4(3)2(2) [34] proposed by Cox and Matthews and later improved by Kassam and Trefethen [35]. This scheme has nonstiff order 4 and stiff order 2, and the embedded scheme has nonstiff order 3 and stiff order 2. We embed the lower-order scheme as shown in table 4. Here, $\varphi_{i}=\varphi_{i}(h\mathbf{L})$ are the basis functions (10) and $\varphi_{i,j}=\varphi_{i}(c_{j}h\mathbf{L})$. As with embedded IF methods, ERK4(3)2(2) follows the FSAL rule.

The second one is ERK4(3)3(3) proposed by Krogstad [29] and is shown in table 5. Its Butcher table is slightly different from that of ERK4(3)2(2) but with it gives better convergence and stability. This scheme has nonstiff order 4 and stiff order 3. Note, the stiff order of neither ERK4(3)2(2) and ERK4(3)3(3) matches its nonstiff order, but they are very popular and for moderate stiff systems. It is also claimed [26] that it is hard to do much better than these two methods. Therefore, we consider these schemes for comparison. The embedded scheme in ERK4(3)3(3) is nonstiff 3rd order, stiff 3rd order.

The third scheme ERK4(3)4(3) is formulated by Hochbruck and Ostermann [38], whose Butcher table is shown in table 6. It is both nonstiff and stiff 4th-order accurate. Since the last node coefficient $c_{5}$ is $1/2$ not $1$, FASL approach fails to embed a one-order-lower scheme in this case. Fortunately, we observe that $c_{2}=c_{3}=c_{5}$, so by setting

[TABLE]

with other $b_{i}$ unchanged and choosing appropriate $x,y$, we hope to embed a 3rd-order scheme. Table 1 lists 5 stiff order conditions that should be satisfied in order to obtain a 3rd-order scheme. The 3rd stiff order condition is already satisfied. Setting $x+y=1$ ensures the 1st, 2nd and 4th stiff order conditions. Finally, setting $x=y=1/2$ make the embedded scheme satisfy a weakened but sufficient [38] 5th stiff order condition. The embedded scheme has stiff order 3. We verify that it is also nonstiff order 3. Note, compared to the FSAL approach, this embedding does not require one additional internal stage; thus saves one evaluation of the nonlinear function $\mathbf{N}(t,y)$.

The last scheme we consider is ERK5(4)4(4) formulated by Luan and Ostermann [56] shown in table 7. It is both nonstiff and stiff 5th-order accurate. Following the FSAL rule, we embed a nonstiff and stiff 4th-order scheme. The general nonstiff order conditions are given in [37] using bi-colored trees, but conditions only up to order 4 are listed explicitly. In this paper, we derive the 5th-order conditions which is presented in A for the completeness of the discussion.

For time adaptive method, from the Butcher table which consists of matrix functions of $h\mathbf{L}$, it seems that the stepsize-adaptive strategy is not efficient due to the cost associated with refilling the Butcher table every time the step size $h$ is updated. However, for cubic-quintic complex Ginzburg-Landau equations, the linear part $\mathbf{L}$ is diagonal, thus evaluation of $\varphi_{j}(h\mathbf{L})$ becomes an arithmetic calculation, which has linear complexity. Even for systems with non-diagonal linear parts, techniques such as Krylov-subspace methods can be deployed to accelerate matrix function evaluation. (See  [24] for more details.)

Another implementation issue associated with ERK methods is that direct evaluation of $\varphi_{j}(z)$ in (11) suffers from loss of accuracy when $z$ is small. It is believed that the contour integral method proposed by Kassam and Trefethen [35] can resolve this problem effectively. We take this approach in our implementations.

3.3 SS4(3)

We introduce one representative embedded SS scheme in this subsection. A lower-order scheme can be embedded in (12) by using a different set of coefficients $a_{i}$ and $b_{i}$. Unlike IF or ERK methods, stepsize-adaptive SS methods do not require recalculation of any run-time coefficients, so time-step adaptation can be implemented with nearly no additional cost. Recently, Auzinger and his coauthors [51, 45, 52] proposed and optimized over 30 different embedded SS schemes with real and complex coefficients $a_{i}$ and $b_{i}$. Four different strategies are suggested to estimate the local integration error, among which, the palindromic-pair strategy tends to have minimal local integration error. See [45] for the details. Here, we focus on the palindromic-pair scheme,

[TABLE]

[TABLE]

[TABLE]

The name comes from equation (14) which says that coefficients $b_{i}$ are totally determined by $a_{i}$, i.e., $b_{i}=a_{s+1-i}$. States $y_{1}$ and $y_{2}$ mirror each other by switching roles of linear and nonlinear operators. They approximate $y(t_{n+1})$ with the same order of accuracy and their leading coefficients of the local error have the same magnitude but opposite signs, so $y_{n+1}$ is the local extrapolation of $y_{1}$ and $y_{2}$ with one more order of accuracy. $\bar{y}_{n+1}$ serves as an error estimator.

In this paper, we focus on SS4(3), a $3$-stage palindromic-pair scheme with real coefficients, which is 4th-order accurate with an embedded 3rd-order scheme. Auzinger called it PP3/4A in [52]. The coefficients of SS4(3) are

[TABLE]

The linear part is integrated exactly in an SS scheme (12), and we solve (13) $s$ times during each single step. For some systems such as the cubic complex Ginzburg-Landau equation and the nonlinear Schrödinger equation, (13) can be solved explicitly [21], but such explicit formula does not exist for the cubic-quintic complex Ginzburg-Landau equation. We use numerical schemes to solve (13), in particular, the 4th-order Runge-Kutta scheme. Since there are many evaluations of the nonlinear function $\mathbf{N}(t,y)$ in a single step in (15) and (16), we restrict ourselves to only considering SS4(3) in this paper.

4 Time-step adaptation and performance metrics

Each stepsize-adaptive scheme mentioned in sect. 3 consists of a higher-order scheme and an embedded lower-order scheme. Though the local error is estimated for the lower-order scheme, following the local extrapolation strategy, we take the higher-order scheme to integrate the system forward. Integration accuracy is locally maintained by rejecting the current step size if the estimated local error is larger than the specified tolerance or accepting the state calculated by the current step size if otherwise. When the step size is rejected, a smaller step size is chosen to recompute the next state. If the calculated state is accepted then step size is scaled up for future computation. Let rtol be the relative tolerance for the local error, and we maintain

[TABLE]

at each integration step. Here, $E_{n+1}$ is the estimated local error for $y_{n+1}$. See expressions of $E_{n+1}$ for all seven schemes in table 8. Then the attempted new step size is

[TABLE]

$p$ is the order of the local error of the embedded scheme. $v$ is a safe factor and is set to $0.9$. Updating step size for each single step is not efficient because of the frequent recalculation of $h\mathbf{L}$ and other dependent coefficients. Also, to avoid step size oscillation, we update step size only when $s<1$ or when the difference between $h_{attempt}$ and $h$ is large enough. So we adopt the lazy adaptation strategy as in  [50]:

[TABLE]

The rule for updating step size is then $h_{new}=\mu h$.

In order to compare the performance of different schemes, the following metrics are used.

Nn: Number of evaluations of nonlinear function $\mathbf{N}(t,y)$ during the whole integration process. For complex systems, evaluations of $\mathbf{N}(t,y)$ take the majority part of the total integration time, thus we use Nn as one of the main metrics to compare different methods. nn denotes the number of evaluations of $\mathbf{N}(t,y)$ in a single step, whose values for seven different schemes are listed in table 8. Note that, in SS4(3), we use the 4th-order Runge-Kutta scheme to solve the nonlinear propagation equation (13). Thus its nn  entry is 24.

• 2.

Nab: Number of calculations of coefficients $a_{ij}$ or $b_{i}$ during the whole integration process. nab denotes the number of distinct $a_{ij}$ and $b_{i}$ entries in a Butcher table. Elements in Butcher tables of IF and ERK methods are exponential functions like $e^{h\mathbf{L}}$, which need to be recalculated whenever the step size is updated. Moreover, coefficients $a_{ij}$ and $b_{i}$ in ERK methods are evaluated by contour integrals, which need more time to calculate than those in IF methods. Thus we only consider evaluations of $a_{ij}$ and $b_{i}$ in ERK methods. Table 8 lists the nab values of four ERK schemes. Note, nab of ERK4(3)2(2) is 4 not 6 thanks to an implementation strategy from [35].

• 3.

Global relative accuracy: By using a very small step size, one can obtain a solution relatively close to the “true” solution. The global relative accuracy of each stepsize-adaptive scheme is then calculated.

• 4.

Wt: wall-clock time used for the integration, which is measured on a desktop equipped with 6 Intel i7-4930K 3.40GHz cores and 32G memory.

5 Numerical experiments and comparisons

We first show the performance of ERK4(3)2(2) for the three different soliton solutions displayed in Figure 1, and then compare the performance of the seven stepsize-adaptive schemes specifically for the exploding soliton solution. Figure 2 shows the integration results of ERK4(3)2(2). The spacing in these fence plots indicates the relative magnitudes of the step sizes. The pulsating soliton in panel (a) is integrated almost with a constant step size as in Figure 1(a), which is anticipated since there is only one time scale in the dynamics of this soliton. On the other hand, time-step adaptation slows down dramatically the integration of the high-spike parts for the extremely pulsating soliton in panel (b), and the performance of other six stepsize-adaptive schemes is similar for this soliton. This observation indicates that stepsize-adaptive methods efficiently integrate extremely pulsating solitons. For the exploding soliton in panel (c), ERK4(3)2(2) slows down the exploding parts moderately.

The initial condition that generates the exploding soliton in Figure 2(c) is

[TABLE]

which is a Gaussian wave in the middle of the domain composed with a small perturbation on the left side. Figure 3 shows the integration result in the form of heat map for the constant time-stepping method, ERK4(3)2(2) and SS4(3) respectively. Different from Figure 2(c), the heat maps scale the time axis (y-axis) to indicate the magnitude of step sizes. Two asymmetric explosions appear during the integration time window $t\in[0,20]$.

The explosions take small fractions of the total integration time. The soliton profile does not change for the rest. Figure 4 shows the estimated local integration error for Figure 3(a) and the step sizes used by ERK4(3)2(2) in Figure 3(b). During the fast-exploding parts, the estimated local error bursts substantially, spanning 2 to 3 orders of magnitude. The fast-exploding parts are the main cause of the accuracy lose. Step size is reduced when explosions happen and return to the normal level after explosions end. As shown in Figure 4(b), IF4(3), ERK4(3)2(2), ERK4(3)3(3) and ERK4(3)4(3) have almost the same adaptation pattern, while SS4(3) behaves more aggressively in the slow-moving parts. This is the reason for the explosions in Figure 3(c) are more stretched than those in Figure 3(b). Moreover, in Figure 4(b) the two holes of SS4(3) are slightly shifted to the left side compared to other 4th-order schemes, because SS4(3) uses large step sizes during slow-moving parts.

In Figure 3, the fast-exploding parts in panel (b) and (c) are stretched slightly compared with panel (a) and SS4(3) is than ERK4(3)2(2). We further explore this in sect. 6.

We compare the performance of the seven stepsize-adaptive schemes. Figure 5 shows the performance of all the seven schemes for integrating the exploding soliton. Panel (a) shows that rtol effectively controls the relative global error of all schemes. A smaller rtol usually produces a more accurate result. However, such tendency saturates for $\texttt{rtol}<10^{-11}$. SS4(3) has the largest relative global error since SS4(3) uses larger step sizes during slow-moving parts as shown in Figure 4(b). The integrator should spend more time on the fast-exploding parts, and with an appropriate rtol one can achieve a required global accuracy. Panel (c) shows the number of evaluations of the nonlinear function Nn versus rtol. The two 5th-order schemes IF5(4) and ERK5(4)5(4) have the least Nn because they use far fewer steps even though there are more evaluations of the nonlinear function in each step. Except for SS4(3), all other 4th-order schemes share a similar behavior of Nn. SS4(3) has a much larger Nn because there are 24 evaluations of $\mathbf{N}(t,y)$ in a single step for SS4(3). Panel (e) plots the wall time elapsed in seconds versus rtol. IF4(3), IF5(4) and SS4(3) have the same tendency as in panel (c). However, the relation saturates for ERK methods for a large rtol, which is most significant for ERK5(4)5(4). For an ERK method, the time used to refilling its Butcher table takes a larger percentage when rtol increases as shown in panel (f).

For ERK methods, refilling a Butcher table involves recalculation of $a_{ij}$ and $b_{j}$ which constitutes a large part of the total time. As shown in panel (b), Nab of ERK5(4)5(4) increases substantially when rtol gets beyond $10^{-11}$, and for ERK4(3)2(2), ERK4(3)3(3) and ERK4(3)4(3), Nab becomes slightly larger when rtol increases to $10^{-7}$. The time used for refilling Butcher tables dominates Wt at large rtol for ERK methods. This phenomenon raises a question: why does Nab increase when rtol becomes large enough? Nab is proportional to the number of times that a step size is rejected during the integration process. When rtol becomes larger, the step size is larger and thus there is an increased possibility of step-size oscillation in the integration process. Panel (d) shows that the number of rejections increases as rtol increases. The percentage of the time used to recalculate $h\mathbf{L}$-dependent coefficients increases, as shown in panel (f). For ERK methods, the time spent on refilling the Butcher table of ERK5(4)5(4) almost takes the whole computation time when rtol reaches $10^{-9}$ as shown in panel (f).

The stepsize-adaptive schemes effectively control the local error during the integration process. The two 5th-order schemes IF5(4) and ERK5(4)5(4) have the best performance for exploding solitons, but the performance of ERK5(4)5(4) deteriorates when rtol becomes too large.

6 Comoving-frame improvement for ERK methods

In Sect. 5, we show that stepsize-adaptive schemes slow down the integration of fast-exploding parts, but the slow-moving parts still take the majority of computation time as we compare panel (b), (c) with panel (a) in figure 3. Also, figure 4(b) shows that there is a plateau for the step sizes used in the slow-moving parts. What prevents one from using a larger step size is the fast-phase rotation of the complex field $A(x,t)$. Figure 6 shows the real part of the zeroth Fourier mode $\operatorname{Re}(a_{0})$ of $A(x,t)$ during the integration period in figure 3(a). The phase of $a_{0}$ rotates with a high frequency even though the profile $|A(x,t)|$ changes slowly in the slow-moving parts. Thus, if fast-phase rotation of complex field $A(x,t)$ can be handled effectively, one can further accelerate the slow-moving parts. We propose to use a comoving frame which has a similar rotating frequency as the original system to improve the results. In particular, we show that the performance of ERK methods can be improved in this frame work.

6.1 Dynamics in a comoving frame

To overcome the fast-phase rotation difficulty, we integrate the system in a comoving frame whose frequency is similar to the rotating frequency of $A(x,t)$. Let

[TABLE]

where $\Omega$ is the rotating frequency of this comoving frame. $A(x,t)$ and $\tilde{A}(x,t)$ is the state in the static and comoving frame respectively. Substituting (18) into the one-dimensional version of (1), we get

[TABLE]

By changing the real coefficient $\mu$ to the complex one $\mu-i\Omega$, we obtain the integrator in the comoving frame. Thus comoving frame introduces nearly no additional computational cost compared with the integrator in the static frame. To find the frequency of this comoving frame, one can simply measure the rotating frequency of $A(x,t)$ in the whole domain $x\in[0,L]$ for a certain slow-moving part, then obtain an average rotating rate. However, this approach is hard to automate and other issues such as the phase-wrapping effect, i.e., aliasing, can complicate this process.

In this paper, we utilize the underlining dynamically invariant structure of this exploding phenomenon. Figure 3(a) illustrates that the basic structure of the dynamics is a slow-moving soliton which undergoes intermittent explosions. If this soliton is viewed as a traveling wave, then an exploding part can be regarded as one homoclinic orbit of this traveling wave. We express an invariant solution of form

[TABLE]

Here, $A_{0}(x)=A(x,0)$ is a localized field. Constants $c$ and $\omega$ are spatial translation and phase velocity respectively. Definition (20) originates from the consideration of the two continuous symmetries of cubic-quintic complex Ginzburg-Landau equation, namely, equation (1) is invariant under spatial translation $A(x,t)\to A(x+\ell,t)$ and phase rotation $A(x,t)\to e^{i\phi}A(x,t)$. Soliton explosions are the result of the rapid growth of perturbations in the unstable directions of a traveling wave. The collapse of explosions is due to the dispersion effect. For more descriptive details, see [57]. Therefore, we set the frequency $\Omega$ of the comoving frame to the rotating frequency $\omega$ of this traveling wave, and integrate $\tilde{A}(x,t)$ instead of $A(x,t)$ for a better performance.

We find traveling waves in the Fourier mode space, in which equation (20) becomes

[TABLE]

where $q_{c}=2\pi c/L$. Taking time derivative on both sides, we get

[TABLE]

Here velocity field $\dot{a}_{k}$ is given in (2). Equation (21) defines an underdetermined system with respect to variables $(a_{k},c,\omega)$, whose roots are traveling waves. Given relatively good initial guesses, Newton-based methods converge quadratically to the traveling wave solution. In practice, we use Levenberg-Marquardt algorithm [58, 59] for its good performance in solving underdetermined systems. More details can be found in [57]. The traveling wave obtained is shown in figure 7(a) with

[TABLE]

This traveling wave lives in the symmetric subspace and thus has no spatial shift, but its phase rotates rapidly. By integrating in the comoving frame $\Omega=\omega$, we obtain $\operatorname{Re}(\tilde{a}_{0})$ shown in figure 7(b). $\tilde{a}_{0}$ is the zeroth Fourier mode of $\tilde{A}(x,t)$. Compared to figure 6, we see that the fast-phase rotation is effectively reduced for the slow-moving parts, while the explosion parts still have fast-phase dynamics.

We emphasize that the traveling wave found for one specific set of parameters of the cubic-quintic complex Ginzburg-Landau equation can be used as an initial guess to find traveling waves in the parameter space, so finding an appropriate frequency of the comoving frame for different parameters can be easily automated.

6.2 Performance in the comoving frame

Figure 8(a) shows the integration result of ERK4(3)2(2) in the comoving frame. Compared to figure 3(b)(c), the slow-moving part is sufficiently accelerated. Time steps used during integration process are shown in figure 8(b). There are several interesting observations comparing figure 8(b) with figure 4(b). First, for all seven methods, step sizes used during the fast-exploding parts are almost the same $h\simeq 10^{-4}$ in both static and comoving frame. This is reasonable because the comoving frame cannot reduce the rapid phase rotation during explosions.

Second, step sizes have increased substantially for ERK methods in the comoving frame for the slow-moving parts, while there is no change for IF4(3)4, IF5(4) and SS4(3). The comoving frame only promotes the performance of ERK methods, not IF or SS methods. This intuitively contradictory result comes from the difference between the intermediate states among different methods. The comoving frame changes the linear part from $\mathbf{L}$ to $\mathbf{L}-i\Omega$ in (19). The intermediate state (6) of an IF method

[TABLE]

and the two Butcher tables, table 2 and table 3, show that coefficients $e^{hc_{i}\mathbf{L}}$ and $a_{ij}e^{h\alpha_{ij}}$ shift by $e^{-i\Omega c_{i}h}$ and $e^{-i\Omega\alpha_{ij}h}$ respectively. Assume $Y_{j}$, $j<i$, has shift $e^{-i\Omega c_{j}h}$, then $a_{ij}e^{h\alpha_{ij}\mathbf{L}}\mathbf{N}(t_{n}+c_{j}h,Y_{j})$ has shift $e^{-i\Omega c_{i}h}$ because $\alpha_{ij}+c_{j}=c_{i}$. The intermediate state at $t_{n}+c_{i}h$ changes from $Y_{i}$ to $Y_{i}e^{-i\Omega c_{i}h}$. Such a phase change only introduces a phase shift for the local error estimation in table 8. Therefore, IF methods are invariant under $\mathbf{L}\to\mathbf{L}-i\Omega$. For SS methods, both (15) and (16) have the same phase shift, so transformation $\mathbf{L}\to\mathbf{L}-i\Omega$ only introduces a phase rotation for the local error estimation and thus does not change the behavior of SS methods either. However, for ERK methods (9), coefficients $a_{ij}(h\mathbf{L})$ and $b_{i}(h\mathbf{L})$ are functions of $\varphi_{j}(h\mathbf{L})$, which does not have an explicit phase rotation relation under transformation $\mathbf{L}\to\mathbf{L}-i\Omega$. So, intermediate state $Y_{i}$ is not transformed to $Ye^{-i\Omega c_{i}h}$. The comoving frame modifies the local error estimation for ERK methods. For the slow-moving parts, the rapid phase rotation is effectively reduced in the comoving frame and intermediate states tend to have smaller differences in phase. Figure 8(b) illustrates how comoving frame accelerates the integration of the slow-moving parts.

We repeat the numerical experiments of figure 5 in a comoving frame shown in Figure 9. The performance of IF4(3), IF5(4) and SS4(3) does not change in the comoving frame compared to that in the static frame. On the other hand, there are several differences for the ERK methods. Comparing figure 9(a) with figure 5(a), the global accuracy deteriorates in the comoving frame since larger step sizes are used for the slow-moving parts. Smaller rtol should be chosen in order to achieve the required global accuracy in the comoving frame for ERK methods. This is reasonable if one cares more about the percentage of time spent on the fast-exploding parts. Figure 9(c) shows that the number of evaluation of $\mathbf{N}(t,y)$ reduced significantly compared to the data in figure 5(c). The total integration time decreased as shown in figure 9(e). The number of times that a step size is rejected becomes 2 to 4 times larger in the comoving frame as shown in figure 9(d) compared to figure 5(d). For ERK5(4)5(4), the time for refilling the Butcher table takes the majority computation time even when rtol is as small as $10^{-13}$, shown in figure 9(f). For a large rtol, ERK5(4)5(4) is not as efficient as other schemes, as indicated by figure 9(e).

In summary, the performance of ERK methods improves in a comoving frame. To spend more integration time on the fast-exploding parts, one should use ERK4(3)2(2), ERK4(3)3(3), or ERK4(3)4(3), but not ERK5(4)5(4), since the last one spends too much time refilling its Butcher table.

7 Two-dimensional numerical experiments

In this section, we consider exploding solitons in two-dimensional cubic-quintic complex Ginzburg-Landau equation. All seven stepsize-adaptive schemes are exploited and compared in both static and comoving frame. The following initial condition, which represents a Gaussian wave at the center of the grid with a small perturbation in the southwest direction, is used to generate a sequence of exploding solitons.

[TABLE]

We integrate the system for time window $t\in[0,20]$, during which there are three explosions. The snapshots of one such explosion are shown in figure 10 (a)$\sim$(e). The profile of the soliton augments asymmetrically during this exploding process. Also, there are several sharp cracks at the center of the soliton in panel (c), which move and change from panel (c) to (d). All these fast-changing structures require a smaller step size to maintain a certain local integration accuracy compared to the step size needed before and after the explosion, i.e., panel (a) and (e) respectively.

Figure 11 shows the performance of the seven adaptive time-stepping schemes. Panel (a) shows the step sizes used during the integration. There are three dips in this figure, which represent that the integrators slow down during three exploding instances. IF5(4) and ERK5(4)5(4) use larger step sizes than those of 4th-order schemes. All 4th-order schemes use similar step sizes. Panel (b)$\sim$(f) show the performance measured by different metrics. The tendencies shown in panel (b)$\sim$(f) are very similar to those in figure 5. Among the four ERK methods, ERK5(4)5(4) has the largest number of times of calculating the elements of its Butcher table as shown in panel (b). Two facts account for this: First, there are 23 different elements in the Butcher table of ERK5(4)5(4), which are far more than other ERK methods. Second, as shown in panel (d), attempted trials of ERK5(4)5(4) are more likely to get rejected compare to other ERK methods. Panel (c) shows the number of evaluations of the nonlinear function. Similar to figure 5(c), the two 5th-order schemes have least evaluations of the nonlinear function. SS4(3) has the largest number of evaluations because it computes the nonlinear function 24 times in a single step. Panel (e) and (f) show the wall time of integration and the percentage of time spent on recalculating coefficients. The two-dimensional simulation is about 3- or 4-order more expensive than the one-dimensional simulation comparing panel (e) with figure 5(e). In the two-dimensional integration, ERK5(4)5(4) spends most of the time recalculating the coefficients. it is less efficient than IF5(4).

Similar to the one-dimensional case, explosions in the two-dimensional case can be visualized as homoclinic orbits of an unstable traveling wave of this system. A two-dimensional traveling wave has form

[TABLE]

with $c_{x}$, $c_{y}$, and $\omega$ respectively the translation velocity in the $x$, $y$ direction, and the phase velocity. By using a shooting method [60], we find a traveling wave living in the symmetric subspace with

[TABLE]

Figure 10(f) shows the profile of this traveling wave. For the slow-moving parts, the system has a similar phase-rotation rate as this traveling wave. Therefore, similar to (18) and (19), we can define a comoving frame for the two-dimensional cubic-quintic complex Ginzburg-Landau equation.

Figure 12 displays the performance of the seven stepsize-adaptive schemes in the comoving frame. As the one-dimensional case, only the performance of ERK schemes is affected by the comoving frame. Comparing panel (a) in figure 11 and that in figure 12, the step sizes used during the slow-moving parts are about 2 to 3 times larger for ERK schemes in the comoving frame, which is manifest for the part after the third explosion. Also, for ERK schemes in panel (c), the number of times to evaluate the nonlinear function decreased insignificantly. However, the rejection rates increased for ERK schemes as shown in panel (d). As a consequence, the total integration time only decreased by a small amount as shown in panel (e). The improvement of performance in the comoving frame is not striking. The reason is that the phase velocity (23) is not large; thus the fast-phase rotation problem is not severe in the static frame. Therefore, to integrate the two-dimensional cubic-quintic complex Ginzburg-Landau equation with the parameters chosen in this paper, the best choice out of the seven schemes is IF5(4).

8 Conclusions

We explored different stepsize-adaptive schemes to integrate soliton solutions in one- and two-dimensional cubic-quintic complex Ginzburg-Landau equation. We put an emphasis on the exploding solitons which have different time scales. The slow-moving part and the fast-exploding part should use different step sizes in order to integrate the system efficiently. By embedding lower order schemes in IF, ERK and SS methods, local integration error is controlled effectively. The step size is adapted to maintain a relative accuracy in each single integration step. For solitons with extremely pulsating amplitudes, time-step adaptation works well in the static frame. While for exploding solitons, to better handle the fast-phase rotation difficulty, we integrated the system in a comoving frame whose rotating frequency is similar to that of the soliton solution. We show that integration in the comoving frame can further accelerate the slow-moving parts for ERK methods.

In the one-dimensional case, the two 5th-order methods IF5(4) and ERK5(4)5(4) have the best performance. IF5(4) is easy to implement because it has a simple Butcher table structure. ERK5(4)5(4) can benefit from the comoving frame and remarkably slow down the fast-exploding parts. Since ERK5(4)5(4) spends a long time refilling its Butcher table when local error tolerance rtol becomes large, we prefer ERK5(4)5(4) when rtol is small but choose IF5(4) otherwise. In the two-dimensional case, we find that IF5(4) has the best performance. When the phase rotation is not very severe, a comoving frame may not improve the results much. When one need to implement a stepsize-adaptive scheme quickly, then the 4th-order schemes are suitable. ERK4(3)2(2) may be the best choice among all 4th-order methods considered in this paper because of its simple Butcher table structure.

The main focus of this paper is to explore different methods to experiment cubic-quintic complex Ginzburg-Landau equation. Nonetheless, these numerical schemes can be applied to other stiff (and non-stiff) systems that exhibit intermittent behaviors.

9 Acknowledgements

We are grateful to P. Cvitanović for providing insightful arguments about the phase rotation phenomenon in one-dimensional cubic-quintic complex Ginzburg-Landau equation. X.Ding is supported by a grant from G. Robinson, Jr.. S.H. Kang is supported by Simons Foundation grant 282311.

Appendix A

Table 9 displays the nonstiff 5th-order conditions. There are 37 equations. The first 26 come from adding a black or white node to the root node of the 4th-order bi-colored trees [37]. The remaining 11 equations are enumerated by counting the 2-fork type (8 equations), the 3-fork type (2 equations) and the 4-fork type (1 equation). For the general result to an arbitrary order, see ref. [37] (Theorem 2.1).