{'fraction_non_alphanumeric': 0.04463737253777796, 'fraction_numerical': 0.01429214955567386, 'mean_word_length': 4.226883561643835, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-order splitting methods for separable non-autonomous parabolic equations

We consider the numerical integration of non-autonomous separable parabolic equations using high order splitting methods with complex coefficients (methods with real coefficients of order greater than two necessarily have negative coefficients). We propose to consider a class of methods in which one set of the coefficients are real and positive numbers, and to split properly the system in the extended phase space where the time is taken as a new coordinate. This allows us to evaluate all time-dependent operators at real values of the time, leading to schemes which are stable and simple to implement. If the system can be considered as the perturbation of an exactly solvable problem and the flow of the dominant part is advanced using the real coefficients, it is possible to build highly efficient methods for these problems. We show the performance of this class of methods on several numerical examples and present some new improved schemes.


Introduction
We consider the numerical integration of non-autonomous separable parabolic equations using high order splitting methods with complex coefficients. This class of methods has been recently used for the numerical integration of the autonomous case, showing good performances [5,12,17]. Splitting methods with real coefficients of order greater than two necessarily have negative coefficients and can not be used for solving these problems [4,15,21,23]. However, solutions with complex coefficients with positive real part exist, and some of these methods can provide a high performance in spite the equations have to be solved in the complex domain. Previous works with applications among other in celestial mechanics and quantum mechanics where splitting methods with complex coefficients are considered already exist [2,3,13,19,20,22,23,24].
A straightforward application of splitting methods with complex coefficients to non-autonomous problems require the evaluation of the time-dependent functions on the operators at complex times, and the corresponding flows in the numerical scheme are, in general, not well conditioned. In this work we propose to consider a class of splitting methods in which one set of the coefficients belong to the class of real and positive numbers. This can allow to evaluate all time-dependent operators at real values of the time, leading to schemes which are stable and simple to implement.
If the system can be considered as the perturbation of an exactly solvable problem (or easy to numerically solve) and the flow of the dominant part is advanced using the real coefficients, it is possible to build highly efficient methods for these problems.

The problem
Let us consider the non-autonomous separable PDE u(x, t) ∈ R D , and where the (possibly unbounded) operators A, B and A + B generate C 0 semi-groups for positive t over a finite or infinite Banach space. Equations of this form are encountered in the context of parabolic partial differential equations, an example being the inhomogeneous non-autonomous heat equation where t ≥ 0, x ∈ R d or x ∈ T d and ∆ denotes the Laplacian with respect to the spatial coordinates, x. Another example corresponds to reaction-diffusion equations of the form ∂u ∂t = D(t)∆u + B(t, u), (1.3) where D(t) is a matrix of diffusion coefficients (typically a diagonal matrix) and B accounts for the reaction part. In general, A(t, u), B(t, u) can also depend on x, ∇, etc., which are omitted for clarity in the presentation.
For simplicity, we write the non-linear equation (1.1) in the (apparently) linear form where L A , L B are the Lie operators associated to A, B, i.e.
which act on functions of u If the problem is autonomous, the formal solution is given by u(t) = e t(L A(u) +L B(u) ) u 0 , which is a short way to write If the subproblems have exact solutions or can efficiently be numerically solved, it is usual to consider splitting methods as numerical integrators. If we denote by e hL A(u) , e hL B(u) the exact h-flows for each problem in (1.5) (and for a sufficiently small time step, h) the simplest method within this class is the Lie-Trotter splitting which is a first order approximation in the time step to the solution, while the symmetrized version is referred to as Strang splitting, and is an approximation of order 2, i.e. S(h) = e h(L A(u)+B(u) ) +O(h 3 ). Upon using an appropriate sequence of steps, high-order approximations can be obtained as Ψ(h) = e hb m+1 L B e hamL A · · · e hb 2 L B e ha 1 L A e hb 1 L B , (1.8) and methods with real coefficients at any order can be obtained [14,22,26]. However, as already mentioned, splitting methods of order greater than two (with real coefficients) have at least one of the coefficients a i negative as well as at least one of the coefficients b i so, the flows e tL A and/or e tL B may not be well defined (this is indeed the case, for instance, for the Laplacian operator) and this prevents the use of methods which embed negative coefficients. For this reason, exponential splitting methods of at most order p = 2 have been considered up to recently.
In order to circumvent this order-barrier, the papers [12] and [17] simultaneously presented a systematic analysis for a class of composition methods with complex coefficients having positive real parts. Using this extension from the real line to the complex plane, the authors of [12] and [17] built up methods of orders 3 to 14 by considering a technique known as triple-jump composition. More efficient high order methods are obtained in [5].
In this work we are interested, however, in the numerical integration of the non-autonomous problem (1.1) where the use of complex coefficients involve additional constraints as we will see. A method of choice for solving numerically (1.1) consists in advancing the solution alternatively along the exact (or numerical) solutions of the two problems The exact flows are, in general, not known. This is the case, for example, if = 0 (and similarly for B(t)).
If the exact solution is not known, it can be replaced by a sufficiently accurate numerical approximation.
This procedure is equivalent to take the time as two new coordinates, t 1 , t 2 with ′ ≡ d dt , and to split the system in the extended space as follows [8] (1.11) A more convenient way to split the system which transforms the non-autonomous problems into autonomous is the following (1.12) Notice that the explicit time-dependency in A and B is frozen in each subproblem and the formal solution corresponds to the exponential of the Lie operators where the time dependency in the operators are frozen on each time interval. Unfortunately, to use splitting methods with complex coefficients for non-autonomous problems requires, in general, to compute A(t, u), B(t, u) for t ∈ C, leading, in general, to badly conditioned algorithms.
In this work we show that splitting method having one set of coefficients real and positive valued, i.e.
allow to build algorithms where the operators A(t, u), B(t, u) are evaluated only for t ∈ R, leading to well defined methods. Several splitting methods with this structure have already been constructed 1 .
We will also explore the case in which B ≪ A , which we refer as a perturbed problem. We first show how this class of methods has to be used in these problems and next we study how to build high order efficient methods for these problems.

Splitting methods for non-autonomous problems
Suppose we have a splitting method with say, a i ∈ R + and b i ∈ C + . To solve the eq. (1.9) we propose to take the time as one new coordinate and split the system as follows [8] , a fourth-order method was obtained with a i ∈ R + . In a similar way, in [5] sixth-order schemes were also explored with a i ∈ R + . The coefficients can be found at: http://www.gicas.uji.es/Research/splitting-complex.html.

Let us denote by Φ [a i h] A
the map associated to the exact solution (or a sufficiently accurate numerical approximation) of the non-autonomous equation and a 0 = 0. Then, the splitting method (1.8) for the non-autonomous equation reads now . Notice that in this scheme, since t 1 is advanced with the coefficients a i (and then it takes real values) the operators A(t, u) and B(t, u) are evaluated on real values of t. On the other hand, if A(t, u) is an unbounded operator and a numerical methods is used to approximate the flow Φ

[h]
A , it must be well defined for 0 ≤ h < h * for some positive h * , and this is not guaranteed for general methods. For example, some commutator-free methods up to fourth-order can be used. Given the equation A = e hL A(tn+h/2,u) corresponds to a symmetric second order method, and h, u , corresponds to a fourth-order method [10,25]. Notice that the Lie operators, since being derivatives, are written in the reverse order than the maps, and this is very important to keep inmind for non-linear non-autonomous problems in order to apply the method correctly (see [9] for more details on the Magnus series expansion and Magnus integrators for nonautonomous non-linear differential equations). This scheme corresponds to the composition of the 1-flow maps for the equations and the solution given by the map corresponds to u 2 (1). These second and fourth-order commutator-free methods can be used for unbounded operators and higher order commutator-free Magnus integrators for unbounded operators are under investigation at this moment [7].

Splitting methods for non-autonomous perturbed systems
In some cases, the system can be considered as the perturbation of an exactly solvable problem. In those cases, it is usually convenient to split into the dominant part and the perturbation and to build methods which take advantage of this relevant property. However, if the problem is non-autonomous and the time-dependency is not treated properly, the performance of the methods designed for perturbed problems deteriorate considerably.
Suppose that B(t, u) ≪ A(t, u) . To make this fact more evident, we replace B by εB with |ε| ≪ 1 2 . In the autonomous case, for example, the Lie-Trotter composition for this split satisfies i.e., it has a local error of order O(εh 2 ).
Since A and B are qualitatively different for perturbed problems, it is usual to consider ABA and BAB compositions.
An m-stage symmetric BAB compositions given by . ., and ABA compositions are given by . .. We will use the following short notation for these methods Notice that eq. (2.6) is a BAB composition which, for the particular case where b 1 = b m+1 = 0 transforms into a ABA composition (but with a different computational cost). It seems then natural to consider separately the following four cases: The cases 2 and 4 require to split the system as follows so t 1 will take real values. This split is similar to the one shown in the previous section by changing the roles of A and B.
In the extended phase space these two systems are equivalent to solve separately the following system written in terms of Lie operators The commutators of the Lie operators A and B, which measure the error of the splitting methods in the extended phase space is which is not proportional to ε due to the term dA(t 1 ,u) dt 1 , and this also happens with higher order commutators.
The cases 1 and 3 are associated to the split This system can be written in the extended phase space as d dt which is proportional to the small parameter ε (see [8] for more details).
Obviously, this split makes sense if one can exactly solve the non-autonomous equation associated to the dominant part du dt = A(t, u) (2.14) at a relatively low computational cost (or one can numerically solve it up to sufficiently high accuracy and at a relatively low computational cost) being the commutator-free Magnus integrators an appropriate choice in most cases. A similar methods was used in [1] for perturbed Schrödinger and Gross-Pitaevskii equations, but in those problems negative real coefficients are allowed.

Order conditions
For consistent symmetric methods we can formally write Following [6], for the composition (2.7) we have that (taking m = s) with a 0 = 0 and c s+1 = 1. The symbol ∼ indicates that, if the low order conditions are satisfied, both terms are proportional so, if the r.h.s of p aba and p abb vanish then p aba = p abb = 0. Here, the polynomial p abaaa corresponds to the dominant error term in fourth-order methods for perturbed problem. This algebraic analysis remains also valid for unbounded operators under appropriate conditions on the operators (see [16] for more details).
If we take a 1 = a m+1 = 0 we obtain a BAB composition, and the equations of p aba and p abb can be easily adjusted to obtain BAB compositions.
Second order symmetric methods which cancel the terms of order h 2p+1 ε for p = 1, 2, . . . , s and for different values of s exist with positive and real coefficients [18]. The error of these methods is of order O(h 2s+1 ε + h 3 ε 2 ) and we say the methods have effective order (2s, 2). For instance, a method which satisfies p aba = p abaaa = 0 has effective order (6,2), and this can be attained with the sequence [18] (b 1 , a 1 , b 2 , a 2 , b 2 , a 1 , b 1 ) = 1 12 , 5 − √ 5 10 , Fourth-order methods require to satisfy p aba = 0, p abb = 0, and this can not be accomplished with a i , b i real and positive valued coefficients. We are then interested on the existence of methods in which a i ∈ R + and b i ∈ C + . To get splitting methods where the coefficients satisfy these constraints we fix the values of the coefficients a i ∈ (0, 1) such that consistency and symmetry is satisfied, and leave the coefficients b i to solve the order conditions. Obviously, since the coefficients are chosen real and positive, the equations only admit complex solutions for the coefficients b i . Among all solutions obtained we will choose solutions with positive real part, i.e. b i ∈ C + from the set of all solutions found (in case these solutions exist).
Let us now analyse the number of free parameters and computational cost of ABA and BAB compositions in order to choose the most appropriate sequence: A symmetric (2k)-stage BAB composition has k coefficients a i and k + 1 coefficients b i while an ABA sequence has k + 1 coefficients a i and k coefficients b i so, the BAB composition has one more free parameter to solve the equations. In addition, since the dominant part is associated to the coefficients a i and requires the numerical solution of a non-autonomous differential equation, it is not usually possible to concatenate the last map in one step with the first one in the following step, and in practice an ABA composition with the same number of stages as a BAB composition can be computationally more costly up to one additional stage. For these reasons (number of free parameters to solve the equations and the computational cost) we only consider in this work BAB compositions.

Fourth-order methods
Fourth-order methods can be obtained with a 4-stage composition which satisfy the consistency conditions a 1 + a 2 = 1/2, 2(b 1 + b 2 ) + b 3 = 1. We can fix the values of a 1 such that a 1 ∈ (0, 1/2) and take, e.g. b 1 , b 2 to solve the equations p aba = p abb = 0. The choice a 1 = 1 4 leads to the solution obtained in [12]. However, we can take a 1 as a free parameter to minimise the dominant error term 3 min a 1 ∈(0,1/2) |Re(p abaaa )| = min 19) The corresponding system of polynomial equations with two unknowns have only two solutions (complex conjugate to each other) for each choice of a 1 and with this process we obtain following method: (2.20) A 5-stage BAB composition has the same number of coefficient b i and for this reason we have not considered it. To vanish the dominant error term at order 6 we need at least a 6-stage composition where the coefficients b i are used to satisfy the conditions p aba = p abb = p abaaa = 0. in addition to consistency.
The goal of this work is not to make an exhaustive search of methods but to show this class of methods are of interest for non-autonomous problems and to indicate how highly efficient methods could be obtained, and the optimal method can depend on the algebraic structure of each problem 4 . Then, just as an illustration we take a 1 = a 2 = a 3 = 1 6 . We have obtained one complex solution (and its complex conjugate) with coefficients:

