Explicit volume-preserving splitting methods for divergence-free ODEs by tensor-product basis decompositions

We discuss the construction of volume-preserving splitting methods based on a tensor product of single-variable basis functions. The vector field is decomposed as the sum of elementary divergence-free vector fields (EDFVFs), each of them corresponding to a basis function. The theory is a generalization of the monomial basis approach introduced in Xue & Zanna (2013, BIT Numer. Math., 53, 265–281) and has the trigonometric splitting of Quispel & McLaren (2003, J. Comp. Phys., 186, 308–316) and the splitting in shears of McLachlan & Quispel (2004, BIT, 44, 515–538) as special cases. We introduce the concept of diagonalizable EDFVFs and identify the solvable ones as those corresponding to the monomial basis and the exponential basis. In addition to giving a unifying view of some types of volume-preserving splitting methods already known in the literature, the present approach allows us to give a closed-form solution also to other types of vector fields that could not be treated before, namely those corresponding to the mixed tensor product of monomial and exponential (including trigonometric) basis functions.


Introduction
Volume preservation is an important property shared by several dynamical systems.The principal example is the velocity field for an incompressible fluid flow and applications include the study of turbulence and of mixing of flows (see Arnold & Khesin, 1998).The preservation of phase-space volume is also a key ingredient in methods that involve stochastics like the hybrid Monte Carlo algorithm (see Neal, 2011).Thus, preservation of volume by a numerical method for differential equations is a desirable property either in conjunction with other properties, like energy preservation or time reversal symmetry, or on its own.The latter case is the subject of this paper.
Designing volume-preserving numerical integrators is a hard task, especially given the existence of no-go theorems within the class of 'standard' methods: in Qin & Zhu (1993) it was proved that no standard numerical method, for example a Runge-Kutta scheme, is volume preserving for all linear vector fields.A similar result was proved in Kang & Shang (1995).Further, Iserles et al. (2007) and Chartier & Murua (2007) proved that no B-series method can be volume preserving for all possible divergence-free vector fields.The only methods in this class that can be volume preserving are the composition of exact flows of divergence-free subsystems.Thus, the design of efficient methods that preserve volume is still an open problem in geometric integration (see McLachlan & Scovel, 1998).
Despite the existence of no-go theorems, volume-preserving numerical methods can be constructed by other techniques, for instance using generating functions (see e.g.Scovel, 1991;Shang, 1994;Quispel, 1995;Xue & Zanna, 2014) or by splitting methods.One of the earliest volume-preserving splitting methods is the splitting method by Kang & Shang (1995), decomposing the vector field into the sum 2 of 18 A. ZANNA of essentially two-dimensional Hamiltonian fields, which are then solved by a (typically implicit) symplectic method.
Because of the difficulty of addressing the general space of divergence-free vector fields, recent efforts have concentrated on smaller, yet still interesting, function spaces, for instance, the space of polynomial fields.An earlier paper on the splitting of polynomial vector fields is McLachlan & Quispel (2004), dealing mainly with the Hamiltonian case.In Xue & Zanna (2013), we presented a novel approach for arbitrary polynomial divergence-free vector fields, by expanding the divergence equation in terms of the monomial basis, that allowed us to develop explicit volume-preserving splitting methods.It turns out that the approach can be generalized to an arbitrary tensor basis framework, which embraces several known splitting methods as special cases, namely the monomial splitting of Xue & Zanna (2013), the splitting in shears of McLachlan & Quispel (2004) and the trigonometric splitting of Quispel & McLaren (2003).The main idea is simple: to expand the divergence equation in terms of a tensor product of one-dimensional basis functions.Each basis function gives rise to a uniquely determined (up to an off-diagonal contribution) elementary divergence-free vector field (EDFVF).Away from some critical points, these EDFVFs are proved to be integrable.Further, we identify solvable EDFVFs, recovering the monomial basis and the exponential basis (and the mixed tensor product of these).The resulting composition methods are thus explicit and volume preserving.
Finally, we also consider expansions in EDFVFs for divergence-free fields with measure μ(x) = 1, and show that the noncanonical case can be viewed as a time-reparametrization of the canonical one.Thus, the results and methods presented in this paper can be applied to the case ∇ • (μf) = 0.

