Stepsize-adaptive integrators for dissipative solitons in cubic-quintic complex Ginzburg-Landau equations
SStepsize-adaptive integrators for dissipative solitons in cubic-quintic complexGinzburg-Landau equations
X. Ding a, ∗ , S. H. Kang b a Center for Nonlinear Science, School of Physics, Georgia Institute of Technology, Atlanta, GA 30332-0430, USA b School of Mathematics, Georgia Institute of Technology, Atlanta, GA 30332-0430, USA
Abstract
This paper is a survey on exponential integrators to solve cubic-quintic complex Ginzburg-Landau equations andrelated sti ff problems. In particular, we are interested in accurate computation near the pulsating and explodingsoliton solutions where di ff erent time scales exist. We explore stepsize-adaptive variations of three types of ex-ponential 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 forcompleteness, and consider seven di ff erent stepsize-adaptive integrating schemes to solve the cubic-quintic complexGinzburg-Landau equation. Moreover, we propose using a comoving frame to resolve fast phase rotation for bet-ter performance. We present thorough comparisons and experiments in the one- and two-dimensional cubic-quinticcomplex Ginzburg-Landau equations. Keywords: stepsize-adaptive, cubic-quintic complex Ginzburg-Landau, dissipative solitons, exponential integrator,comoving frame
1. Introduction
The complex Ginzburg-Landau equation [1, 2] is one of the most frequently studied nonlinear equations in physicsand applied mathematics community. Derived as a general amplitude equation near the onset of instability in thecontext of pattern formation [1], it has applications in various fields of physics ranging from nonlinear optics [3] tosuperconductivity [4]. In order to preserve the invariance of the equation under a phase change A → Ae i φ , complexGinzburg-Landau equations only have odd-order nonlinear terms. The form with a cubic term is frequently used tostudy spatiotemporal chaos and intermittent traveling waves. Recently, complex Ginzburg-Landau equation with anadditional quintic term A t = µ A + ( D r + iD i ) ∆ A + ( β r + i β i ) | A | A + ( γ r + i γ i ) | A | A (1)has attracted attention for its peculiar dissipative soliton solutions in one- and two-dimensional cases [3, 5–9]. Here, A ( x , t ) or A ( x , y , t ) is a complex field. ∆ stands for the Laplace operator ∂ xx or ∂ xx + ∂ yy . A dissipative soliton is aself-localized structure that arises in spatially extended dissipative systems. It maintains its shape temporally similarto a constantly propagating solitary wave packet in an optical medium. Dissipative solitons appear as a result of thebalance between dispersion and nonlinearity and the balance between the gain and loss of energy. Dispersion spreadsthe shape while nonlinearity focuses it. A nontrivial internal energy flow and the dependence on an external energysource di ff erentiate dissipative solitons from solitons in integrable systems. Apart from a phase shift, solitons in inte-grable systems remain unchanged during soliton-soliton interactions; while dissipative solitons are free from energyor momentum conservation during scattering and annihilation. For a more comprehensive discussion on dissipativesolitons, see ref. [10]. While dissipative solitons lack most of the very special properties of soliton solutions of inte-grable Hamiltonian systems, they are generic and physically important, both because such structures are observed for ∗ Corresponding author
Email addresses: [email protected] (X. Ding), [email protected] (S. H. Kang)
Preprint submitted to Computer Physics Communications November 1, 2017 a r X i v : . [ phy s i c s . c o m p - ph ] O c t ide ranges of equation parameters [11, 12], and because many types of (pulsating) soliton solutions are observed inlaser experiments [13–15].By exploring the parameter space of (1), Akhmediev et al. [3, 6] identified several types of dissipative soliton so-lutions, 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 insoliton solutions in both one- and two-dimensional cubic-quintic complex Ginzburg-Landau equation, and the studyof their bifurcations. For example, Soto-Crespo, Akhmediev and Chiang [7] found two coexisting solitons with a highand a low energy respectively. Chang, Soto-Crespo, Vouzas and Akhmediev [16, 17] found a new class of pulsatingsolitons with large ratios of maximal to minimal energies as shown in figure 1(b). The energy may change more thantwo orders of magnitude in each period. Tsoy and Akhmediev [18] studied bifurcations from stationary to pulsatingsolitons based on reducing the infinite-dimensional dynamics (1) to a five-dimensional model. Meanwhile, Mancasand Choudhury [19] obtained a three-dimensional model of (1) by a variational method in the study of pulsating andsnake solitons. Among all the discoveries, solitons which undergo intermittent explosions stand out. As shown infigure 1(c), an exploding soliton moves uniformly for the most time with occasional substantial changes, after whichit quickly restores back to the pre-explosion profile. Both symmetric and asymmetric explosions [20] are recorded andthe 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 di ff usionrate [8, 9], which is unexpected for a deterministic system. x t (a) x t (b) x t | A | (c) Figure 1: Di ff erent soliton profiles. The z -coordinate represents the magnitude of the field | A ( x , t ) | . In call cases, µ = − . D i = . γ r = − . D r = . β r = . β i = γ i = − .
08. (b) A soliton with extremely pulsating amplitudes. D r = . β r = β i = . γ i = −
5. (c) An exploding soliton. D r = . β r = β i = . γ i = − . 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 e ffi ciently integrate dissipative solitons. There is little concernabout solitons that propagate uniformly such as plain or composite solitons. But for pulsating or exploding solitonswhose profiles change periodically or intermittently with interleaved fast and slow dynamics, the fast-changing partssuch as the high-amplitude parts in figure 1(b) and the exploding parts in figure 1(c) require smaller integration timesteps. Thus stepsize-adaptive integration schemes are preferred over constant time-stepping schemes. The purpose ofthis paper is to introduce several stepsize-adaptive schemes to integrate dissipative solitons in cubic-quintic complexGinzburg-Landau equation.For complex Ginzburg-Landau equations, finite di ff erence method [21], Fourier pseudo-spectral method [22] andChebyshev-Tau spectral method [23] are applied as spatial discretization methods. In this paper, we explore theFourier pseudo-spectral approach, for its spectral accuracy when the geometry of the solution is smooth and regular.Periodic boundary condition x ∈ [0 , L ] or ( x , y ) ∈ [0 , L x ] × [0 , L y ] is enforced for the one- and two-dimensional casesrespectively. In the Fourier space, equation (1) takes form of˙ a k = [ µ − ( D r + iD i ) q k ] a k + ( β r + i β i ) F k ( | A | A ) + ( γ r + i γ i ) F k ( | A | A ) (2)for the one-dimensional case, and˙ a mn = [ µ − ( D r + iD i )( q x m + q y n )] a mn + ( β r + i β i ) F mn ( | A | A ) + ( γ r + i γ i ) F mn ( | A | A ) (3)2or the two-dimensional case. Here, q = π/ L , q x = π/ L x and q y = π/ L y . a k and a mn are the k th Fourier modeof A ( x , t ) and mn th Fourier mode of A ( x , y , t ) respectively. F k ( · ) and F mn ( · ) denote the k th and mn th component ofthe one- and two-dimensional discrete Fourier transform. The linear part is represented by a diagonal matrix in thefrequency domain. The nonlinear part is first transformed back to the physical domain by the inverse fast Fouriertransform (FFT), integrated in the physical domain and then transformed again to the frequency domain. This is whythis strategy takes name “pseudo-spectral”. In soliton simulations, the size of the domain is set to L = L x = L y = Fourier modes are used respectivelyin (2) and (3). The number of Fourier modes has been doubled in numerical experiments but without noticeableimprovement in accuracy.For the time derivative, the linear parts of (2) and (3) have quadratic structures and thus are sti ff . Also, to modelpulsating and exploding soliton solution accurately, popular single-step integrators such as Runge-Kutta methods re-quires 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 di ff erent types of stepsize-adaptive exponential integrators.Our motivation is to tackle cubic-quintic complex Ginzburg-Landau equation, and to the best of our knowledge, thisis the first work to apply these methods to dissipative soliton integration in cubic-quintic complex Ginzburg-Landauequations. The main contributions of this paper are as follows. First, We not only convert two popular integratingfactor methods and four popular exponential Rung-Kutta methods to their stepsize-adaptive versions, but also considerone stepsize-adaptive split-step method with symmetric coe ffi cients. We note that the Runge-Kutta methods in theinteraction picture used in the quantum mechanics community are equivalent to the integrating factor methods thatwill be explained in sect. 2.1. Second, we formulate an embedded lower-order method in a 4th-order exponentialRunge-Kutta method in table 6 without an additional internal stage and are the first to study the embedded 5th-orderexponential Runge-Kutta method which will be introduced in sect. 3. Third, we utilize invariant dynamical structuresin 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 changeof 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 complexGinzburg-Landau equations e ff ectively.This paper is organized as follows. In sect. 2, we review three di ff erent types of exponential integrators: integratingfactor (IF) methods, exponential Rung-Kutta (ERK) methods and split-step (SS) methods. In sect. 3, we discuss howto embed a lower-order scheme in an exponential integrator and introduce seven representative schemes. The strategyto update step size and the metrics to evaluate the performance of an integrator is discussed in sect. 4. Sect. 5 is devotedto the numerical experiments on the one-dimensional cubic-quintic complex Ginzburg-Landau equation. We comparethe performance of constant time-stepping schemes with that of stepsize-adaptive schemes. In sect. 6, we introduce theidea of using a comoving frame to alleviate the fast-phase rotation problem. Sect. 7 is devoted to the discussion of theperformance 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 y (cid:48) ( t ) = f ( t , y ) = L y + N ( t , y ) , (4)which is usually spatially discretized from a parabolic PDE such as the cubic-quintic complex Ginzburg-Landauequation (1) considered in this paper. Here, matrix L has a large norm or an unbounded spectrum which makes thesystem sti ff , and thus an ordinary explicit method is forced to use small time steps to ensure a stable integrationprocess. To cope with the sti ff ness, exponential integrators treat the linear part exactly and approximate the nonlinearpart by expansion series. There are basically three di ff erent types of single-step exponential integrators: integrating3actor (IF) methods, exponential Runge-Kutta (ERK) methods, and split-step (SS) methods. Both IF and ERK areRunge-Kutta based methods, which fill out Butcher tables by exponential functions of L and derive order conditionssimilar to ordinary Runge-Kutta methods. ERK methods are very popular in the applied mathematics community. Seeref. [24] for a detailed survey of IF and ERK methods. SS methods split a single integration step into several substepsand integrate the system forward by considering only the linear and nonlinear e ff ect interchangeably. SS methods arefrequently used in physics community especially in the field of nonlinear optics. In this section, we give an overviewof IF, ERK and SS methods applied to solving semilinear problem (4). The integrating factor (IF) method is also called Lawson method [25] which dates back to 1967. It alleviates thesti ff ness of the linear part in (4) by a change of variables, i.e., z ( t ) = e − t L y ( t ), resulting in a nonlinear system of thenew variable z (cid:48) ( t ) = e − t L N ( t , e t L z ) . (5)Here, e − t L is the integrating factor. This transformation e ff ectively stabilizes this system since the linear part isintegrated explicitly and the Jacobian of the new system e − t L · ∂ N ( t , y ) ∂ y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) y = e t L z · e t L has the same spectrum as that of the original one ∂ N ( t , y ) /∂ y , which is assumed to be nonsti ff . 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 implementedfor (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 (cid:80) i − j = a i j f ( t n + c j h , Y j ) with y n + = y n + h (cid:80) si = b i f ( t n + c i h , Y i ) to equation (5), and obtain Y i = e hc i L y n + h i − (cid:88) j = a i j e h α ij L N ( t n + c j h , Y j ) , y n + = e h L y n + h s (cid:88) i = b i e h β i L N ( t n + c i h , Y i ) , (6)where β i = − c i and α i j = 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 wasfirst used to integrate the time-dependent Gross-Pitaevskii equation by Caradoc-Davies [27]. The transformation usedin such case was z ( t ) = e − ( t − t (cid:48) ) L y ( t ) with t (cid:48) = t n + h n during period t n to t n + . Here h n = t n − − t n . The name originatesform the exponential transformation from the Schr¨odinger picture to the interaction picture and implies that the latterwill simplify the calculation. Later, this transformation is applied to the (generalized) nonlinear Schr¨odinger equationby 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 coe ffi cients when the linear term L hasa large norm, and they do not preserve fixed points of the original system (4). (For an improvement to generalize IFmethods, see [29].) However, even with these drawbacks, IF methods have a merit of an easy implementation, and weconsider IF methods for cubic-quintic complex Ginzburg-Landau equations. Originally used in computational electrodynamics [30], exponential Rung-Kutta (ERK) methods are widely usedin other fields of physics [31–33] and often referred to as exponential time di ff erencing methods [34, 35]. An ERKmethod integrates the linear part exactly as IF methods. However, instead of using a change of variables, it resorts toan exact integration of (4) from t n to t n + by the variation-of-constants formula y ( t n + ) = e h L y ( t n ) + h (cid:90) e (1 − θ ) h L N ( t n + θ h , y ( t n + θ h )) d θ . (7)4ere, h = t n + − t n is the time step. In order to approximate the integral part above, a polynomial interpolation at s non-confluent quadrature nodes c , · · · , c s is applied to N ( t n + θ h , y ( t n + θ h )), i.e., N ( t n + θ h , y ( t n + θ h )) = s (cid:88) i = (cid:89) j (cid:44) i θ − c j c i − c j N ( t n + c i h , y ( t n + c i h )) . So we get y ( t n + ) = e h L y ( t n ) + h s (cid:88) i = b i ( h L ) N ( t n + c i h , y ( t n + c i h )) , b i ( h L ) = (cid:90) e (1 − θ ) h L (cid:89) j (cid:44) i θ − c j c i − c j d θ . (8)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 Y i = e hc i L y n + h i − (cid:88) j = a i j ( h L ) N ( t n + c j h , Y j ) , y n + = e h L y n + h s (cid:88) i = b i ( h L ) N ( t n + c i h , Y i ) . (9)This can also be obtained by replacing a i j e h α ij L and b i e h β i L in (6) with functions a i j ( h L ) and b i ( h L ). Here, coe ffi cientfunctions a i j ( h L ) and b i ( h L ) are chosen as linear combinations of functions ϕ j ( h L ), ϕ j ( z ) = (cid:90) e (1 − θ ) z θ j − ( j − d θ , j ≥ . Here b i ( h L ) in (8) is indeed a linear combination of ϕ j ( h L ). Furthermore, to gain a quick intuition for the choice of ϕ j ( h L ), we can expand the nonlinear part in (7) with respect to θ and get y ( t n + ) = e h L y ( t n ) + h ∞ (cid:88) r = h r N ( r ) (cid:90) e (1 − θ ) h L θ r r ! d θ = e h L y ( t n ) + h ∞ (cid:88) r = h r N ( r ) ϕ r + ( h L ) . Here, N ( r ) is the r th derivative of N with respect to θ . If we define ϕ ( z ) = e z , then ϕ j ( z ) has a recursion relation ϕ j + ( z ) = ϕ j ( z ) − / j ! z , ϕ j (0) = lim z → ϕ j ( z ) = j ! , j ≥ ϕ ( z ) = e z − z , ϕ ( z ) = e z − z − z , ϕ ( z ) = e z − z / − z − z . (11)Order conditions of ERK methods were developed in two directions of nonsti ff order conditions and sti ff orderconditions . Friedli [36] first derived nonsti ff order conditions up to order 5 by matching Taylor series expansions ofthe exact and the numerical solutions, which were later extended to an arbitrary order by using B-series [37]. Byexpanding a i j ( h L ) and b i ( h L ) in power series of h L and truncating them to a certain order, nonsti ff order conditionsare obtained by matching the coe ffi cients of h with the same order on both sides of (9).However, this process is limiting when h L has a large norm, which implies nonsti ff order conditions are blind tothe sti ff ness of the problem. To account for the sti ff ness, sti ff 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 asa perturbation of the exponential Euler method. Later, the authors generalized the sti ff order conditions to an arbitraryorder by using exponential B-series [40]. One of the benefits of sti ff order conditions is that they preserve the equilibriaof the original system [38]. Sti ff order conditions up to order 3 are listed in table 1.Given a scheme with a certain nonsti ff order, order reduction [38] may appear if this scheme has a lower sti ff order, However, sti ff order conditions are rather restrictive, and under favorable conditions [38], schemes can show a5 able 1: Sti ff order conditions for ERK up to order 3. Here, z = h L . J denotes an arbitrary square matrix. ψ j , i ( h L ) = (cid:80) i − k = a ik ( h L ) c j − k / ( j − − c ji ϕ j ( c i h L ). For the table up to order 5, see [39]. Order Index Sti ff order condition1 1 (cid:80) si = b i ( z ) = ϕ ( z )2 2 (cid:80) si = b i ( z ) c i = ϕ ( z )3 (cid:80) i − j = a i j ( z ) = c i ϕ ( c i z ) , i = , . . . , s (cid:80) si = b i ( z ) c i = ϕ ( z )5 (cid:80) si = b i ( z ) J ψ , i ( z ) = ff order conditions. For quite a few physicallyinteresting PDEs [35], order reduction does not appear. The order behavior of an ERK method applied to a specificsystem is subtle. Thus, in this paper, we sort ERK methods by their nonsti ff orders, and formulate stepsize-adaptiveschemes whose sti ff orders match their nonsti ff orders. The main idea of split-step (SS) methods is that if the velocity field of a physical system can be decomposed asa sum of several separable sub-processes, then integration in one step can be approximated by several consecutivesubsteps. In each of the substeps only one sub-process takes e ff ect. SS methods were first proposed in the 1950s byBagrinovskii and Godunov [41]. It was also formulated by Strang [42] as an alternating-direction di ff erence 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 SSmethod, one step of integration is decomposed into 2 s substeps as follows, y n + = φ N ( b s h ) ◦ e a s h L ◦ · · · ◦ φ N ( b h ) ◦ e a h L ◦ y n (12)here, ◦ is a composition operator, which means that the integration result of one sub-process is the input to the nextsub-process. φ N denotes the integration operator induced only by the nonlinear part: y (cid:48) ( t ) = N ( t , y ) . (13)The local error in each step can be expressed in terms of commutator [ L , N ] = LN − NL . For the order conditiontheory of SS methods, see [43–46].
3. Embedded exponential integrators
In this paper, we explore time adaptive versions of numerical schemes considered in the previous sections forcubic-quintic complex Ginzburg-Landau equations. Numerical schemes such as SS [8, 47], 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 e ffi cient to integrate pulsating or explodingsoliton solutions which have di ff erent time scales as indicated in figure 1. In literature, step doubling and embed-ded methods are the two main approaches in stepsize control for ordinary Runge-Kutta methods. For exponentialintegrators, performance is a stringent concern and thus embedded methods are more preferable with its less inducedoverhead. Guided by this idea, Whalen, Brio and Moloney [50] incorporated time-step adaptation into several IFand ERK schemes by embedding lower-order schemes. W. Auzinger and his authors [45, 51, 52] made time-stepadaptation 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 basedmethods, where the two numbers a ( b ) indicate that the scheme is a th-order accurate with an embedded b th-orderscheme. 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 nonsti ff order a and sti ff order c , and the embedded scheme has nonsti ff order b and sti ff order d . For these six IF or ERK based schemes, we follow the first-same-as-last (FSAL) rule to embed6ower-order schemes. That is, the last stage is evaluated at the same point as the first stage of the next step. The lastscheme is SS4(3) which is based on an SS method which is 4th-order accurate and there is an embedded 3rd-orderscheme. To estimate the local integration error in a general s -stage IF method described by formula (6), we embed alower-order scheme of form ¯ y n + = e h L y n + h (cid:80) s + i = ¯ b i e h β i L N ( t n + c i h , Y i ). The estimated local error can be expressed as E n + = ¯ y n + − y n + = h s + (cid:88) i = (¯ b i − b i ) e h β i L N ( t n + c j h , Y j ) . Here, the embedded scheme has one more stage, that the summation is from 1 to s +
1. We consider two classi-cal embedded Runge-Kutta schemes proposed by Dormand and Princes [53, 54] whose Butcher tables are tuned tominimize the truncation error coe ffi cients. One is 4th-order accurate and the other is 5th-order accurate, both haveone-order-lower embedded schemes. Their corresponding IF schemes are IF4(3) in table 2 and IF5(4) in table 3. Wenote that IF4(3) was used by Balac and Mah´e [55] in the interaction picture context. Following FSAL rule, in theirButcher tables, the last intermediate stage is the same as the stage of evaluating y n + . The expressions for their localerror estimation are listed in table 8. Table 2: Butcher table of IF4(3). Here, z = h L .
12 12 e z / e z / e z e z / e z / b i e z e z / e z / ¯ b i e z e z / e z / Table 3: Butcher table of IF5(4). Here, z = h L .
15 15 e z / e z /
10 940 e z / e z / − e z / e z /
289 193726561 e z / − e z /
45 644486561 e z / − e z / e z − e z / e z /
10 49176 e z / − e z / e z e z /
10 125192 e z / − e z / b i e z e z /
10 125192 e z / − e z / ¯ b i e z e z /
10 393640 e z / − e z / As mentioned in [50], caution should be taken when ordinary embedded Runge-Kutta methods are imported intoIF methods. If α i j < e α ij L is troublesome for some systemswith unbounded negative linear parts. We consider four di ff erent embedded ERK methods. The first one is ERK4(3)2(2) [34] proposed by Cox andMatthews and later improved by Kassam and Trefethen [35]. This scheme has nonsti ff order 4 and sti ff order 2, and7he embedded scheme has nonsti ff order 3 and sti ff order 2. We embed the lower-order scheme as shown in table 4.Here, ϕ i = ϕ i ( h L ) are the basis functions (10) and ϕ i , j = ϕ i ( c j h L ). As with embedded IF methods, ERK4(3)2(2)follows the FSAL rule. Table 4: Butcher table of ERK4(3)2(2). Here, ϕ i = ϕ i ( h L ), ϕ i , j = ϕ i ( c j h L ).
12 12 ϕ , ϕ , ϕ , ( ϕ , −
1) 0 ϕ , ϕ − ϕ + ϕ ϕ − ϕ ϕ − ϕ ϕ − ϕ b i ϕ − ϕ + ϕ ϕ − ϕ ϕ − ϕ ϕ − ϕ b i ϕ − ϕ + ϕ ϕ − ϕ ϕ − ϕ ϕ − ϕ The second one is ERK4(3)3(3) proposed by Krogstad [29] and is shown in table 5. Its Butcher table is slightlydi ff erent from that of ERK4(3)2(2) but with it gives better convergence and stability. This scheme has nonsti ff order4 and sti ff order 3. Note, the sti ff order of neither ERK4(3)2(2) and ERK4(3)3(3) matches its nonsti ff order, but theyare very popular and for moderate sti ff systems. It is also claimed [26] that it is hard to do much better than these twomethods. Therefore, we consider these schemes for comparison. The embedded scheme in ERK4(3)3(3) is nonsti ff ff Table 5: Butcher table of ERK4(3)3(3).
12 12 ϕ ,
212 12 ϕ , − ϕ , ϕ , ϕ , − ϕ , ϕ , ϕ − ϕ + ϕ ϕ − ϕ ϕ − ϕ ϕ − ϕ b i ϕ − ϕ + ϕ ϕ − ϕ ϕ − ϕ ϕ − ϕ b i ϕ − ϕ + ϕ ϕ − ϕ ϕ − ϕ ϕ − ϕ The third scheme ERK4(3)4(3) is formulated by Hochbruck and Ostermann [38], whose Butcher table is shown intable 6. It is both nonsti ff and sti ff ffi cient c is 1 / c = c = c , so by setting¯ b = xb , ¯ b = yb , ¯ b = b i unchanged and choosing appropriate x , y , we hope to embed a 3rd-order scheme. Table 1 lists 5 sti ff order conditions that should be satisfied in order to obtain a 3rd-order scheme. The 3rd sti ff order condition is alreadysatisfied. Setting x + y = ff order conditions. Finally, setting x = y = / ffi cient [38] 5th sti ff order condition. The embedded scheme has sti ff order 3. We verify that it is also nonsti ff order 3. Note, compared to the FSAL approach, this embedding does notrequire one additional internal stage; thus saves one evaluation of the nonlinear function 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 bothnonsti ff and sti ff ff and sti ff ff order conditions are given in [37] using bi-colored trees, but conditions only up to order 4 are listedexplicitly. In this paper, we derive the 5th-order conditions which is presented in Appendix A for the completenessof the discussion.For time adaptive method, from the Butcher table which consists of matrix functions of h L , it seems that thestepsize-adaptive strategy is not e ffi cient due to the cost associated with refilling the Butcher table every time thestep size h is updated. However, for cubic-quintic complex Ginzburg-Landau equations, the linear part L is diagonal,thus evaluation of ϕ j ( h L ) becomes an arithmetic calculation, which has linear complexity. Even for systems withnon-diagonal linear parts, techniques such as Krylov-subspace methods can be deployed to accelerate matrix function8 able 6: Butcher table of ERK4(3)4(3).
12 12 ϕ ,
212 12 ϕ , − ϕ , ϕ , ϕ , − ϕ , ϕ , ϕ ,
412 12 ϕ , − a , − a , a , a , ϕ , − a , b i ϕ − ϕ + ϕ − ϕ + ϕ ϕ − ϕ ¯ b i ϕ − ϕ + ϕ ϕ − ϕ ϕ − ϕ − ϕ + ϕ a , = ϕ , − ϕ , + ϕ , − ϕ , evaluation. (See [24] for more details.)Another implementation issue associated with ERK methods is that direct evaluation of ϕ j ( z ) in (11) su ff ers fromloss of accuracy when z is small. It is believed that the contour integral method proposed by Kassam and Trefethen [35]can resolve this problem e ff ectively. We take this approach in our implementations. We introduce one representative embedded SS scheme in this subsection. A lower-order scheme can be embeddedin (12) by using a di ff erent set of coe ffi cients a i and b i . Unlike IF or ERK methods, stepsize-adaptive SS methodsdo not require recalculation of any run-time coe ffi cients, so time-step adaptation can be implemented with nearly noadditional cost. Recently, Auzinger and his coauthors [45, 51, 52] proposed and optimized over 30 di ff erent embeddedSS schemes with real and complex coe ffi cients a i and b i . Four di ff erent strategies are suggested to estimate the localintegration 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,( a , b , . . . , a s , b s ) = ( b s , a s , . . . , b , a ) , (14) y = φ N ( b s h ) ◦ e a s h L ◦ · · · ◦ φ N ( b h ) ◦ e a h L ◦ y ( t n ) , (15) y = e b s h L ◦ φ N ( a s h ) ◦ · · · ◦ e b h L ◦ φ N ( a h ) ◦ y ( t n ) , (16) y n + =
12 ( y + y ) , ¯ y n + =
12 ( y − y ) . The name comes from equation (14) which says that coe ffi cients b i are totally determined by a i , i.e., b i = a s + − i . States y and y mirror each other by switching roles of linear and nonlinear operators. They approximate y ( t n + ) with thesame order of accuracy and their leading coe ffi cients of the local error have the same magnitude but opposite signs,so y n + is the local extrapolation of y and y with one more order of accuracy. ¯ y n + serves as an error estimator.In this paper, we focus on SS4(3), a 3-stage palindromic-pair scheme with real coe ffi cients, which is 4th-orderaccurate with an embedded 3rd-order scheme. Auzinger called it PP3 /
4A in [52]. The coe ffi cients of SS4(3) are( a , a , a ) = (0 . , − . , . . The linear part is integrated exactly in an SS scheme (12), and we solve (13) s times during each single step. Forsome systems such as the cubic complex Ginzburg-Landau equation and the nonlinear Schr¨odinger equation, (13)can be solved explicitly [21], but such explicit formula does not exist for the cubic-quintic complex Ginzburg-Landauequation. We use numerical schemes to solve (13), in particular, the 4th-order Runge-Kutta scheme. Since thereare many evaluations of the nonlinear function N ( t , y ) in a single step in (15) and (16), we restrict ourselves to onlyconsidering SS4(3) in this paper. 9 able 7: Butcher table of ERK5(4)5(4).
12 12 ϕ ,
212 12 ϕ , − ϕ , ϕ ,
314 14 ϕ , − ϕ , ϕ ,
412 12 ϕ , − ϕ , + ϕ , − ϕ , + ϕ , ϕ , − ϕ ,
515 15 ϕ , − ϕ , − a a , ϕ , − a ,
423 23 ϕ , + a − a − a − a a a ϕ , − a − a − a a a a ϕ − b − b − b b b b b i ϕ − b − b − b b b b b i ϕ − b − b − b b b b a = ϕ , − ϕ , , a = a − ϕ , + ϕ , , a = a + ϕ , − ϕ , ,φ = a − ϕ , + ϕ , − ϕ , + ϕ , + ϕ , + ϕ , , a = ϕ , − ϕ , − φ , a = − ϕ , + ϕ , + φ , a = − ϕ , + ϕ , + φ , b = ϕ − ϕ + ϕ , b = − ϕ + ϕ − ϕ , b = ϕ − ϕ + ϕ
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 extrapolationstrategy, we take the higher-order scheme to integrate the system forward. Integration accuracy is locally maintainedby rejecting the current step size if the estimated local error is larger than the specified tolerance or accepting thestate calculated by the current step size if otherwise. When the step size is rejected, a smaller step size is chosen torecompute 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 || E n + || ∞ < rtol · || y n + || ∞ at each integration step. Here, E n + is the estimated local error for y n + . See expressions of E n + for all seven schemesin table 8. Then the attempted new step size is h attempt = s · h , s = v (cid:32) rtol · || y n + || ∞ || E n + || ∞ (cid:33) / p . 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 foreach single step is not e ffi cient because of the frequent recalculation of h L and other dependent coe ffi cients. Also, toavoid step size oscillation, we update step size only when s < ff erence between h attempt and h is large10nough. So we adopt the lazy adaptation strategy as in [50]: µ = . s < . s . ≤ s < . .
85 0 . ≤ s <
11 1 ≤ s < . s . ≤ s < s ≥ h new = µ h .In order to compare the performance of di ff erent schemes, the following metrics are used. • Nn : Number of evaluations of nonlinear function N ( t , y ) during the whole integration process. For complexsystems, evaluations of N ( t , y ) take the majority part of the total integration time, thus we use Nn as one ofthe main metrics to compare di ff erent methods. nn denotes the number of evaluations of N ( t , y ) in a singlestep, whose values for seven di ff erent schemes are listed in table 8. Note that, in SS4(3), we use the 4th-orderRunge-Kutta scheme to solve the nonlinear propagation equation (13). Thus its nn entry is 24. • Nab : Number of calculations of coe ffi cients a i j or b i during the whole integration process. nab denotes thenumber of distinct a i j and b i entries in a Butcher table. Elements in Butcher tables of IF and ERK methodsare exponential functions like e h L , which need to be recalculated whenever the step size is updated. Moreover,coe ffi cients a i j and b i in ERK methods are evaluated by contour integrals, which need more time to calculatethan those in IF methods. Thus we only consider evaluations of a i j 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 strategyfrom [35]. • 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. • Wt : wall-clock time used for the integration, which is measured on a desktop equipped with 6 Intel i7-4930K3.40GHz cores and 32G memory. Table 8: Characteristics of the seven embedded schemes. ‘order’: the order of the local error of the embedded scheme. ‘ nn ’: the number ofevaluations of N ( t , y ) in a single step. ‘ nab ’: the number of distinct a ij and b i entries in a Butcher table. ‘ E n + ’: the expression for the estimatedlocal error. scheme order nn nab E n + IF4(3) 4 5 - h [ N ( t n + h , Y ) − N ( t n + h , Y )]IF5(4) 5 7 - h (cid:110) − N ( t n , Y ) + N ( t n + h , Y ) − N ( t n + h , Y ) + N ( t n + h , Y ) − N ( t n + h , Y ) + N ( t n + h , Y ) (cid:111) ERK4(3)2(2) 4 5 4 b h [ N ( t n + h , Y ) − N ( t n + h , Y )]ERK4(3)3(3) 4 5 8 b h [ N ( t n + h , Y ) − N ( t n + h , Y )]ERK4(3)4(3) 4 5 11 b h (cid:104) N ( t n + h , Y ) − ( N ( t n + h / , Y ) + N ( t n + h / , Y )) (cid:105) ERK5(4)5(4) 5 9 23 b h [ N ( t n + h , Y ) − N ( t n + h , Y )]SS4(3) 4 24 - ( y − y ) / . Numerical experiments and comparisons We first show the performance of ERK4(3)2(2) for the three di ff erent soliton solutions displayed in Figure 1, andthen 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 magni-tudes 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-stepadaptation 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 indicatesthat stepsize-adaptive methods e ffi ciently integrate extremely pulsating solitons. For the exploding soliton in panel(c), ERK4(3)2(2) slows down the exploding parts moderately. x t (a) x t (b) x t | A | (c) Figure 2: Integration of solitons in figure 1 by ERK4(3)2(2) with relative tolerance rtol = − . The spacing between consecutive profilesindicates the relative magnitude of integration step size. The initial condition that generates the exploding soliton in Figure 2(c) is A ( x , = . − (cid:32) xL − (cid:33) + . − (cid:32) xL − (cid:33) , (17)which is a Gaussian wave in the middle of the domain composed with a small perturbation on the left side. Figure 3shows the integration result in the form of heat map for the constant time-stepping method, ERK4(3)2(2) and SS4(3)respectively. Di ff erent from Figure 2(c), the heat maps scale the time axis (y-axis) to indicate the magnitude of stepsizes. Two asymmetric explosions appear during the integration time window t ∈ [0 , (a) x t (b) x t (c) x t Figure 3: Heat maps of integration of initial condition (17) for t ∈ [0 , | A ( x , t ) | . (a) constant time-steppingmethod with step size h = · − . (b) ERK4(3)2(2) and (c) SS4(3) both with relative tolerance rtol = − . fractions of the total integration time. The soliton profile does not change for the rest. Figure 4 shows the estimated12ocal integration error for Figure 3(a) and the step sizes used by ERK4(3)2(2) in Figure 3(b). During the fast-explodingparts, the estimated local error bursts substantially, spanning 2 to 3 orders of magnitude. The fast-exploding parts arethe main cause of the accuracy lose. Step size is reduced when explosions happen and return to the normal levelafter 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 thesame adaptation pattern, while SS4(3) behaves more aggressively in the slow-moving parts. This is the reason forthe explosions in Figure 3(c) are more stretched than those in Figure 3(b). Moreover, in Figure 4(b) the two holes ofSS4(3) are slightly shifted to the left side compared to other 4th-order schemes, because SS4(3) uses large step sizesduring slow-moving parts. In Figure 3, the fast-exploding parts in panel (b) and (c) are stretched slightly compared (a) t -17 -14 -11 -8 e s t i m a t e d l o c a l e rr o r IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3) (b) t -4 -3 h IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3)
Figure 4: (a) Estimated local errors for the constant time-stepping result in Figure 3(a). (b) Step sizes used during t ∈ [0 ,
20] for ERK4(3)2(2) inFigure 3(b). 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 theseven schemes for integrating the exploding soliton. Panel (a) shows that rtol e ff ectively controls the relative globalerror of all schemes. A smaller rtol usually produces a more accurate result. However, such tendency saturates for rtol < − . SS4(3) has the largest relative global error since SS4(3) uses larger step sizes during slow-moving partsas 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 fewersteps even though there are more evaluations of the nonlinear function in each step. Except for SS4(3), all other4th-order schemes share a similar behavior of Nn . SS4(3) has a much larger Nn because there are 24 evaluations of 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) andSS4(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 alarger percentage when rtol increases as shown in panel (f).For ERK methods, refilling a Butcher table involves recalculation of a i j and b j which constitutes a large part of thetotal time. As shown in panel (b), Nab of ERK5(4)5(4) increases substantially when rtol gets beyond 10 − , and forERK4(3)2(2), ERK4(3)3(3) and ERK4(3)4(3), Nab becomes slightly larger when rtol increases to 10 − . The timeused 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 stepsize is rejected during the integration process. When rtol becomes larger, the step size is larger and thus there is anincreased possibility of step-size oscillation in the integration process. Panel (d) shows that the number of rejectionsincreases as rtol increases. The percentage of the time used to recalculate h L -dependent coe ffi cients increases, asshown in panel (f). For ERK methods, the time spent on refilling the Butcher table of ERK5(4)5(4) almost takes thewhole computation time when rtol reaches 10 − as shown in panel (f).13 -13 -11 -9 -7 rtol -11 -8 -5 -2 r e l a t i v e e rr o r (a) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -13 -11 -9 -7 rtol N a b (b) ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4) -13 -11 -9 -7 rtol N n (c) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -13 -11 -9 -7 rtol nu m b e r o f r e j e c t i o n s (d) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3) -13 -11 -9 -7 rtol W t (e) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -13 -11 -9 -7 rtol -6 -4 -2 t i m e r a t i o (f) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3)
Figure 5: The performance of the seven stepsize-adaptive schemes for integrating the exploding soliton with initial condition Eq. (17).
The stepsize-adaptive schemes e ff ectively 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 ofERK5(4)5(4) deteriorates when rtol becomes too large. 14 . Comoving-frame improvement for ERK methods In Sect. 5, we show that stepsize-adaptive schemes slow down the integration of fast-exploding parts, but theslow-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 onefrom using a larger step size is the fast-phase rotation of the complex field A ( x , t ). Figure 6 shows the real part of thezeroth Fourier mode Re( a ) of A ( x , t ) during the integration period in figure 3(a). The phase of a rotates with a highfrequency even though the profile | A ( x , t ) | changes slowly in the slow-moving parts. Thus, if fast-phase rotation ofcomplex field A ( x , t ) can be handled e ff ectively, one can further accelerate the slow-moving parts. We propose to usea 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. t −300−200−1000100200300 R e ( a ) Figure 6: Real part of the zeroth Fourier mode of A ( x , t ) versus time. To overcome the fast-phase rotation di ffi culty, we integrate the system in a comoving frame whose frequency issimilar to the rotating frequency of A ( x , t ). Let A ( x , t ) = ˜ A ( x , t ) e i Ω t , (18)where Ω is the rotating frequency of this comoving frame. A ( x , t ) and ˜ A ( x , t ) is the state in the static and comovingframe respectively. Substituting (18) into the one-dimensional version of (1), we get˜ A t = ( µ − i Ω ) ˜ A + ( D r + iD i ) ˜ A xx + ( β r + i β i ) | ˜ A | ˜ A + ( γ r + i γ i ) | ˜ A | ˜ A . (19)By changing the real coe ffi cient µ to the complex one µ − i Ω , we obtain the integrator in the comoving frame. Thuscomoving 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 wholedomain x ∈ [0 , L ] for a certain slow-moving part, then obtain an average rotating rate. However, this approach is hardto automate and other issues such as the phase-wrapping e ff ect , 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 thistraveling wave. We express an invariant solution of form A ( x , t ) = A ( x + ct ) e i ω t . (20)15ere, A ( x ) = A ( x ,
0) is a localized field. Constants c and ω are spatial translation and phase velocity respec-tively. Definition (20) originates from the consideration of the two continuous symmetries of cubic-quintic complexGinzburg-Landau equation, namely, equation (1) is invariant under spatial translation A ( x , t ) → A ( x + (cid:96), t ) and phaserotation A ( x , t ) → e i φ A ( x , t ). Soliton explosions are the result of the rapid growth of perturbations in the unstabledirections of a traveling wave. The collapse of explosions is due to the dispersion e ff ect. For more descriptive details,see [57]. Therefore, we set the frequency Ω of the comoving frame to the rotating frequency ω of this traveling wave,and integrate ˜ 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 a k ( t ) = exp( i ω t + ikq c t ) a k (0) , k = , ± , ± , . . . , where q c = π c / L . Taking time derivative on both sides, we get˙ a k ( t ) − ikq c a k − i ω a k = , k = , ± , ± , . . . . (21)Here velocity field ˙ a k is given in (2). Equation (21) defines an underdetermined system with respect to variables( a k , c , ω ), whose roots are traveling waves. Given relatively good initial guesses, Newton-based methods convergequadratically to the traveling wave solution. In practice, we use Levenberg-Marquardt algorithm [58, 59] for its goodperformance in solving underdetermined systems. More details can be found in [57]. The traveling wave obtained isshown in figure 7(a) with c = , ω = . . This traveling wave lives in the symmetric subspace and thus has no spatial shift, but its phase rotates rapidly. Byintegrating in the comoving frame
Ω = ω , we obtain Re(˜ a ) shown in figure 7(b). ˜ a is the zeroth Fourier mode of˜ A ( x , t ). Compared to figure 6, we see that the fast-phase rotation is e ff ectively reduced for the slow-moving parts,while the explosion parts still have fast-phase dynamics. x | A | (a) t −300−200−1000100200300 R e ( ˜ a ) (b) Figure 7: (a) Profile of the traveling wave. (b) Real part of the zeroth Fourier mode of ˜ A ( x , t ) versus time. We emphasize that the traveling wave found for one specific set of parameters of the cubic-quintic complexGinzburg-Landau equation can be used as an initial guess to find traveling waves in the parameter space, so findingan appropriate frequency of the comoving frame for di ff erent parameters can be easily automated. Figure 8(a) shows the integration result of ERK4(3)2(2) in the comoving frame. Compared to figure 3(b)(c), theslow-moving part is su ffi ciently 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, step16izes used during the fast-exploding parts are almost the same h (cid:39) − in both static and comoving frame. This isreasonable because the comoving frame cannot reduce the rapid phase rotation during explosions. (a) x t (b) t -5 -4 -3 -2 h IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3)
Figure 8: (a) Spatiotemporal plot of | A ( x , t ) | produced by ERK4(3)2(2) in the comoving frame. (b) Step size h vs t with rtol = − in thecomoving frame. Second, step sizes have increased substantially for ERK methods in the comoving frame for the slow-movingparts, while there is no change for IF4(3)4, IF5(4) and SS4(3). The comoving frame only promotes the performanceof ERK methods, not IF or SS methods. This intuitively contradictory result comes from the di ff erence between theintermediate states among di ff erent methods. The comoving frame changes the linear part from L to L − i Ω in (19).The intermediate state (6) of an IF method Y i = e hc i L y n + h i − (cid:88) j = a i j e h α ij L N ( t n + c j h , Y j )and the two Butcher tables, table 2 and table 3, show that coe ffi cients e hc i L and a i j e h α ij shift by e − i Ω c i h and e − i Ω α ij h respectively. Assume Y j , j < i , has shift e − i Ω c j h , then a i j e h α ij L N ( t n + c j h , Y j ) has shift e − i Ω c i h because α i j + c j = c i .The intermediate state at t n + c i h changes from Y i to Y i e − i Ω c i h . Such a phase change only introduces a phase shiftfor the local error estimation in table 8. Therefore, IF methods are invariant under L → L − i Ω . For SS methods,both (15) and (16) have the same phase shift, so transformation L → L − i Ω only introduces a phase rotation for thelocal error estimation and thus does not change the behavior of SS methods either. However, for ERK methods (9),coe ffi cients a i j ( h L ) and b i ( h L ) are functions of ϕ j ( h L ), which does not have an explicit phase rotation relation undertransformation L → L − i Ω . So, intermediate state Y i is not transformed to Ye − i Ω c i h . The comoving frame modifiesthe local error estimation for ERK methods. For the slow-moving parts, the rapid phase rotation is e ff ectively reducedin the comoving frame and intermediate states tend to have smaller di ff erences in phase. Figure 8(b) illustrates howcomoving 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 ofIF4(3), IF5(4) and SS4(3) does not change in the comoving frame compared to that in the static frame. On the otherhand, there are several di ff erences for the ERK methods. Comparing figure 9(a) with figure 5(a), the global accuracydeteriorates in the comoving frame since larger step sizes are used for the slow-moving parts. Smaller rtol should bechosen in order to achieve the required global accuracy in the comoving frame for ERK methods. This is reasonableif one cares more about the percentage of time spent on the fast-exploding parts. Figure 9(c) shows that the numberof evaluation of N ( t , y ) reduced significantly compared to the data in figure 5(c). The total integration time decreasedas shown in figure 9(e). The number of times that a step size is rejected becomes 2 to 4 times larger in the comovingframe as shown in figure 9(d) compared to figure 5(d). For ERK5(4)5(4), the time for refilling the Butcher tabletakes the majority computation time even when rtol is as small as 10 − , shown in figure 9(f). For a large rtol ,17 -13 -11 -9 -7 rtol -11 -8 -5 -2 r e l a t i v e e rr o r (a) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -13 -11 -9 -7 rtol N a b (b) ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4) -13 -11 -9 -7 rtol N n (c) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -13 -11 -9 -7 rtol nu m b e r o f r e j e c t i o n s (d) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3) -13 -11 -9 -7 rtol W t (e) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -13 -11 -9 -7 rtol -6 -4 -2 t i m e r a t i o (f) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3)
Figure 9: Same as figure 5 but in the comoving frame.
ERK5(4)5(4) is not as e ffi cient 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 timeon 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), sincethe last one spends too much time refilling its Butcher table.18 . Two-dimensional numerical experiments (a) 12 (b) 12(c) 123 (d) 12(e) 12 (f) 12 Figure 10: (a) ∼ (e) Snapshots of one explosion in the two-dimensional cubic-quintic complex Ginzburg-Landau equation using initial condition(22). The system is integrated by ERK4(3)2(2). These five snapshots correspond to the state at t = , , , ,
12 respectively. (f) The profile of anunstable traveling wave in this system. In all figures, the color represents the magnitude of | A ( x , y , t ) | . System parameters are µ = − . D i = . γ r = − . D r = . β r = β i = . γ i = − . t -4 -3 -2 h (a) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3) -11 -9 -7 rtol N a b (b) ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4) -11 -9 -7 rtol N n (c) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -11 -9 -7 rtol nu m b e r o f r e j e c t i o n s (d) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3) -11 -9 -7 rtol W t (e) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -11 -9 -7 rtol -6 -4 -2 t i m e r a t i o (f) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3)
Figure 11: Performance of adaptive time-stepping schemes in the static frame for the two-dimensional cubic-quintic complex Ginzburg-Landauequation. (a) Step sizes used during the integration process when rtol = − . (b) ∼ (f) Performance measured by di ff erent metrics when rtol varies. In this section, we consider exploding solitons in two-dimensional cubic-quintic complex Ginzburg-Landau equa-tion. All seven stepsize-adaptive schemes are exploited and compared in both static and comoving frame. Thefollowing initial condition, which represents a Gaussian wave at the center of the grid with a small perturbation in the20 t -4 -3 -2 h (a) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3) -11 -9 -7 rtol N a b (b) ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4) -11 -9 -7 rtol N n (c) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -11 -9 -7 rtol nu m b e r o f r e j e c t i o n s (d) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3) -11 -9 -7 rtol W t (e) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3)ERK4(3)4(3)ERK5(4)5(4)SS4(3) -11 -9 -7 rtol -6 -4 -2 t i m e r a t i o (f) IF4(3)IF5(4)ERK4(3)2(2)ERK4(3)3(3) ERK4(3)4(3)ERK5(4)5(4)SS4(3)
Figure 12: The same as figure 11 but in a comoving frame. southwest direction, is used to generate a sequence of exploding solitons. A ( x , y , = . − (cid:32) xL − (cid:33) − (cid:32) yL − (cid:33) + . − (cid:32) xL − (cid:33) − (cid:32) yL − (cid:33) . (22)21e integrate the system for time window t ∈ [0 , ∼ (e). The profile of the soliton augments asymmetrically during thisexploding process. Also, there are several sharp cracks at the center of the soliton in panel (c), which move andchange from panel (c) to (d). All these fast-changing structures require a smaller step size to maintain a certain localintegration 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 usedduring the integration. There are three dips in this figure, which represent that the integrators slow down during threeexploding instances. IF5(4) and ERK5(4)5(4) use larger step sizes than those of 4th-order schemes. All 4th-orderschemes use similar step sizes. Panel (b) ∼ (f) show the performance measured by di ff erent metrics. The tendenciesshown in panel (b) ∼ (f) are very similar to those in figure 5. Among the four ERK methods, ERK5(4)5(4) has thelargest number of times of calculating the elements of its Butcher table as shown in panel (b). Two facts account forthis: First, there are 23 di ff erent elements in the Butcher table of ERK5(4)5(4), which are far more than other ERKmethods. Second, as shown in panel (d), attempted trials of ERK5(4)5(4) are more likely to get rejected compare toother ERK methods. Panel (c) shows the number of evaluations of the nonlinear function. Similar to figure 5(c), thetwo 5th-order schemes have least evaluations of the nonlinear function. SS4(3) has the largest number of evaluationsbecause it computes the nonlinear function 24 times in a single step. Panel (e) and (f) show the wall time of integrationand the percentage of time spent on recalculating coe ffi cients. The two-dimensional simulation is about 3- or 4-ordermore expensive than the one-dimensional simulation comparing panel (e) with figure 5(e). In the two-dimensionalintegration, ERK5(4)5(4) spends most of the time recalculating the coe ffi cients. it is less e ffi cient than IF5(4).Similar to the one-dimensional case, explosions in the two-dimensional case can be visualized as homoclinic orbitsof an unstable traveling wave of this system. A two-dimensional traveling wave has form A ( x , y , t ) = A ( x + c x t , y + c y t ) e i ω t , with c x , c y , and ω respectively the translation velocity in the x , y direction, and the phase velocity. By using a shootingmethod [60], we find a traveling wave living in the symmetric subspace with c x = , c y = , ω = . . (23)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 thetwo-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 a ff ected 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 forERK schemes in the comoving frame, which is manifest for the part after the third explosion. Also, for ERK schemesin panel (c), the number of times to evaluate the nonlinear function decreased insignificantly. However, the rejectionrates increased for ERK schemes as shown in panel (d). As a consequence, the total integration time only decreasedby a small amount as shown in panel (e). The improvement of performance in the comoving frame is not striking. Thereason 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 parameterschosen in this paper, the best choice out of the seven schemes is IF5(4).
8. Conclusions
We explored di ff erent 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 di ff erenttime scales. The slow-moving part and the fast-exploding part should use di ff erent step sizes in order to integrate thesystem e ffi ciently. By embedding lower order schemes in IF, ERK and SS methods, local integration error is controllede ff ectively. The step size is adapted to maintain a relative accuracy in each single integration step. For solitons withextremely pulsating amplitudes, time-step adaptation works well in the static frame. While for exploding solitons, tobetter handle the fast-phase rotation di ffi culty, we integrated the system in a comoving frame whose rotating frequency22s similar to that of the soliton solution. We show that integration in the comoving frame can further accelerate theslow-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 thecomoving frame and remarkably slow down the fast-exploding parts. Since ERK5(4)5(4) spends a long time refillingits Butcher table when local error tolerance rtol becomes large, we prefer ERK5(4)5(4) when rtol is small butchoose IF5(4) otherwise. In the two-dimensional case, we find that IF5(4) has the best performance. When the phaserotation is not very severe, a comoving frame may not improve the results much. When one need to implementa stepsize-adaptive scheme quickly, then the 4th-order schemes are suitable. ERK4(3)2(2) may be the best choiceamong all 4th-order methods considered in this paper because of its simple Butcher table structure.The main focus of this paper is to explore di ff erent methods to experiment cubic-quintic complex Ginzburg-Landau equation. Nonetheless, these numerical schemes can be applied to other sti ff (and non-sti ff ) systems thatexhibit intermittent behaviors.
9. Acknowledgements
We are grateful to P. Cvitanovi´c 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 A.9 displays the nonsti ff ReferencesReferences [1] M. C. Cross, P. C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65 (1993) 851–1112. doi:10.1103/RevModPhys.65.851 .[2] I. S. Aranson, L. Kramer, The world of complex Ginzburg-Landau equation, Rev. Mod. Phys. 74 (2002) 99–143. doi:10.1103/RevModPhys.74.99 .[3] N. Akhmediev, J. M. Soto-Crespo, G. Town, Pulsating solitons, chaotic solitons, period doubling, and pulse coexistence in mode-lockedlasers: Complex Ginzburg-Landau equation approach, Phys. Rev. E 63 (2001) 056602. doi:10.1103/PhysRevE.63.056602 .[4] V. L. Ginzburg, Nobel Lecture: On superconductivity and superfluidity (what i have and have not managed to do) as well as on the “physicalminimum” at the beginning of the XXI century, Rev. Mod. Phys. 76 (2004) 981–998. doi:10.1103/RevModPhys.76.981 .[5] W. Chang, A. Ankiewicz, N. Akhmediev, J. M. Soto-Crespo, Creeping solitons in dissipative systems and their bifurcations, Phys. Rev. E 76(2007) 016607. doi:10.1103/PhysRevE.76.016607 .[6] J. M. Soto-Crespo, N. Akhmediev, A. Ankiewicz, Pulsating, creeping, and erupting solitons in dissipative systems, Phys. Rev. Lett. 85 (2000)2937–2940. doi:10.1103/PhysRevLett.85.2937 .[7] J. M. Soto-Crespo, N. Akhmediev, K. S. Chiang, Simultaneous existence of a multiplicity of stable and unstable solitons in dissipativesystems, Phys. Lett. A 291 (2001) 115–123. doi:10.1016/S0375-9601(01)00634-X .[8] C. Cartes, J. Cisternas, O. Descalzi, H. R. Brand, Model of a two-dimensional extended chaotic system: Evidence of di ff using dissipativesolitons, Phys. Rev. Lett. 109 (2012) 178303. doi:10.1103/PhysRevLett.109.178303 .[9] J. Cisternas, O. Descalzi, T. Albers, G. Radons, Anomalous di ff usion of dissipative solitons in the cubic-quintic complex Ginzburg-Landauequation in two spatial dimensions, Phys. Rev. Lett. 116 (2016) 203901. doi:10.1103/PhysRevLett.116.203901 .[10] N. Akhmediev, A. Ankiewicz, Dissipative Solitons, Springer, Berlin, 2005.[11] N. Akhmediev, V. V. Afanasjev, Novel arbitrary-amplitude soliton solutions of the cubic-quintic complex Ginzburg-Landau eequation, Phys.Rev. Lett. 75 (1995) 2320–2323. doi:10.1103/PhysRevLett.75.2320 .[12] N. N. Akhmediev, V. V. Afanasjev, J. M. Soto-Crespo, Singularities and special soliton solutions of the cubic-quintic complex Ginzburg-Landau equation, Phys. Rev. E. 53 (1996) 1190. doi:10.1103/PhysRevE.53.1190 .[13] S. T. Cundi ff , J. M. Soto-Crespo, N. Akhmediev, Experimental evidence for soliton explosions, Phys. Rev. Lett. 88 (2002) 073903. doi:10.1103/PhysRevLett.88.073903 . able A.9: Nonsti ff Tree condition1 (cid:80) β r , α j , r c j = (cid:80) β r , c r = (cid:80) β r , α j , r α k , j c j c k = (cid:80) β r , α j , r c r c j = (cid:80) β r , α j , r α k , j c j = (cid:80) β r , α j , r c r = (cid:80) β r , α j , r α k , j c k = (cid:80) β r , α j , r c j = (cid:80) β r , α j , r c j = (cid:80) β r , c r = (cid:80) β r , α j , r α k , j α p , k c p = (cid:80) β r , α j , r α k , j c k = (cid:80) β r , α j , r α k , j α p , k = (cid:80) β r , α j , r α k , j = Tree condition15 (cid:80) β r , α j , r α k , j c k = (cid:80) β r , α j , r c j = (cid:80) β r , α j , r α k , j = (cid:80) β r , α j , r = (cid:80) β r , α j , r α k , j c k = (cid:80) β r , α j , r c j = (cid:80) β r , α j , r α k , j = (cid:80) β r , α j , r = (cid:80) β r , α j , r c j = (cid:80) β r , c r = (cid:80) β r , α j , r = (cid:80) β r , = Tree condition27 (cid:80) β r , α j , r c j c r = (cid:80) β r , α j , r α k , j c k c r = (cid:80) β r , α j , r α k , j c r = (cid:80) β r , α j , r c j α k , r c k = (cid:80) β r , α j , r α k , r = (cid:80) β r , α j , r c j α k , r = (cid:80) β r , c r α j , r c j = (cid:80) β r , c r α j , r = (cid:80) β r , c r α j , r c j = (cid:80) β r , c r α j , r = (cid:80) β r , c r = [14] J. M. Soto-Crespo, M. Grapinet, P. Grelu, N. Akhmediev, Bifurcations and multiple-period soliton pulsations in a passively mode-locked fiberlaser, Phys. Rev. E 70 (2004) 066612. doi:10.1103/physreve.70.066612 .[15] A. F. J. Runge, N. G. R. Broderick, M. Erkintalo, Observation of soliton explosions in a passively mode-locked fiber laser, Optica 2 (2015)36–39. doi:10.1364/optica.2.000036 .[16] W. Chang, J. M. Soto-Crespo, P. Vouzas, N. Akhmediev, Extreme amplitude spikes in a laser model described by the complex Ginzburg–Landau equation, Opt. Lett. 40 (2015) 2949–2952. doi:10.1364/OL.40.002949 .[17] W. Chang, J. M. Soto-Crespo, P. Vouzas, N. Akhmediev, Extreme soliton pulsations in dissipative systems, Phys. Rev. E 92 (2015) 022926. doi:10.1103/PhysRevE.92.022926 .[18] E. N. Tsoy, N. Akhmediev, Bifurcations from stationary to pulsating solitons in the cubic-quintic complex Ginzburg-Landau equation, Phys.Lett. A 343 (2005) 417–422. doi:10.1016/j.physleta.2005.05.102 .[19] S. C. Mancas, S. R. Choudhury, A novel variational approach to pulsating solitons in the cubic-quintic Ginzburg–Landau equation, Theoret.Math. Phys. 152 (2007) 1160–1172. doi:10.1007/s11232-007-0099-8 .[20] N. Akhmediev, J. M. Soto-Crespo, Strongly asymmetric soliton explosions, Phys. Rev. E 70. doi:10.1103/PhysRevE.70.036613 .[21] S. Wang, L. Zhang, An e ffi cient split-step compact finite di ff erence method for cubic–quintic complex Ginzburg–Landau equations, Comput.Phys. Commun. 184 (2013) 1511–1521. doi:10.1016/j.cpc.2013.01.019 .[22] A. Taleei, M. Dehghan, Time-splitting pseudo-spectral domain decomposition method for the soliton solutions of the one- and multi-dimensional nonlinear Schr¨odinger equations, Comput. Phys. Commun. 185 (2014) 1515–1528. doi:10.1016/j.cpc.2014.01.013 .[23] H. Wang, An e ffi cient Chebyshev–Tau spectral method for Ginzburg–Landau–Schr¨odinger equations, Comput. Phys. Commum. 181 (2010)325–340. doi:10.1016/j.cpc.2009.10.007 .[24] M. Hochbruck, A. Ostermann, Exponential integrators, Acta Numerica 19 (2010) 209–286. doi:10.1017/S0962492910000048 .
25] J. D. Lawson, Generalized Runge-Kutta processes for stable systems with large Lipschitz constants, SIAM J. Numer. Anal. 4 (1967) 372–380. doi:10.1137/0704033 .[26] H. Montanelli, N. Bootland, Solving sti ff PDEs in 1D, 2D and 3D with exponential integrators (2016).URL http://arxiv.org/abs/1604.08900 [27] B. M. Caradoc-Davies, Vortex dynamics in Bose-Einstein condensates, Ph.D. thesis, Univ. Otago, Dunedin, New Zealand (2000).[28] J. Hult, A fourth-order Runge-Kutta in the interaction picture method for simulating supercontinuum generation in optical fibers, J. LightwaveTechnol. 25 (2007) 3770–3775. doi:10.1109/JLT.2007.909373 .[29] S. Krogstad, Generalized integrating factor methods for sti ff PDEs, J. Comput. Phys. 203 (2005) 72–88. doi:10.1016/j.jcp.2004.08.006 .[30] A. Taflove, S. C. Hagness, Computational Electrodynamics: The Finite-di ff erence Time-domain Method, Artech House, 2005.[31] R. Holland, Finite-di ff erence time-domain (FDTD) analysis of magnetic di ff usion, IEEE Trans. Electromagn. Compat. 36 (1994) 32–39. doi:10.1109/15.265477 .[32] P. G. Petropoulos, Analysis of exponential time-di ff erencing for FDTD in lossy dielectrics, IEEE Trans. Antennas Propag. 45 (1997) 1054–1057. doi:10.1109/8.585755 .[33] C. Schuster, A. Christ, W. Fichtner, Review of FDTD time-stepping schemes for e ffi cient simulation of electric conductive media, Microw.Opt. Technol. Lett. 25 (2000) 16–21. doi:10.1002/(SICI)1098-2760(20000405)25:1<16::AID-MOP6>3.0.CO;2-O .[34] S. M. Cox, P. C. Matthews, Exponential time di ff erencing for sti ff systems, J. Comput. Phys. 176 (2002) 430–455. doi:10.1006/jcph.2002.6995 .[35] A.-K. Kassam, L. N. Trefethen, Fourth-order time-stepping for sti ff PDEs, SIAM J. Sci. Comput. 26 (2005) 1214–1233. doi:10.1137/S1064827502410633 .[36] A. Friedli, Verallgemeinerte Runge-Kutta Verfahren zur L¨oesung steifer Di ff erentialgleichungssysteme, in: R. Bulirsch, R. D. Grigorie ff ,J. Schr¨oder (Eds.), Numerical Treatment of Di ff erential Equations: Proceedings of a Conference, Held at Oberwolfach, July 4–10, 1976,Springer, Berlin, 1978, pp. 35–50. doi:10.1007/BFb0067462 .[37] H. Berland, B. Owren, B. Skaflestad, B-series and order conditions for exponential integrators, SIAM J. Numer. Anal. 43 (2005) 1715–1727. doi:10.1137/040612683 .[38] M. Hochbruck, A. Ostermann, Explicit exponential Runge–Kutta methods for semilinear parabolic problems, SIAM J. Numer. Anal. 43(2005) 1069–1090. doi:10.1137/040611434 .[39] V. T. Luan, A. Ostermann, Sti ff order conditions for exponential Runge–Kutta methods of order five, in: G. H. Bock, P. X. Hoang, R. Ran-nacher, P. J. Schl¨oder (Eds.), Modeling, Simulation and Optimization of Complex Processes - HPSC 2012: Proceedings of the Fifth Inter-national Conference on High Performance Scientific Computing, March 5-9, 2012, Hanoi, Vietnam , Springer, Cham, 2014, pp. 133–143. doi:10.1007/978-3-319-09063-4_11 .[40] V. T. Luan, A. Ostermann, Exponential B-series: The sti ff case, SIAM J. Numer. Anal. 51 (2013) 3431–3445i. doi:10.1137/130920204 .[41] K. A. Bagrinovskii, S. K. Godunov, Di ff erence schemes for multi-dimensional problems, Dokl. Acad. Nauk. 115 (1957) 431–433.[42] G. Strang, On the construction and comparison of di ff erence schemes, SIAM J. Numer. Anal. 5 (1968) 506–517. doi:10.1137/0705041 .[43] H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A 150 (1990) 262–268. doi:10.1016/0375-9601(90)90092-3 .[44] W. Auzinger, W. Herfort, Local error structures and order conditions in terms of Lie elements for exponential splitting schemes, OpusculaMath. 34 (2014) 43–255. doi:10.7494/opmath.2014.34.2.243 .[45] W. Auzinger, H. Hofst¨atter, D. Ketcheson, O. Koch, Practical splitting methods for the adaptive integration of nonlinear evolution equations.Part I: construction of optimized schemes and pairs of schemes, BIT Numer. Math. 57 (2016) 55–74. doi:10.1007/s10543-016-0626-9 .[46] S. Blanes, F. Casas, A. Farr´es, J. Laskar, J. Makazaga, A. Murua, New families of symplectic splitting methods for numerical integration indynamical astronomy, Appl. Numer. Math. 68 (2013) 58 – 72. doi:10.1016/j.apnum.2013.01.003 .[47] T. Mohammedi, A. Aissat, An accurate Fourier splitting scheme for solving the cubic-quintic complex Ginzburg–Landau equation, Superlat-tices Microstruct. 75 (2014) 424–434. doi:10.1016/j.spmi.2014.08.007 .[48] F. B´erard, C.-J. Vandamme, S. C. Mancas, Two-dimensional structures in the quintic Ginzburg–Landau equation, Nonlinear Dyn. 81 (2015)1413–1433. doi:10.1007/s11071-015-2077-2 .[49] A.-K. Kassam, L. N. Trefethen, Solving reaction-di ff usion equations 10 times faster, Numerical Analysis Group Research Report 1192,Oxford University (2003).[50] P. Whalen, M. Brio, J. V. Moloney, Exponential time-di ff erencing with embedded Runge-Kutta adaptive step control, J. Comput. Phys. 280(2015) 579–601. doi:10.1016/j.jcp.2014.09.038 .[51] O. Koch, C. Neuhauser, M. Thalhammer, Embedded exponential operator splitting methods for the time integration of nonlinear evolutionequations, Appl. Numer. Math. 63 (2013) 14–24. doi:10.1016/j.apnum.2012.09.002 .[52] W. Auzinger, O. Koch, , M. Quell, Adaptive high-order splitting methods for systems of nonlinear evolution equations with periodic boundaryconditions, Numer. Algorithms (2016) 1–23 doi:10.1007/s11075-016-0206-8 .[53] J. R. Dormand, P. J. Prince, New Runge-Kutta algorithms for numerical simulation in dynamical astronomy, Celest. Mech. 18 (1978) 223–232. doi:10.1007/BF01230162 .[54] J. R. Dormand, P. J. Prince, A family of embedded Runge-Kutta formulae, J. Comput. Appl. Math. 15 (1986) 203–211. doi:10.1016/0771-050X(80)90013-3 .[55] S. Balac, F. Mah´e, Embedded Runge-Kutta scheme for step-size control in the interaction picture method, Comput. Phys. Commun. 184(2013) 1211–1219. doi:10.1016/j.cpc.2012.12.020 .[56] V. T. Luan, A. Ostermann, Explicit exponential Runge–Kutta methods of high order for parabolic problems, J. Comput. Appl. Math. 256(2014) 168–179. doi:10.1016/j.cam.2013.07.027 .[57] X. Ding, P. Cvitanovi´c, Exploding relative periodic orbits in cubic-quintic complex Ginzburg–Landau equation, in preparation (2017).[58] K. Levenberg, A method for the solution of certain non-linear problems in least squares, Quart. Appl. Math. 2 (1944) 164–168. doi:10.1090/qam/10666 .
59] D. M. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, J. Soc. Indust. Appl. Math. 11 (1963) 431–441. doi:10.1137/0111030 .[60] G. J. Chandler, R. R. Kerswell, Invariant recurrent solutions embedded in a turbulent two-dimensional Kolmogorov flow, J. Fluid Mech. 722(2013) 554–595. doi:10.1017/jfm.2013.122 .[61] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Di ff erential Equations, 2ndEdition, Springer, Berlin, 2006. doi:10.1007/3-540-30666-8 .[62] S. Balac, High order embedded Runge-Kutta scheme for adaptive step-size control in the interaction picture method, J. KSIAM 17 (2013)238–266. doi:10.12941/jksiam.2013.17.238 .[63] M. P. Calvo, A. M. Portillo, Variable step implementation of ETD methods for semilinear problems, Appl. Math. Comput. 196 (2008)627–637. doi:10.1016/j.amc.2007.06.025 .[64] O. Descalzi, H. R. Brand, Transition from modulated to exploding dissipative solitons: Hysteresis, dynamics, and analytic aspects, Phys. Rev.E 82 (2010) 026203. doi:10.1103/PhysRevE.82.026203 .[65] J. Cisternas, O. Descalzi, C. Cartes, The transition to explosive solitons and the destruction of invariant tori, Cent. Eur. J. Phys. 10 (2012)660–668. doi:10.2478/s11534-012-0023-1 .[66] J. Cisternas, O. Descalzi, Intermittent explosions of dissipative solitons and noise-induced crisis, Phys. Rev. E 88 (2013) 022903. doi:10.1103/PhysRevE.88.022903 .[67] J. M. Ghidaglia, B. H´eron, Dimension of the attractors associated to the Ginzburg-Landau partial di ff erential equation, Physica D 28 (1987)282–304. doi:10.1016/0167-2789(87)90020-0 ..