Numerical examples
To analyse the performance of the new methods we first consider a simple non-autonomous ODE as a test bench of the methods and next we apply the methods to a linear non-autonomous PDE and a non-linear non-autonomous PDE. We compare the performance of the methods with complex coefficients versus other methods which involve real coefficients. We choose the (6,2) splitting method (2.17) which is a method of second order. As a fourth-order method we consider extrapolation (which involves substraction of quantities) where the Strang splitting symmetric second order method is used as the basic scheme to raise the order. To be more precise, we consider where we denote A i = A(t n + ih, u), B i = B(t n + ih, u). This scheme can be considered as the standard Strang decomposition applied to the non-autonomous system, but if we split it as shown in (1.12). If we take S(h) as the basic method, high order methods by extrapolation can be obtained and they only involve positive time steps. A fourth-order method is given by the composition which in our case it can be written as Φ [4] (h) = 4 3 e 4 Higher order and more efficient methods require a considerably deeper analysis, and methods belonging to this class as well as more general methods are being considered by the authors of Ref. [5].
The following schemes with real coefficients are then considered: • Strang: The second-order symmetric Strang splitting method (as a reference method); • (6,2): The symmetric splitting method of effective order (6,2) whose coefficients are given in (2.17); • (EXT4): The fourth-order extrapolation method (3.3); and the following schemes with real coefficients and a i ∈ R + are considered: • (RC4): The 4-stage fourth-order method from [12]; • (O4): The 4-stage fourth-order method built in [5], whose coefficients are available at http://www.gicas.uji.es/Research/splitting-complex.html, and referred as "Order 4 (optimized)"; • (SM4): The new optimized 4-stage fourth-order method given in (2.20); • (SM(6,4)): The new 6-stage fourth-order method whose coefficients are given in (2.22); The numerical approximations u n obtained by a given method, Ψ(h), which involve complex coefficients are computed as u n = ℜ(Ψ(h)u n−1 ), i.e. we project on the real axis after completing each time step. To measure the performance of the methods we compute the error of each method at the end of the time integration (we take as the exact solution a numerical approximation computed to a high precision) and we take as the cost of the method the number of evaluations of Φ