The basis-decomposition approach
We consider the problem of solving the differential equation (2.1) subject to the divergence-free condition by splitting the vector field f as the sum of EDFVFs that are integrable by analytical methods.As the exact flow of each of the elementary vector fields is itself volume preserving, the resulting composition method is volume preserving.
In McLachlan et al. (2008), it was proposed to decompose any divergence-free vector field f(x) into a diagonal and off-diagonal part, where, componentwise, From the definition of divergence, it is clear that only the coefficients of the diagonal part are directly involved in the divergence-free condition (2.2); therefore, vector fields with zero diagonal part are automatically divergence-free.Further, there exists a generic volume-preserving splitting technique to solve for these off-diagonal vector fields, namely, the splitting into canonical shears: one integrates one variable at a time, keeping the remaining variables constant; see McLachlan & Quispel (2004) and McLachlan et al. (2008).Vector fields with nonzero diagonal are more difficult to treat; thus,

of 18
without loss of generality, we assume that the function f has zero off-diagonal part, i.e.
Proposition 2.1 Let f(x) be a given vector field, with divergence p(x) = ∇ • f(x), and assume that the latter can be expanded in a set of basis functions {φ j (x)} j∈J , For each j ∈ J , there exists a uniquely determined vector field F j (x) such that is also divergence-free and, moreover, the decomposition f(x) = j F j (x) is unique.
Proof.Since the basis functions are independent, for all indices j ∈ J whenever f obeys (2.2).Consider next vector fields F j , constructed so that where (a j ) i is the component of ∂ i f i (x) along φ j (x).Note that the F j s are uniquely defined up to a function h i (x) independent of x i , which we can choose to be zero, as it would only contribute to the off-diagonal part.By construction, F j (x) is such that ∇ • F j (x) = p j φ j (x) = 0 since p j = 0. Hence, F j is divergence-free.Finally, by collecting all terms in f by picking out those which contribute to each p j , we are able to split f as the sum of EDFVFs, f = j∈J F j .
A practical way to find the components of the vectors a j in (2.6) is the following: let J [0] = {} be the initial (empty) list of indices for the basis functions.Starting with i = 1, for i = 1, . . ., n, expand ∂ i f i in the given basis.If new basis functions are added to the expansion, the corresponding indices are added to the list of indices J [i] .For each index j ∈ J [i] , the coefficient (a j ) i is determined.When i = n, J [n] = J is the final list of indices of the expansion and the a j are completely determined.
It might happen that some of the J [i] s, hence J , are infinite.This is the case when some of the ∂ i f i s expand in infinite terms.In practice, it is desirable to work with a finite number of basis elements: we assume that each ∂ i f i can be projected to the same finite-dimensional space spanned by a finite number of basis functions.Correspondingly, we obtain a projection of the vector field f onto a finite-dimensional set of divergence-free vector fields.The dimension of such a set is ideally chosen by the problem owner not to alter the dynamic of the system of interest in a significant manner.Alternatively, we assume 4 of 18 A. ZANNA that the remainder of the finite-dimensional projection can be treated by some other volume-preserving numerical method.Thus, without further ado, we will assume J to be finite.Definition 2.2 Given the basis {φ j (x)} j∈J , the vector field F j (x) obeying (2.4) will be called an EDFVF associated to the basis function φ j (x).
The main idea of the paper can be summarized as follows: decompose the divergence function using the basis functions {φ j (x)} j∈J , p(x) = j∈J p j φ j (x).
For each j ∈ J , • identify the terms in f(x) that contribute to p j to obtain F j (x) (EDFVF); • integrate ẋ = F j (x) exactly (if possible) or by some volume preserving method.
The composition resulting from the flows of the F j is then volume preserving by construction.The main goal is to find suitable basis functions, in the sense that the corresponding vector fields must be easy to integrate exactly.

