{'fraction_non_alphanumeric': 0.04365135122476699, 'fraction_numerical': 0.017188890450096014, 'mean_word_length': 4.442773387713484, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-Order Schemes of Exponential Time Differencing for Stiff Systems with Nondiagonal Linear Part

Exponential time differencing methods is a power toll for high-performance numerical simulation of computationally challenging problems in condensed matter and chemical physics, where mathematical models often possess fast oscillating or decaying modes -- in other words, are stiff systems. Practical implementation of these methods for the systems with nondiagonal linear part of equations is exacerbated by infeasibility of an analytical calculation of the exponential of a nondiagonal linear operator; in this case, the coefficients of the exponential time differencing scheme cannot be calculated analytically. We suggest an approach, where these coefficients are numerically calculated with auxiliary problems. We rewrite the high-order Runge--Kutta type schemes in terms of the solutions to these auxiliary problems and practically examine the accuracy and computational performance of these methods for a heterogeneous Cahn--Hilliard equation and a sixth-order spatial derivative equation governing pattern formation in the presence of an additional conservation law.

In condensed matter and chemical physics, many systems are governed by the mathematical models which can be called 'stiff' systems-the systems where some modes are fast oscillating or decaying. The dynamics of wave function in atom lattices with frozen parametric disorder in potential is fundamentally contributed by the fast oscillating modes and this dynamics gives rise to the phenomenon of Anderson localization [1][2][3][4][5][6][7] . The spinodal decomposition [8][9][10][11][12][13] of two-component mixtures and heat-mass-transfer in active media [14][15][16][17][18][19][20][21] is governed by partial differential equations with high-order spatial derivatives; these derivatives drive a fast decay of short wavelength perturbations. The systems near a SNIPER (saddlenode infinite period) bifurcation of chemical oscillations 22 are subject to thermal noise and can be subject to other noise sources; at the bifurcation point, these systems can be described by the same mathematical models as the quadratic integrate-and-fire neurons with noise and synaptic coupling 23 . The numerical simulation of the Fokker-Planck equation of these chemical oscillation systems is again the case of a stiff system. Generally, fast oscillating or decaying modes impose severe limitation on the time stepsize, which has to be very small. For oscillating modes it is needed to maintain a demanded accuracy of simulation 7 . Fast decaying modes are dying-out and do not influence the dynamics of the chemical/physical system; therefore, there is no need for a high precision of their simulation. However, for a finite stepsize, fast decay can produce numerical instability if the stepsize in not small enough. As a result, the numerical simulation of stiff system requires the employment of ad hoc approaches and implicit schemes, which require a sophisticated individual mathematical preliminary work, is costly or infeasible in higher dimensions, and are helpless for fast oscillatory modes in conservative systems. Otherwise, the simulation of these problems should be performed with a tiny time stepsize and becomes extremely CPU (central processing unit) time consuming. An alternative approach is a relatively new class of methods of exponential time differencing (ETD) [24][25][26][27][28] .
In ETD methods 24 the linear part of equations, which gives the origin to the fast modes, is solved exactly. For the nonlinear part, one constructs an exact solution for a quasistatic approximation of these terms during one time step in the case of a first-order scheme. For an nth-order scheme, one effectively constructs the exact solution for the case where the nonlinear part in one time step is an (n − 1)-order polynomial of time 24 . This approach provides high resolution for temporal dynamics of fast modes and prevents numerical instabilities even for a quite large values of time stepsize.
The application of ETD methods is straightforward when the linear part is diagonal, which is abundant for spectral methods for the problems with periodic boundary conditions and Fourier basis. However, for spectral methods with other basis functions (e.g., Chebyshev polynomials) the linear part can become nondiagonal. Moreover, for problems with frozen heterogeneity of parameters 23,[29][30][31][32][33] , the linear part becomes essentially nondiagonal. In this case the analytical calculation of the coefficients of ETD schemes (considered in Sec. I) can become not just difficult but completely impossible.
In this paper we develop an approach to calculation of the coefficients of high-order schemes of exponential time differencing for systems with a nondiagonal linear part. In this approach the coefficients are calculated by means of numerical integration of auxiliary problems. This integration is carriedout with a plain standard method, like the predictor-corrector one, and a very small time stepsize τ 1 , ensuring numerical stability, but over a short time interval-the time step τ of the ETD scheme. The basic concept was suggested in [34] and here we develop and apply the approach to the Runge-Kutta type schemes 24 of 3rd and 4th order. For the sake of illustration, we consider the numerical simulation of the discretespace versions of partial differential equations: Cahn-Hilliard equation 8,9 and Matthews-Cox equation 16,17 . The performance gain varies from two to many orders of magnitude. In multiple dimensions 34 , the ETD methods yield significant gain in performance, but in some cases an extensive usage of memory for storage of scheme coefficients is required. An extensive usage of memory (in some cases) and the fact, that the bulk of the simulation process is linear algebra operations, make our approach reasonably suitable for parallel-and supercomputing.
The paper is organized as follows. In Sec. I, we provide three exponential differencing schemes of the Runge-Kutta type and rewrite them in terms of the solutions of the auxiliary problems. The procedure for constructing the solutions of the auxiliary problems is provided. In Sec. II, we considers examples of the application of the ETD schemes to two paradigmatic mathematical models: Cahn-Hilliard equation (Sec. II A) and sixth-order spatial derivative Matthews-Cox equation (Sec. II B). The accuracy and CPU time performance of the ETD schemes is analysed. Sec. III finalizes the paper with conclusion.

I. EXPONENTIAL TIME DIFFERENCING
We deal with the problem of numerical integration of the equations system of the following form: where u(t) is an N-component vector, L is an [N × N]-matrix with time-independent elements, f(u,t) is the nonlinear part of equations. The decomposition of the equations into the linear part with time-independent coefficients, L · u, and the nonlinear part is not unique and guided by the merit of convenience: f(u,t) can even contain a time-independent linear part. However, it is desirable that all the dominating sources of numerical instability should be collected into L · u as fully as possible.
In this paper we consider implementation of three exponential differencing schemes of the Runge-Kutta type 24 ; for equation system (1) these schemes take the following form: (ETD2RK) two-step scheme: where τ is the time stepsize of the numeric scheme, a is a preliminary approximation of u for the next time step t + τ, I is the unitary matrix, the exponential of a matrix is defined by the series exp A = I + A + A 2 /2! + A 3 /3! + . . .; (ETD3RK) three-step scheme: where a and b are the preliminary approximations of u for t + τ/2 and t + τ, respectively; (ETD4RK) four-step scheme: The error of the ETD2RK scheme on one time step is −τ 3f /12 [24], the ETD3RK scheme error is ∝ τ 4 d 3 f/dt 3 , the ETD4RK scheme error is ∝ τ 5 d 4 f/dt 4 . In the numeric schemes above, the reciprocal matrix L −1 is introduced for the brevity of expressions using the exponential; wherever the matrix L −n appears in equations, it is multiplied by expressions the series representation of which starts from the matrix L m , m ≥ n. For instance, L −1 · (e Lτ − I) is actually a formal expression for the series τI + τ 2 L/2 + τ 3 L 2 /3! + τ 4 L 3 /4! + . . . . Therefore, the zero eigenvalues of the matrix L, which would make L −1 → ∞, are not an issue; the divergence of L −1 only prevents the usage of a shorter form of equations. Albeit, in this paper, we are handling the case of a nondiagonal shape of matrix L, for which analytical calculation of the matrix e Lτ can be impossible or problematic.
More convenient for our approach is to recast the numeric schemes in terms of the solutions of the following auxiliary problems 34 . The general solution to the Cauchy probleṁ is given by the matrix Q(τ) ≡ e Lτ such that Hence, one obtains the definition of the matrix Q(τ), which is convenient for numerical simulations: The solutions to the problemṡ are given by the matrices M n (τ) ≡ τ 0 e L(τ−t) t n−1 dt : Hence, one obtains the definition of the matrix M n (τ), which is convenient for numerical simulations: Using the recurrent relationship M n+1 (τ) = L −1 · (−τ n I + nM n (τ)), one can obtain As proposed in [34], for an essentially nondiagonal shape of matrix L, the evaluation of matrices Q and M n defined by equations (6) and (8) can be conducted via the direct numerical integration of problems (5) and (7) with a very small time stepsize but on a short time interval -one step of the ETD scheme τ.
In terms of matrices M n the schemes (2), (3), and (4) acquire the following form: (ETD2RK): where Q τ , and for the full stepsize τ the corresponding subscript is omitted; For a broad class of active media systems the pattern formation in thin layers is governed by the Cahn-Hilliard equation 14,15,18 (CHE), which as well describes the spinodal decomposition of two-component mixtures 8,9,13 . For the generality of consideration, we will also allow for an advective transfer along the layer 9,18 . In the one-dimensional case (the patterns are homogeneous along the second direction in the layer plane), CHE for the field u(x,t) with advection velocity v reads where q(x) is the local deviation of the bifurcation parameter from the instability threshold of the homogeneous infinite layer (the patterns are excited for q > 0).
For an example, we consider CHE in a domain 0 < x < L with trivial boundary conditions: The local excitability parameter q(x) is assumed to be positive within certain excitation zone and negative beyond it: For a localized excitation q(x), strong enough advection v results in an oscillatory behavior (e.g., see [33]). For v = 1 and q(x) given by (13), the oscillatory solution of Eq. (12) is presented in Fig. 1.

Numeric scheme
The spatially-discrete version of CHE (12) for vector u(t) ≡ {u j (t)| j = 1, 2, ..., N}, where u j (t) = u( jh x ,t), h x = L/N, can be written as Eq. (1) with where q j = q( jh x ) (13). The trivial boundary conditions we consider imply that one can formally set u −1 = u 0 = u N+1 = u N+2 = 0 in the right-hand parts of (14). The problem (1), (14) can be numerically simulated with the basic Euler scheme or the predictor-corrector (PC) one. For both of these schemes, the forth-order x-derivative term imposes the limitation on the time stepsize, τ 1 < h 4 x /8; for a bigger time stepsize the direct numerical simulation becomes unstable. More sophisticated Runge-Kutta type schemes or any other high-order methods are not needed here as the PC scheme with such a small time stepsize already warrants an excessive accuracy of the numerical integration in time. Henceforth, we will employ the PC scheme: For problem (1), (14), the PC scheme (15) was employed to integrate the auxiliary problems (5) and (7) for t ∈ [0, τ] to evaluate matrices Q τ 2 , Q (6) and M 1, τ 2 , M n , n = 1, 2, 3 (8), respectively. The time stepsize for the direct numerical integration of the auxiliary problems was τ 1 = 0.1h 4 x . In Fig. 2, a sample structure of these matrices can be seen.

Accuracy and performance of numerical simulation with exponential time differencing schemes
The ETD2RK, ETD3RK, and ETD4RK schemes with matrices Q τ 2 , Q, M 1, τ 2 , M n , n = 1, 2, 3 calculated as described in Sec. II A 1 were employed for the numerical simulation of CHE (12) with v = 1 and q(x) given by (13) in the domain of length L = 10 (see Fig. 1). In Figs. 3 and 4, the dependencies of the accuracy and the simulation performance versus the stepsize τ of the ETD schemes are presented.
In Fig. 3, one can compare the performance of ETD methods to that of the PC method (blue solid circles). The numeric error of a specific scheme is evaluated as the deviation of the result from the result calculated with the PC scheme (15) with τ 1 = 0.01h 4 x , which is by a factor of 10 smaller then the time stepsize 0.1h 4 x needed for practical computations. The accuracy of the PC method is excessive, but it is dictated by the maximal admissible value of the time stepsize h 4 x /8. In turn, the ETD methods can provide decent accuracy even for large values of τ, where they perform much faster (see the right panel) than the PC method. For example, with a required accuracy of 10 −6 (see crests in Fig. 3), one can pick-up as large time stepsize as τ ≈ 0.005 for the ETD2RK scheme, 0.02 for ETD3RK, and 0.04 for ETD4RK. The biggest performance gain is achieved with ETD4RK scheme and exceeds a factor of 100 for all three schemes. For an increased accuracy of the spatial discretization (see Fig. 4 where N is doubled compared to the case of Fig. 3), the performance gain provided by ETD schemes significantly increases.
Note, the CPU time for the preliminary calculations of Q and M n increases linearly with τ and can become nonsmall (see open symbols in the right panels of Figs. 3 and 4). This CPU time can be significant for the problems where simulation over short or moderate intervals of time t is sufficient. For these problems, an optimal performance is achieved for the values of τ, which provide approximately equal CPU times for the preliminary calculation of Q and M n and the run of the ETD scheme (see the crossings of the lines marked by solid and open symbols in Figs. 3 and 4). The overall performance gain is decreased, but still large. The ETD2RK scheme becomes the most efficient one in this case. However, the most computationally demanding problems are those of complex dynamics [9][10][11][19][20][21][35][36][37] (including spatiotemporal chaos) and dynamics in the presence of frozen parametric disorder 29-33 , where the CPU time for the preliminary calculations can be neglected.

B. 1D pattern formation with a conservation law
The conservation laws (of the chemical specie mass, etc.) are structurally stable (persistent) properties of systems and influence the general form of equations governing the pattern formation. Matthews and Cox 16 studied the model equation governing pattern formation with an additional conservation law: where the terms inside the brackets are just those of the wellstudied Swift-Hohenberg equation 38 and an additional advective v-term is introduced 17  ary conditions: The local excitability parameter q(x) is assumed to be the same as for CHE (13). For a localized excitation q(x), strong enough advection v results in an oscillatory behavior. For v = 1 and q(x) given by (13), the oscillatory solution of Eq. (16) is presented in Fig. 1.

Numeric scheme
The spatially-discrete version of MCE (16) for vector u(t) ≡ {u j (t)| j = 1, 2, ..., N}, where u j (t) = u( jh x ,t), h x = L/N, can be written as Eq. (1) with where q j = q( jh x ) (13). The trivial boundary conditions we consider imply that one can formally set u −2 = u −1 = u 0 = u N+1 = u N+2 = u N+2 = 0 in the right-hand parts of (17). The problem (1), (17) can be numerically simulated with the basic Euler scheme or the predictor-corrector one. For both of these schemes, the sixth-order x-derivative term imposes the limitation on the time stepsize, τ 1 < h 6 x /32; for a bigger time stepsize the direct numerical simulation becomes unstable. Henceforth, we will employ the PC scheme (15).

Accuracy and performance of numerical simulation with exponential time differencing schemes
The ETD2RK, ETD3RK, and ETD4RK schemes with matrices Q τ 2 , Q, M 1, τ 2 , M n , n = 1, 2, 3 calculated as described in Sec. II B 1 were employed for the numerical simulation of MCE (16) with v = 1 and q(x) given by (13) in the domain of length L = 10 (see Fig. 1). In Figs. 6 and 7, the dependencies of the accuracy and the simulation performance versus the stepsize τ of the ETD schemes are presented. One can see even more dramatic performance gain than for the numerical simulation of CHE (12).
Generally, for the partial differential equations with the highest-order spatial derivative ∂ m /∂ x m , the performance gain for the ETD methods is ∝ h m−2 x , as anticipated in [34]. However, in Figs. 2 and 5, one can see that matrix elements rapidly decay away from the diagonal. For nonlarge τ, the absolute values of a significant fraction of elements are below ∼ 10 −16 which is the level of the double precision computer accuracy. The elements smaller than the required error level can be set to zero without any damage to the scheme accuracy. The optimization of a program code for the sparseness of matrices Q and M n allows for a noticeable performance acceleration. For a given stepsize τ, the highest derivative ∂ m /∂ x m creates a bigger spread of nonsmall values of matrix elements (compare Fig. 2 to Fig. 5, where the stepsize τ is even somewhat smaller). Therefore, this additional optimization becomes less beneficial for high-order derivatives, where the gain is actually large without any optimization.

III. CONCLUSION
We have suggested a practical approach for the implementation of the exponential time differencing schemes of the Runge-Kutta type 24 for numerical simulation of stiff systems with nondiagonal linear part. In this approach, the numeric schemes are written in terms of the solutions of auxiliary problems which are preliminary integrated numerically by means of an explicit basic method with a very small time step but over a short interval of time-one step of the ETD scheme. As a result, one can not only use the ETD methods for the systems where an analytical calculation of the exponential of a nondiagonal linear operator is infeasible, but also employ the same program code for different equation systems, the only part of code which is subject to change is the subroutine with the equations to be simulated.
The employment of this approach for Cahn-Hilliard equation (12) and sixth-order spatial derivative Matthews-Cox equation (16) demonstrated that one can have a performance gain by two orders of magnitude for the former and by three orders for the latter. Moreover, for partial differential equations, the performance gain grows as one increases the re-quired accuracy of spatial resolution; the basic gain for an equation with the highest-order spatial derivative ∂ m /∂ x m is ∝ h m−2 x . For nonlarge m, one can introduce code optimization accounting for the sparseness of the matrices Q and M n , which gives some further performance acceleration on top of the basic one (this acceleration is by far not that large as the basic gain). For the simulation of a long-term evolution, as required for complex spatiotemporal dynamics or Anderson localization phenomena, the maximal performance gain is achieved with the 4th-order Runge-Kutta type ETD scheme; for a short-time simulation, the usage of the 2th-order ETD scheme is somewhat more efficient.
The suggested approach can be extended into multiple dimensions in a straightforward way 34 . However, the basic version of the code for a partial differential equation, say, in 2d requires matrices Q and M n of size [N × N] with N = N x × N y , where N x and N y are the number of nodes (or modes) in the x-and y-directions. With the optimization for the sparseness of Q and M n , one can store a reduced amount [N × N spr ] of elements of each matrix, where N spr = 2mτ min N y |c x | ln[h x /(err τ)], N x |c y | ln[h y /(err τ)] for the highest-order spatial derivative term c x ∂ m u/∂ x m + c y ∂ m u/∂ y m and the acceptable relative error per a unit time err. For high-precision simulations in multiple dimensions the usage of the ETD schemes becomes memory-demanding. Elevated memory requirements and the fact, that the bulk of the computation process is the multiplication of large matrices and vectors, make these high-performance simulation methods naturally suitable for parallel-and super-computing. ACKNOWLEDGMENTS We thank A. Pikovsky for drawing our attention to the ETD methods when we encountered the problems a conclusive numerical study of which with conventional methods was practically impossible; we appreciate his comments on the presented approach.

Conflict of Interest
The authors have no conflicts to disclose.

DATA AVAILABILITY STATEMENT
The data that supports the findings of this study are available within the article in the graphic form. The data sheets for the graphs and the program codes in FORTRAN are available on request from the authors.