[h]
A which usually carries most of the computational cost.
Example 1 Let us consider the non-autonomous and non-linear perturbed equation When Ω is a constant, the system describes the motion of a charged particle in a magnetic field perturbed by s electrostatic plane waves, each with the same wavenumber and amplitude, but with different temporal frequencies ω j [11]. This equation can be written as a first order system of equations which we split as follows The linear part has, in general, no solution in closed form and we approximate its flow using the the 4th-order commutator-free Magnus integrator (2.3) which for this problem reads 5 where Ω i = Ω(t n + c i h) and the exponential of each matrix can be easily computed taking into account that (3.6) The evolution for the perturbation is immediate since both q and t are frozen.
Notice that if the computational cost is dominated by the evaluation of the time-dependent functions and the rotation matrix, then since a i ∈ R the overall cost does not change considerably either if b i is real or complex.
For the numerical experiments we take Ω(t) = 1+ 1 2 cos( 3 2 t) and the same initial conditions and parameters as given in [11]: (q 0 , p 0 , t 0 ) = (0, 11.2075, 0), s = 3 ω j = 7j. We integrate for t ∈ [0, 2π] and measure the error at the final time. All the computations are done for ε = 1/4 and ε = 1/10. Fig. 1 shows the error versus the number of evaluations for different methods. We clearly observe the superiority of the methods which consider complex coefficients versus the lower order splitting methods with real coefficients or extrapolation when high accuracy is desired as well as the high performance of the new methods. with u 0 (x) = sin(2πx) and periodic boundary conditions in the space domain [0, 1]. We take α(t) = 1 4 + µ cos(wt), V (x, t) = 1 10 (3(1 − e −t ) + sin(2πx)) and discretize in space . . , N with δx = 1/N, thus arriving at the differential equation where U = (U 1 , . . . , U N ) = (u 1 , . . . , u N ) ∈ R N . The Laplacian ∆ has been approximated by the matrix A of size N × N given by 6 and B(t) = diag(V (x 1 , t), . . . , V (x N , t)). We take µ = 1/6, w = 2, N = 100 points and compare different composition methods by computing the corre-  Fig. 2 where the superiority of the splitting methods with complex coefficients is also manifest.
The splitting considered here corresponds to solving, on one hand, the linear equation with A given by (3.9) and on the other hand, the nonlinear ordinary differential equation with initial condition u(x, 0) = u 0 (x). After discretization in space, we arrive  at the differential equation where U = (U 1 , . . . , U N ) = (u 1 , . . . , u N ) ∈ R N , A is a circulant matrix of size N × N as in the the previous linear case and F (U, t) is now defined by Here we consider splitting technique (2.11) as (3.12) where γ 1 = γ(t 1 ). Since γ(t) is frozen at real values of t, it must be considered as a constant, and the scalar equations can be solved analytically u(x, h) = u 0 (x) e γ 1 h 1 + u 0 (x)(e γ 1 h −1) , which is well defined for small complex time h. We proceed in the same way as for the previous linear case, starting with u 0 (x) = sin(2πx).
We choose β = 1, µ = 1/6, w = 2, N = 100 and compute the error at the final time t = 1 by applying the same composition methods as in the linear case. The results are collected in Fig. 3.