Tensor-product bases
Consider the case when the basis functions are the tensor product of one-dimensional function bases, (3.1) Observe that the choice corresponds to the case of the monomial basis for polynomials, discussed in detail in Xue & Zanna (2013).
Proposition 3.1 Let φ j (x) be in the form (3.1), for some arbitrary tensor-product basis.The EDFVF F j (2.5) corresponding to the multi-index j = (j 1 , . . ., j n ) can be written as where the a j ∈ R n obey the condition

of 18
Proof.The result is a corollary of Proposition 2.1, specific to the case of tensor bases.Note, in particular, that upon taking the divergence, the coefficients (a j ) i obey the divergence-free condition a j T 1 = p j = 0; thus (3.3) holds.
It is worth mentioning that the second line in (3.2) makes sense when φ j i (x i ) = 0.This motivates the following definition.
Definition 3.2 The points in which φ j i (x i ) = 0 are called critical points of the basis function φ j i .Similarly, the points in which the φ j i (x i ) dx i = 0 are the critical points of the primitives of the basis functions.The collection of both is referred to as the critical points of the EDFVF F j , or simply critical points.
From now on, for ease of notation, we will discard the dependence of the vector a j on the index j, unless the index is needed for the sake of clarity.
i is an integral of the system, i.e.
for all t.In particular, away from the critical points, the system has n − 1 independent integrals and is integrable.
Proof.Let b l be as above.Note that as a consequence of the orthogonality of the vectors.Hence, it follows that log n i=1 As far as integrability is concerned, it is sufficient to show that the integrals are independent, i.e. that the linear combination For this purpose, observe that for any l = 1, . . ., n − 1.Hence, the linear system (3.5) can be written in the form where Ξ and J are the diagonal matrices with elements and dimensions n × n and (n − 1) × (n − 1), respectively.It is immediate to see that the matrix [b 1 , . . ., b n−1 ] has full rank.The matrix J has full rank as long as each of the integrals I l is nonzero, i.e. each of the quantities φ j i (x i ) dx i = 0. Similarly, the matrix Ξ has full rank, provided that φ j i (x i ) = 0. Away from the critical points of the basis functions and their primitives, the system has full rank and the integrals are independent.Thus, the system is integrable.
Definition 3.4 (Diagonalizable EDFVFs) Consider the EDFVF F j corresponding to the multi-index j = (j 1 , . . ., j n ), away from its critical points.We say that F j (x) is diagonalizable if φ j (x) has a closed-form expression.
Diagonalizability means that φ j (x) can be treated as a function of t only.Thus, each x i obeys a differential equation with variable coefficients, indicating that one might find classes of diagonalizable EDFVFs whose solution can be given in closed form.For this purpose, we study the behaviour of the basis function φ j (x) as a function of time.
Theorem 3.5 (Solvable diagonalizable EDFVFs) Consider the EDFVF F j away from its critical points.The tensor basis function We commence by analysing the constant case, where c i = 0 (the case c i = 0 corresponds to constant basis functions and is uninteresting).In this case, the basis functions obey the differential equation which has solution where the dash denotes differentiation with respect to the variable z.Dividing both sides by ΦΦ and integrating, we obtain log where k i is an integration constant.By exponentiation, we obtain and here we distinguish two cases: , where l i is another constant of integration.This corresponds to the basis function φ j i (x i ) = l i k i e k i x i and it is useful to choose l i k i = 1.The corresponding basis function is then (3.12) For this special choice of basis functions, c i = 1 for all i and the divergence-free condition (3.3) imply that the right-hand side of (3.8) is equal to zero; hence A. ZANNA Next, we consider the case when c i = 1.Integrating (3.11) between 0 and z, we obtain where, again, k i is a constant of integration.The constants l i depend on the choice of the value of Φ at z = 0. Their choice will be discussed later.Thus, we have To recover the basis functions, we differentiate Φ(z) to obtain The integration constants k i can be chosen such that k i (1 With this choice, we have In particular, if l i = 0 and m i = j i ∈ Z, we recover the monomial basis studied at great length in Xue & Zanna (2013).The case l i = 0 is just a shifted monomial basis.In particular, the above formula indicates that it is also possible to use noninteger values of m i .Note that the value m i = −1 is not allowed.This value of the exponent (corresponding to 1/x = log x) has also been discussed in Xue & Zanna (2013).
The general case g(t, y), ∂ y g = 0, is not solvable for all a satisfying (3.3).The argument is the following.Assume that By differentiation with respect to x i , we obtain It is readily seen that the left-hand side depends solely on x i , while the right-hand side also depends on all the other variables x l , l = i.Hence, the above equation cannot have a nontrivial solution, unless for very specific choices of a.
Remark 3.6 In fact, the previous theorem classifies all cases of diagonalizable tensor basis functions having a closed-form solution: the monomials and the exponential basis.Other choices might still be integrable, but they might not necessarily give a closed-form solution for all a satisfying (3.3).As an example, the cosine basis φ j i (x i ) = cos j i x i and the sine basis φ j i (x i ) = sin j i x i do not give rise to a solvable tensor basis for all a, as (3.7) is not satisfied.Yet, we know that the corresponding EDFVF is integrable by Theorem 3.3 away from critical points; see also Fig. 1.Combinations of sines and cosines can be treated with solvable EDFVFs by transforming into an exponential basis (Fourier basis), as described in the next section.
Remark 3.7 Note that the two classes of tensor basis functions satisfying Theorem 3.5, monomials and exponentials, are precisely the classes of functions that are closed under differentiation, i.e. any partial derivative of a monomial/exponential is still a monomial/exponential. Further, they are also closed under
multiplication: the product of monomials is still a monomial, the product of exponentials can still be written as a single exponential.This has the important consequence that the commutator of these vector fields will still be a vector field of the same type, for which the same algebraic formula holds.This is of relevance for the construction of higher-order integrators.Typically, symmetric compositions, using only the splitting of the vector field, are used.Here, we might increase the order by also using the commutators.

Connection with the splitting in shears of Quispel and McLachlan
Quispel and McLachlan introduced the technique of splitting (2.1) in divergence-free vector fields of the type ẋ = af (j T x), j T a = 0 (3.13) (splitting into shears), with f an arbitrary scalar function; see McLachlan & Quispel (2004) and McLachlan et al. (2008).Here, j need not be a vector of integers.It is easy to verify that j T x = const along the flow of (3.13); hence the solution is a simple translation parallel to the vector a.Since ẋ is parallel to a, this also implies that b T x = const for any b orthogonal to a.
Proposition 3.8 The splitting into shears (3.13) is equivalent to a diagonal tensor basis splitting, with exponential basis choice (3.12).
Proof.Set f (j T x) = Ke j T x (which is always possible, since j T x is constant along the flow).Then, we can rewrite (3.13) as ẋi = Ka i j i 1 j i e j T x , i = 1, . . ., n.
(3.14) 10 of 18 A. ZANNA Above, we have tacitly assumed that j i = 0 for i = 1, . . ., n.This is no real restriction as, if j i = 0 for some i, the remaining variables do not depend on j i and the system is reducible.Note that the condition a T j = 0 is now equivalent to ãT 1 = 0, where ã = Ka ⊗ j (tensor product).Thus, we see that (3.14) is precisely of the form (3.2) with basis functions φ j i (x i ) = e j i x i and coefficient vector ã obeying (3.3).Thus, Theorem 3.3 holds, and for any vector b orthogonal to ã, the quantity In particular, by choosing b = 1 we recover the condition j T x = const.For the other choices of b, taking the logarithm, we see that i bi (j i x i − j i ); hence i bi j i x i , is a constant of integration.Next, set b i = bi j i .We have thus proved that d dt b T x = 0.Moreover, the condition bT ã = 0 is equivalent to b T a = 0, and the condition bT 1 = 0 is equivalent to b T j = 0, and thus there is complete equivalence of the two formulations.

Trigonometric polynomials and Fourier bases
We next turn our attention to the case of trigonometric polynomials, which often arise when the vector field is expanded in a Fourier series.A typical application is the study of the dynamics of incompressible inviscid fluid flows.Some practical examples, often written in their real form using trigonometric polynomials, are the Arnold-Beltrami-Childress flow (ABC flow, Ricca, 2001), a prototype for the study of turbulence, and ẋ1 = − cos(π x 1 ) sin(π x 2 ) + ε sin(2π x 1 ) sin(π x 3 ), ẋ2 = sin(π x 1 ) cos(π x 2 ) + ε sin(2π x 2 ) sin(π x 3 ), a model for the study of the mixing in a laminar vortex flow with rigid boundaries, occurring in many geophysical, industrial and biophysical applications; see Solomon & Mezić (2003).For our analysis, it is more convenient to work in the complex setting.Trigonometric polynomials can always be written in a complex exponential form using the Euler formulas (4.5).

The method of expansion in a Fourier series
Consider again a vector field ẋ = f(x) and assume that each of the function components f i (x) is expanded in a Fourier series.The result will be a right-hand side that is the product of one-dimensional Fourier series in each variable, i.e.
where i 2 = −1 defines the imaginary unit.The divergence is given as

of 18
We note at once that if the vector field f is expressed in a Fourier basis, then the divergence is also expressed in a Fourier basis.To put this in the formalism of the earlier sections, we identify , so that the relation between the vectors a j in Proposition 3.1 and c j is the normalization Corollary 4.1 Under the assumptions above, f divergence-free implies c j T j = 0 for any multi-index j, or, equivalently, a j T 1 = 0, in the nD Fourier expansion of f.
We state the above result as a corollary because it was proved by Quispel & McLaren (2003) by direct computation on trigonometric polynomials.It is a direct consequence of Theorem 3.3.
Similarly to the approach presented for polynomials, we split the vector field into EDFVFs, each F j associated to a tensor basis element φ j .Because of Theorem 3.5, the basis function is solvable.In fact, we have d dt j T x = j T ẋ = (j T c j )e ij T x = 0, since c j and j are orthogonal (cf. the divergence-free condition).This implies that (4.1) decomposes in EDFVFs of the form ẋ = F j (x) = c j e ij T x , whose solution is given explicitly as A naive implementation of this method would require complex arithmetic.However, if the function f is real, then the coefficients of the Fourier transform obey the condition c −j = cj (the bar denotes complex conjugation).This implies that the EDFVFs F j and F −j share the same integrals and can be combined together.Thus, we need to sum over just half of the multi-indices and, instead of (4.1), we can consider ẋ = Fj (x) = c j e ij T x + cj e −ij T x .(4.4) Writing c j = α j + iβ j , and using the identities it is readily seen that (4.4) becomes which can be integrated exactly as 12 of 18 A. ZANNA

Tensor product of Fourier series and the monomial basis
So far, the pure monomial choice reduces to the methods in Xue & Zanna (2013) and the pure trigonometric choice reduces to the methods of Quispel & McLaren (2003).However, because of the generality of Proposition 2.1 we are also able to treat the mixed polynomial-trigonometric case, which one could not deal with systematically before.We start from the divergence function, which we assume can be written in the form where, by an appropriate relabelling of the variables, we can assume Here, the subscript m is associated to the projection onto the first p components (monomial part), while the subscript f refers to the projection onto the remaining n − p components (Fourier part).Dropping the dependence of the coefficient a on the index j, the corresponding EDFVF is x j m e ij T f x , a T 1 = 0 (5.1) (the ⊗ denoting the tensor (componentwise) product of two vectors).Equivalently, we could write the EDFVF using the unscaled coefficients c, subject to c T ((j + 1) m + ij f ) = 0.
(5.3) From Theorem 3.5, it follows that the basis function φ j (x) = x j m e ij T f x obeys the differential equation Now, for i = 1, . . ., p, we have ẋi = c i x i φ j (t), with solution (5.4) while, for i = p + 1, . . ., n we have ẋi = c i φ j (t), and r i as in (5.4).

Avoiding complex arithmetic
In the case when basis elements x j m e ij T f x are used, the coefficients are generally complex and the method, as described before, will require complex arithmetic.This might be a problem, because the logarithm functions in (5.5) will require the tracking of the argument function, for a correct complex part.As most applications involve real coefficients, it is a disadvantage (and computationally more expensive) to use complex arithmetic.
A first attempt would be to try to generalize the techniques for pure trigonometric tensor products described earlier (and in Quispel & McLaren, 2003).Let c = α + iβ.It is readily seen that the divergence-free condition (5.3) is equivalent to the two conditions α T (j + 1) m − β T j f = 0, (5.6) (5.7) These are, in fact, the conditions for the vector fields (5.8) to be volume preserving.The vector field T 1 + T 2 can be recognized as the real part of (5.2) (we assume that the solution is real).Note the resemblance to (4.6), which we treated simultaneously (because the sine and cosine part shared the same integral).Here, the vector field T 1 has integral x (j+1) m sin(j T f x), while the vector field T 2 has integral x (j+1) m cos(j T f x) .These two integrals do not commute, indicating that there is no easy way to solve T 1 + T 2 simultaneously.On the other hand, we observe that A. ZANNA x j m cos(j T f x) and x j m sin(j T f x) are also basis functions (though not tensor product, but linear combinations of tensor bases).Thus, a natural strategy is to attempt solving T 1 and T 2 separately.
We set It is readily observed that (5.9) with (5.10) The two equations in (5.9) can be combined into which can be integrated by separation of variables to give and substituting φ m (t) into φ f , we obtain the differential equation which can be solved by separation of variables.The solution is given in terms of hypergeometric functions, (5.12) see Abramowitz & Stegun (1965) and Rainville (1967).Some special cases of the parameter ρ are worth discussion.
Case ρ = 1 (κ 1 = 2κ 2 = 0).In this case, (5.11) reduces to φ f = K.This gives . Thus, we have for the monomial part.For the Fourier part, we focus on T 1 and observe that y = sin(j T f x) obeys the same differential equation as x i of the monomial part upon replacement of α i by −κ 2 .Thus, it has a similar solution: A similar formula holds for T 2 .Substituting into ẋi and integrating, we have where T 1 : q i = −β i s 0 |s 0 |, T 2 : q i = α i c 0 |c 0 |.
The general case requires the inversion of the special function in (5.12).Note that if ρ > 1, the series terminates (finite sum).

Divergence-free with measure μ(x) = 1
In the above, we have discussed the case of how to split divergence-free methods for the case of canonical variables, i.e. if the vector field f obeys In some cases, the divergence-free condition is given with respect to a measure μ(x), meaning that the vector field ẋ = f(x) obeys the condition

of 18
The scope of this section is to extend the results for the canonical case μ = 1 to the generic case when μ = μ(x).A result very similar to the case μ = 1 is true also in this case.
Proposition 6.1 Let φ j (x) be in the form (3.1) for some arbitrary tensor-product basis.Consider a decomposition ∇ • (μf) = j p j φ j (x).
The EDFVF F j (2.5) corresponding to the multi-index j = (j 1 , . . ., j n ) must be of the form ẋi = (a j ) i 1 μ(x) φ j i (x i ) dx i n l=1 l =i φ j l (x l ) (6.2) . . . , n, (6.3)where a j T 1 = 0, as in (3.3).Furthermore, for any b orthogonal to a j and 1, one has that is a conserved quantity of the system.
The proof is analogous to the case μ = 1.Thus, also, the generic case is integrable, provided that the integrals are linearly independent.Proposition 6.2 The vector field (6.3) is equivalent to the vector field (3.2) upon time reparametrization dt dτ = μ(x).Thus, all the knowledge of the canonical case μ = 1 can in principle be extended to the more general case μ = 1.
Theorem 3.3 (Integrability of EDFVFs of tensor bases) Consider the EDFVF (3.2) associated to the tensor basis function φ j (x) (a j ≡ a).Let b 1 = 1 and b 2 , . . ., b n−1 be such that the set of vectors {a, b 1 , . . ., b n−1 } are mutually orthogonal.Then, for any l = 1, . . ., n − 1, the quantity from which the first part of the statement follows by exponentiation of both sides of the equation.In passing, note that (3.4) holds also for any b in the span of b 1 = 1, b 2 , . . ., b n−1 , i.e. for any b orthogonal to a.