Conclusions
We have considered the numerical integration of non-autonomous separable parabolic equations using high order splitting methods with complex coefficients. A straightforward application of splitting methods with complex coefficients to non-autonomous problems require the evaluation of the timedependent functions on the operators at complex times, and the corresponding flows in the numerical scheme are, in general, not well conditioned. To avoid this trouble, in this work we consider a class of methods in which one set of the coefficients belong to the class of real and positive numbers. Taking the time as a new coordinate and an appropriate splitting of the system allows us to build numerical schemes where all time-dependent operators are evaluated at real values of the time. This technique shows a great interest for perturbed systems, and this problem is analysed in more detail. In this case, the flow of the dominant part has to be advanced using the real coefficients. We have analysed the algebraic structure of the problem and the cost of the algorithms in order to build efficient high order methods, and some few new methods are reported as an illustration. Higher order and more efficient methods require a considerably deeper analysis, and methods belonging to this class as well as more general methods are being considered by the authors of Ref. [5] and will be published elsewhere. Several numerical examples are considered where it is shown the good performance of this class of methods. We have shown that splitting methods with complex coefficients can also be used on non-autonomous non-linear parabolic problems and they can show a good performance so, high order and more efficient schemes following the guidelines presented in this work can be of great interest. We must also remark that order reductions are expected for problems with Dirichlet or Newman boundary conditions on bounded domains, and the performance of high order methods on these problems diminishes, being an interesting problem that needs further investigation.