You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: |
    [ In this paper, we present new, explicit, volume-preserving vector fields for polynomial divergence-free vector fields of arbitrary degree (both positive and negative). The main idea is to decompose the divergence polynomial by means of an appropriate basis for polynomials: the monomial basis. For each monomial basis function, the split fields are then identified by collecting the appropriate terms in the vector field so that each split vector field is volume preserving. We show that each split field can be integrated exactly by analytical methods. Thus, the composition yields a volume preserving numerical method. Our numerical tests indicate that the methods compare favorably to standard integrators both in the quality of the numerical solution and the computational effort. ]{}

    [**Keywords:** Geometric integration; volume preservation; splitting methods]{}
author:
- '[Huiyan Xue and Antonella Zanna[^1]]{}'
title: 'Explicit Volume-Preserving Splitting Methods for Polynomial Divergence-Free Vector Fields'
---

Introduction
============

Divergence-free vector fields occur naturally in incompressible fluid dynamics, and preservation of phase-space volume is also a crucial ingredient in many, if not all, ergodic theorems. Preservation of volume by a numerical method for differential equations is thus a desirable property in the study of dynamical system. Nevertheless, designing volume preserving numerical integrators is a hard task, as the space of divergence-free vector fields seems to be too large. [@qin93vps] proved that no standard numerical method, for example a Runge–Kutta scheme, is volume preserving for all such vector fields. Thus, the design of efficient methods that preserve volume is still a standing open problem in geometric integration [@Mclachlan98asurvey].

Despite no-go theorems for volume preservation within the class of “standard” methods (see also the recent results in [@MR2334044], [@chartier2007preserving], where it is proved that no B-series method can be volume preserving for all possible divergence-free vector fields), it is known that volume preservation can be achieved, either restricting the class of vector fields or using methods other than B-series. [@chartier2007preserving] have shown that volume-preserving B-series can be obtained if the vector field has some specific dependence (for example $f_1=f_1(x_2), f_2 = f_2(x_3), \ldots, f_n = f_n(x_1)$, or if the variables naturally decompose into two sets, $\mathbf{x} = [\mathbf{y}, \mathbf{z}]^T $, obeying $\dot{\mathbf{y}} = \mathbf{g}(\mathbf{z}), \dot{\mathbf{z}} = \mathbf{h}(\mathbf{y})$). These two examples of vector fields correspond to what we call *off-diagonal*. These are generally easier to treat and other methods will be described in the sequel.

Volume-preserving maps can be constructed techniques other than methods that possess B-series, for instance using generating functions [@MR1370428; @MR1290253; @scovel91sni]. This technique involves evaluating definite integrals of the vector field, and, in addition, the method does not preserve fixed points. Another technique is based on splitting methods. One of the earliest volume-preserving splitting methods is indeed the splitting method by [@kang95vpa], decomposing the vector field into the sum of essentially 2-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 to smaller, yet still interesting, functions spaces, for instance the space of polynomial fields. An earlier paper on splitting polynomial vector fields is by [@mclachlan04egi]. That paper had some discussion of the divergence-free case, but mainly dealt with the Hamiltonian case. Investigations of the Hamiltonian case, which involves expressing a scalar polynomial of degree $d$ in $n$ variables as a sum of functions of fewer variables, have shown that good splitting methods exist, but that finding and analyzing them (especially for general $n$ and $d$) is very difficult [@dragt96sma; @mclachlan04egi; @rangarajan96sco]. The volume-preserving case, which involves $n$ polynomials subject to the divergence-free condition, is even harder, although there is a conjecture by [@mclachlan04egi] that they can be expressed as a sum of $n+d$ shears, each a function of $n-1$ variables. The case of linear and quadratic divergence-free vector fields was studied in detail in [@mclachlan09evp], where several explicit volume-preserving splitting methods were introduced. In that paper, two main classes of methods were considered: a) methods that distinguish the diagonal and off-diagonal part and b) methods that do not. By diagonal part we mean all the terms of the vector field such that $\dot x_i$ depends on $x_i$, for $i=1, \ldots, n$. Similarly, the off-diagonal part refers to all the terms of the vector field such that $\dot x_i$ does not depend on $x_i$, $i=1,\ldots, n$. Furthermore, for the class a), several explicit schemes that treated the diagonal part and the off-diagonal part separately were introduced and tested numerically. Numerical tests indicated that methods that treated the diagonal part by splitting it in terms treated by exponentials had a smaller error than methods splitting in shears.

In this paper we present a new approach that allows us to develop explicit volume-preserving methods for arbitrary polynomial divergence-free vector fields, including those with negative degree. The first main insight is to expand the divergence equation, rather than the vector field, in the monomial basis. For each monomial basis element, we identify the elements in the vector field associated to it to construct a divergence-free elementary vector field. The second insight is to recognize that the elementary divergence-free vector fields can all be treated by the same formalism and therefore can be solved explicitly by elementary analytical methods. The split fields are then composed to obtain explicit first order method and second order method (by symmetrization). The resulting composition method is thus explicit and volume-preserving. Being explicit, the proposed method are computationally efficient and possess excellent qualitative properties. We believe that thisis a consequence of volume preservation solely, as the methods are not necessarily time-reversible nor self-adjoint (for instance, the first order method is neither).

The paper is organized as follows. In Section \[sec:background\] we review some background and introduce notation. In Section \[sec:hom\] we present the monomial basis for polynomial volume-preserving vector fields and prove that the vector fields associated to the basis elements can be integrated exactly. The case of polynomials of negative degree is treated in Section \[sec:laurent\]. In Section \[sec:num\] we test the application to volume-preserving cubic Stokes flows. Finally, Section \[sec:conclusion\] is devoted to some concluding remarks and some directions for future research.

Background and notation {#sec:background}
=======================

We consider the ordinary differential equation $$\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}), \qquad \mathbf{x}(0) = \mathbf{x}_0,
    \label{eq:vf}$$ where $\mathbf{x} \in \RR^n$ and $\mathbf{f} :\RR^n \to \RR^n$, $\mathbf{f}(\mathbf{x}) = [f_1(\mathbf{x}), \ldots, f_n(\mathbf{x})]^T$, is subject to the divergence-free condition $$\nabla \cdot \mathbf{f} = \sum_{i=1}^n \partial_{x_i} f_i (\mathbf{x})=0.
    \label{eq:df}$$ An arbitrary vector field $\mathbf{f(x)}$ can always be decomposed as $$\mathbf{f(x)}= \mathbf{f}^{\mathrm{diag}} (\mathbf{x}) +\mathbf{ f}^{\mathrm{offdiag}} (\mathbf{x})$$ where, component-wise, $f_i^{\mathrm{diag}} (\mathbf{x})$ is the collection of terms in $f_i(\mathbf{x})$ that depend on $x_i$ ( i.e. $\partial_{x_i}  f_i^{\mathrm{diag}} (\mathbf{x})\not=0$). Similarly, $\mathbf{ f}^{\mathrm{offdiag}} (\mathbf{x}) $, is given, component-wise, by the collection of terms in $f_i$ that do not depend on $x_i$, i.e. $ \partial_{x_i} f_i^{\mathrm{offdiag}} (\mathbf{x})=0$. We refer to $\mathbf{f}^\mathrm{diag}$ and $\mathbf{ f}^{\mathrm{offdiag}} (\mathbf{x}) $ as the *diagonal part* and the *off-diagonal* part of the vector field $\mathbf{f}$, respectively.

From the definition of divergence, it is clear that only the coefficients of the diagonal part are directly involved in the divergence-free condition , therefore, vector fields with zero diagonal part are automatically divergence-free.

The off-diagonal part of a vector field is generally easier to treat by volume-preserving methods. The method of splitting in canonical $n$-shears is generic (for each $i =1\ldots n$, we solve for $\dot{x}_i$, while $\dot{x}_k=0$, $k\not=i$). For polynomial systems, it is possible to construct splittings in lower-triangular systems (for each $i =1\ldots n$, we solve for $\dot{x}_i$ depending only on $x_1, \ldots, x_{i-1}$, under a suitable permutation of the indices), see [@mclachlan09evp].

Thus, to obtain volume-preserving methods, special care has to be taken when splitting the diagonal part, because splitting across the conditions for zero-divergence will yield methods that are not volume-preserving. The goal is to split the diagonal part in divergence-free vector fields that can be solved exactly. These observations motivate the following assumption.

Unless otherwise stated, we will assume that the given vector field is *diagonal* only (i.e. $\mathbf{f}=\mathbf{f}^\mathrm{diag}$) and that the off-diagonal part is zero ($\mathbf{f}^\mathrm{offdiag}=\mathbf{0}$).

In this paper we consider the case when the functions $f_i(\mathbf{x})$ are polynomials. To treat the general case with $n$ variables and given degree $d$, it is convenient to introduce a multi-index notation.

Let $\mathbf{j} = (j_1, \ldots, j_n)\in \mathbb{N}^n$ be a multi-index and let $\mathbf{x}^\mathbf{j} = x_1^{j_1} x_2^{j_2}\cdots x_n^{j_n}$ (monomial). In addition, denote by $|\mathbf{j}| = j_1+\cdots + j_n$ the degree of the monomial.

Earlier studies of splitting methods for polynomial fields have focussed on homogeneous polynomials. Assume that the $f_i(\mathbf{x})$ are homogeneous polynomials of degree $d$ in the variables $x_1, x_2, \ldots, x_n$, i.e.$$\dot x_i = f_i(\mathbf{x})=\sum_{|\mathbf{j}|=d}  a^i_{\mathbf{j}}\,  \mathbf{x^j}, \qquad i = 1,\ldots, n.
        \label{eq:polyd}$$ The case when $f_i(\mathbf{x})$ is not homogeneous is easily treated by a further splitting in homogeneous terms. Define $N(n,d) = {n+d-1\choose d}$. In [@mclachlan04egi] it was shown that $N(n,d)$ is the number of coefficients $a^i_\mathbf{j}=a^i_{j_1, \ldots, j_n}$ (for a fixed value of $i$) of an homogeneous polynomial $f_i$ of degree $d$ in $n$ variables. The divergence-free condition is a homogeneous polynomial of degree $d-1$ because of derivation. All the coefficients of this polynomial must be identically equal to zero, because of the arbitrariness of the variables. This gives $N(n, d-1)$ conditions for volume preservation. In the polynomial case, there are $N(n, d-1)$ divergence-free conditions on the coefficients. In order to be volume preserving, a splitting must satisfy one or more conditions on the coefficients. Our approach is inspired in part from that argument. The new idea is to consider a polynomial vector field and to look at a monomial-basis expansion of the divergence function. Each basis element will correspond exactly is associated to a condition on the coefficients of the vector field. As long as we split according to basis terms, we are guaranteed to obtain “elementary” monomial fields that are also divergence free. We show that all these “elementary divergence-free vector fields” possess $n-1$ integrals in evolution and are therefore integrable. Furthermore, we give the explicit analytical solution. This allows us to generate splitting methods for arbitrary divergence-free polynomial fields.

We illustrate the general idea by discussing a simple 2-dimensional example.

Consider the vector field $$\begin{aligned}
    \dot x_1 &=& a^1_ \mathbf{j} \, x_1^3 + a^1_\mathbf{k}  \, x_1^2 x_2 + a^1_\mathbf{l}  \, x_1 x_2^2\\
    \dot x_2 &=& a^2_\mathbf{j}  \, x_1^2 x_2 + a^2_\mathbf{k}  \, x_1 x_2^2 + a^2_\mathbf{l}  \, x_2^3\end{aligned}$$ The divergence of this vector field is $$p(x_1,x_2) = (3 a^1_\mathbf{j} + a^2_\mathbf{j}) x_1^2 + (2 a^1_\mathbf{k} + 2 a^2_\mathbf{k}) x_1 x_2 + ( a^1_\mathbf{l}+ 3a^2_\mathbf{l})x_2^2.$$ The divergence free condition becomes the set of equations $$3 a^1_\mathbf{j} + a^2_\mathbf{j}=0, \qquad a^1_\mathbf{k} + a^2_\mathbf{k}, \qquad a^1_\mathbf{l}+ 3a^2_\mathbf{l} =0.$$ As argued above, for each equation we obtain a divergence-free split vector fields $$\begin{array}{lclcl}
    \dot x_1 =a^1_ \mathbf{j} \, x_1^3  &   & \dot x_1 =  a^1_\mathbf{k}  \, x_1^2 x_2  &  & \dot x_1 = a^1_\mathbf{l}  \, x_1 x_2^2\\
    \dot x_2 =   a^2_\mathbf{j}  \, x_1^2 x_2&  &\dot x_2 = a^2_\mathbf{k}  \, x_1 x_2^2  ,      &  &\dot x_2 =  a^2_\mathbf{l}  \, x_2^3 .     
    \end{array}$$ Surely, each of these vector fields can be integrated. For instance, in the first vector field, one can solve for $x_1$ then substitute into the second equation to solve for $x_2$. In the third vector field, the procedure is similar, just interchange the role of $x_1$ and $x_2$. The second vector field can also be integrated, though not exactly by the same method. The situation becomes far more complicated for several variables and higher order polynomial vector fields, so that, at first glance, this procedure does not yield a method that can be easily generalized. However, observe that the three vector fields can be written as $$\begin{array}{lclcl}
    \dot x_1 = a^1_ \mathbf{j} \,  x_1\cdot x_1^2   &   & \dot x_1 = a^1_\mathbf{k}  \, x_1 (x_1 x_2) &  & \dot x_1 = a^1_\mathbf{l}  \,  x_1 \cdot x_2^2 \\
    \dot x_2 = a^2_\mathbf{j}  \, x_2 \cdot x_1^2   &   &\dot x_2 =  a^2_\mathbf{k}  \,  x_2 (x_1 x_2)   &  &\dot x_2 = a^2_\mathbf{l}  \,x_2 x_2^2,    
    \end{array}$$ i.e.  each vector field obeys an equation of the type $\dot x_i = c_i x_i (x_1^{j_1} x_2^{j_2})$ where $x_1^{j_1} x_2^{j_2}$ are the first, the second and third monomial in $p(x_1,x_2)$. In particular, if $ (x_1^{j_1} x_2^{j_2})(t)$ is known, $x_i$ can be obtained by integration.

In what follows, we shall see that this argument is generic: it yields for any number of variables and any degree of the vector field. We shall therefore identify monomials appearing in the divergence-free condition and solve for them explicitly. When the monomials, as function of time, are known, the remaining variables follow by integration of linear equations with variable coefficients.

Elementary divergence-free vector fields: the polynomial case {#sec:hom}
=============================================================

Consider a polynomial vector field $$\dot {x}_i = f_i(\mathbf{x})=\sum_{|\mathbf{k}|=1}^d a^i_{\mathbf{k}} \mathbf{x}^\mathbf{k}, \qquad i=1, \ldots, n,
 \label{eq:homd}$$ of degree $d$. We assume that this polynomial field is *diagonal* in the sense described earlier. Diagonality implies some (mild) restrictions on the nonzero coefficients of , namely that, for a given index $\mathbf{k}$, one has $ a^i_{\mathbf{k}} \not=0$ provided that $k_i\not=0$.

Rather than looking at the vector fields, it is useful to focus on the divergence polynomial and reconstruct the vector field from the divergence.

Consider the (diagonal) polynomial vector field . Assume that the vector field is divergence-free. Let $$P_{d-1}(\mathbf{x})= \sum_{|\mathbf{j}|=0}^{d-1} p_{\mathbf{j}}{\mathbf{x}^\mathbf{j}}.$$ be the divergence of . For each multi-index $\mathbf{j}$, consider the polynomial field $$\dot{x}_i = a^i_{\mathbf{j}+\mathbf{e}_i} x_i \mathbf{x}^\mathbf{j}, \qquad 1\leq i \leq n.
 \label{eq:xi}$$ where $\mathbf{e}_i$ the canonical unit vector in $\RR^n$ with $1$ in the $i$th position and $0$ elsewhere. The polynomial field is divergence free. Moreover, can be uniquely split in sum of vector fields of the form .

Since monomials are basis for the set of polynomials of degree $d$, zero-divergence of implies that $$p_{\mathbf{j}} = 0 \qquad |\mathbf{j}|=0, \ldots, d-1,$$ where $p_{\mathbf{j}}$ is the coefficient of $\mathbf{x}^\mathbf{j} = x_1^{j_1}\cdots x_n^{j_n}$ in the divergence polynomial. For any multi-index $\mathbf{j}$, a contribution to the coefficient $p_{\mathbf{j}}$ comes from $\partial_{x_i}f_i(\mathbf{x})$, for each $1\leq i\leq n$. If the $f_i(x)$ are polynomials, such a term originates from derivation of $x_i \mathbf{x}^\mathbf{j}=\mathbf{x}^{\mathbf{j+e}_{i}}$ with respect to $x_i$. This corresponds to the term $a^i_{ \mathbf{j}+\mathbf{e}_i} x_i \mathbf{x}^\mathbf{j}$ in $f_i(\mathbf{x})$. This procedure picks up all the terms in $f_i(\mathbf{x})$ contributing to $p_{\mathbf{j}}$ and the condition $p_{\mathbf{j}}=0$ guarantees that is divergence free. The uniqueness comes from the fact that $\int \mathbf{x}^\mathbf{j} \D x_i = \frac{1}{j_i+1}\mathbf{x}^{\mathbf{j+e}_{i}} + c h(\mathbf{x})$, where $\partial_{x_i} h(\mathbf{x})=0 $, i.e. $h$ is an arbitrary function of the remaining variables and does not contribute to the diagonal part of the vector field.

In the sequel, we will often use the short-hand notation $$\dot{\mathbf{x}} = F_\mathbf{j} (\mathbf{x})
    \label{eq:Fj}$$ for . Each of these divergence-free vector fields is associated to a monomial basis element, and will be called an *elementary divergence-free vector field* (in short, EDFVF).

Properties of the EDFVF {#sec:integrals}
=======================

Without loss of generality, we can consider a EDFVF $F_\mathbf{j}$ of the form $$\dot{x}_i =  a_i x_i \mathbf{x^j}, \qquad i=1, \ldots, n,
    \label{eq:edfvf}$$ (for simplicity, we have dropped the dependence of $\mathbf{a}$ on the actual index $\mathbf{j}$). The divergence free condition amounts to the algebraic relation $$\mathbf{a}^T (\mathbf{j+1}) =0, \qquad \mathbf{1} = (1, \ldots, 1)^T.
    \label{eq:edfcond}$$

\[th:integrability\] The function $I_1 = \mathbf{x^{j+1}}$ is an integral of the EDFVF . Moreover, if $\mathbf{b}\in \mathbb{R}^n$ is any other vector orthogonal to $\mathbf{a}$, then $\mathbf{x^b}$ is also an integral of . In particular, it follows that has $n-1$ independent integrals of motion, hence the EDFVF is integrable.

Let $\mathbf{b}$ be an arbitrary vector orthogonal to $\mathbf{a}$, i.e. such that $\mathbf{a}^T \mathbf{b}=0$ (the case $\mathbf{b}_1 = \mathbf{j+1}$ is just a particular choice of $\mathbf{b}$). We have $$\begin{aligned}
    \frac{\D}{\D t} \mathbf{x^b} &=& \sum_{i=1}^n b_i \dot{x}_i \mathbf{x^{b-e_i}} = \sum_{i=1}^n a_i b_i x_i \mathbf{x^j}\mathbf{x^{b-e_i}} \\
    &=& \sum_{i=1}^n a_i b_i \mathbf{x^{j+b}} = (\mathbf{a}^T \mathbf{b}) \mathbf{x^{j+b}} =0,\end{aligned}$$ because of the orthogonality condition. Furthermore, since $\mathbf{a} \in \RR^n$, we can find further $n-2$ vectors orthogonal to $\mathbf{a}$ and to $\mathbf{j+1}=\mathbf{b}_1$, i.e. vectors $\mathbf{b}_2, \ldots, \mathbf{b}_{n-1}$ (each corresponding to the integral $I_k= \mathbf{x}^{\mathbf{b}_k}\not=0$), such that $\{ \mathbf{a}, \mathbf{b}_1, \mathbf{b}_2, \ldots, \mathbf{b}_{n-1}\}$ is an orthogonal basis of $\RR^n$. These integrals are functionally independent, i.e. any linear combination $$c_1 \nabla I_1 + \ldots + c_{n-1} \nabla I_{n-1} = 0 
    \label{eq:lindep}$$ has only the trivial solution $$c_1 = c_2 =\ldots = c_{n-1} =0.$$ To show this, note that $$\nabla I_k = \left[ \begin{array}{c} (\mathbf{b}_k)_1 \mathbf{x}^{\mathbf{b}_k -\mathbf{e}_1}\\
    \vdots \\
    (\mathbf{b}_k)_n \mathbf{x}^{\mathbf{b}_k -\mathbf{e}_n}\end{array}\right]
    =
    \left[ \begin{array}{c} (\mathbf{b}_k)_1/x_1 \\
    \vdots \\
    (\mathbf{b}_k)_n/x_n \end{array}\right] I_k.$$ Because of the arbitrariness of $x_1, \ldots, x_n$, condition is equivalent to the linear system $$[\mathbf{b}_1 , \ldots, \mathbf{b}_{n-1}] \mathbf{c} = \mathbf{0}$$ $\mathbf{c}=(c_1, \ldots, c_{n-1})^T$, which admits only the trivial solution as the matrix columns are orthogonal, hence linearly independent. In conclusion, each of the vector $\mathbf{b}_k$, $k=1,\ldots, n-1$, will give rise to an independent integral of motion, hence integrability of the vector field.

![A 3-dimensional system with $\mathbf{j} = (1,1,1)^T$, $\mathbf{a}=(-\frac53, \frac43,\frac13)^T$ and $\mathbf{b}_2 = \mathbf{a} \times (\mathbf{j+1})$ (cross product). The solution of the system lies on the intersection of the surfaces $I_1(t) = \mathbf{x^{j+1}}$ (lighter surface) and $I_2(t) = \mathbf{x}^{\mathbf{b}_2}$ (darker surface).[]{data-label="fig:int"}](integrals.pdf){width="7cm"}

Figure \[fig:int\] shows two such integral surfaces for a 3-dimensional system with $\mathbf{a}=(-\frac53, \frac43,\frac13)^T$, $\mathbf{j} = (1,1,1)^T$ and initial condition $(1,1,1)^T$. The solution (blue thick line) lies at the intersection of the two integrals of motion, $I_1$ (yellow, lighter color) and $I_2$ (purple, darker color), the latter corresponding to $\mathbf{b}_2 = \mathbf{a} \times (\mathbf{j+1})$.

Among the $n-1$ integrals, the first one $I_1 = \mathbf{x^{j+1}}$ is more fundamental than the others, because it is a direct consequence of the divergence-free condition. In particular, it can be always used to transform the EDFVF (which has nonzero diagonal part) into a new system with zero diagonal part.

The EDFVF is equivalent to the off-diagonal divergence free system $$\dot{x}_i =  a_i I_1(t_0) \mathbf{x^{-1+e_i}}, \qquad i=1, \ldots, n,$$ where $I_1(t) = \mathbf{x^{j+1}}$ is an integral of the system.

The above system is obviously off diagonal since the $i$th component of $\mathbf{-1+e}_i$ is zero.

We have demonstrated integrability of the polynomial EDFVF by providing the existence of $n-1$ independent integral of motions. However, it is well known that the explicit knowledge of the integrals of motions does not necessarily give explicit informations about the solution. Below, we give the explicit solution for the EDFVF . This is based on two fundamental observations:

-   $x_i(t)$ can be reconstructed by integration provided that $\mathbf{x^j}(t)$ is known,

-   remarkably, the basis functions $\mathbf{x^j}$ *always* obey the *same type of solvable* differential equation, independently of the multi-index $\mathbf{j}$.

\[th:dcase\] Consider the EDFVF associated to the monomial basis element $\mathbf{x^j}$. The associate basis function $\mathbf{x}^\mathbf{j}(t)$ obeys the differential equation $$\frac{\mathrm{d}}{\mathrm{d} t} \mathbf{x}^\mathbf{j} = c_\mathbf{j} \mathbf{x}^{2\mathbf{j}},
    \label{eq:xjdot}$$ where $c_\mathbf{j}$ is the constant $$c_\mathbf{j} = \mathbf{a}^T\mathbf{j}.
    \label{eq:kj}$$ For sufficiently small $\Delta t\leq \displaystyle{\frac{1}{|c_\mathbf{j}\mathbf{x}^\mathbf{j} (t_0)|}}$, the solution of is $$x_i(t) = \begin{cases} x_i(t_0) \Big(1-c_\mathbf{j} \mathbf{x}^\mathbf{j}(t_0)(t-t_0)\Big)^{-r^i_{\mathbf{j}}} & \hbox{for $c_\mathbf{j}\not=0$}\\
    x_i(t_0) \E^{a_i \mathbf{x^j}(t_0) (t-t_0)} & \hbox{for $c_\mathbf{j}=0$},
    \end{cases}
    \label{eq:xi_sol}$$ for $i=1, \ldots, n$, where $$r^i_{\mathbf{j}} = \frac{a_i}{c_{\mathbf{j}}}.
\label{eq:rij}$$

By direct computation, we have: $$\begin{aligned}
 \frac{\mathrm{d}}{\mathrm{d} t} \mathbf{x}^\mathbf{j} &=& j_1 \dot{x}_{1} \mathbf{x}^{\mathbf{j} -\mathbf{e}_1} + \cdots + j_n \dot{x}_{n} \mathbf{x}^{\mathbf{j} -\mathbf{e}_n}\\
  &=& j_1 a_1 x_{1} \mathbf{x}^\mathbf{j}   \mathbf{x}^{\mathbf{j} -\mathbf{e}_1}  + \cdots + j_n a_n x_{n} \mathbf{x}^\mathbf{j}   \mathbf{x}^{\mathbf{j} -\mathbf{e}_n} = \sum_{i=1}^n j_{i} a_i \mathbf{x}^{2\mathbf{j}}\\
  &=& c_\mathbf{j}\mathbf{x}^{2\mathbf{j}}.\end{aligned}$$

The second passage follows from substituting the components of the vector field $F_{\mathbf{j}}$ in place of $\dot{x}_{i}$ and the last passage follows from the definition of $c_\mathbf{j}$. Remarkably, this differential equation is exactly the same as the first one in the quadratic case [@mclachlan09evp], therefore it can be solved in the same manner, for instance, by separation of variables, yielding $$\mathbf{x}^\mathbf{j}(t) = \frac{\mathbf{x}^\mathbf{j}(t_0)}{1-c_\mathbf{j} \mathbf{x}^\mathbf{j}(t_0)(t-t_0)}
    = {\mathbf{x}^\mathbf{j}(t_0)}\Big( 1-c_\mathbf{j} \mathbf{x}^\mathbf{j}(t_0)(t-t_0)\Big) ^{-1},
    \label{eq:xj}$$ in the interval $[t_0, t_0 + \Delta t]$.

Thereafter, each of the equations is of the form $\dot{x}_i = h(t) x_i(t)$, where $h(t) = a_i \mathbf{x}^\mathbf{j}(t)$ is a known function. As for the quadratic case (see [@mclachlan09evp]), the equation is linear in $x_i$, with variable coefficients, and has solution $x_i(t) = x_i(t_0) \E^{\int_{t_0}^t h(\tau) \mathrm{d} \tau}$. Hence, for $c_\mathbf{j}\not=0$, $$\begin{aligned}
    x_i(t) &=& x_i(t_0) \E^{-\frac{a_i}{c_\mathbf{j}} \ln (1-{c_\mathbf{j}}\mathbf{x}_\mathbf{j}(t_0) (t-t_0))} = x_i(t_0) \Big(\E^{\ln (1-{c_\mathbf{j}}\mathbf{x}_\mathbf{j}(t_0) (t-t_0))}\Big)^{-r^i_{ \mathbf{j}}} \\
    &=& x_i(t_0) \Big(1-c_\mathbf{j} \mathbf{x}^\mathbf{j}(t_0)(t-t_0)\Big)^{-r^i_{\mathbf{j}}},\end{aligned}$$ where $r^i_{\mathbf{j}}$ is defined as in . The case $c_\mathbf{j}=0$ follows in a similar manner. This completes the proof of the statement.

Consider the divergence-free differential equation $$\begin{array}{lcl}
    \dot x_1 &=& x_1^4 + x_3\sin x_1 \cos x_2 + 4 x_2 x_3\\
    \dot x_2 &=& -\frac12 x_1^3 x_2-\frac12 x_3 \cos x_1 \sin x_2 -x_1 \\
    \dot x_3 &=& -\frac12 x_1^3 x_3-\frac14 x_3^2 \cos x_1 \cos x_2 + x_1^2 x_2.
    \end{array}
    \label{eq:poly_ex}$$ This system consists of a diagonal part and a off-diagonal part. In many applications, it is common to approximate such vector fields by a polynomial field of low degree (like linear, quadratic, etc.). This is typically done by using a Taylor expansion and truncating to a certain degree $d$. The choice $d=5$ gives the system $$\dot{\mathbf{x}} = \mathbf{f}^\mathrm{diag} (\mathbf{x}) + \mathbf{f}^\mathrm{offdiag}(\mathbf{x}),
    \label{eq:poly_ex_split}$$ where $$\mathbf{f}^{\mathrm{diag}} (\mathbf{x}) = \left[ \begin{array}{c}
    x_1 x_3 + x_1^4 -\frac16 x_1^3 x_3 -\frac12 x_1 x_2^2 x_3\\
    -\frac12 x_2 x_3-\frac12 x_1^3 x_2 +\frac1{12} x_2^3 x_3 +\frac14 x_1^2 x_2 x_3\\
    -\frac14 x_3^2-\frac12 x_1^3 x_3 +\frac18 x_1^2 x_3^2 +\frac18 x_2^2 x_3^2
    \end{array}
    \right],
    \quad
    \mathbf{f}^\mathrm{offdiag}(\mathbf{x}) = \left[ \begin{array}{c}
        4 x_2 x_3 \\
        -x_1 \\
        x_1^2 x_2
    \end{array}
    \right].$$ The off-diagonal part can be treated using canonical shears. For the diagonal, observe that the divergence-free condition for is a polynomial in $x_3, x_1^3, x_1^2 x_3, x_2^2 x_3$. These correspond to the multi-indices $\mathbf{j}=(0,0,1)$, $\mathbf{k} = (3,0,0)$, $\mathbf{l} = (2,0,1)$, and $\mathbf{m} = (0,2,1)$. Each of these multi-indices is associated to an EDFVF. For instance, for $\mathbf{l}$, one has $c_{\mathbf{l}} =  -\frac5{24}$ and gives $r^1_{ \mathbf{l}} = \frac45, r^2_{\mathbf{l}} = -\frac65,  r^3_{\mathbf{l}} =-\frac35$. Assuming that the solution is known at $t_k$, we obtain $$\begin{aligned}
    x_1(t) &=& x_1(t_k) (1+ \frac{5}{24} x_1^2(t_k)x_3(t_k) (t-t_k))^{-4/5}, \\
    x_2(t) &=& x_2(t_k) (1+ \frac{5}{24} x_1^2(t_k)x_3(t_k) (t-t_k))^{6/5},\\
    x_3(t) &=& x_3(t_k) (1+ \frac{5}{24} x_1^2(t_k)x_3(t_k) (t-t_k))^{3/5},\end{aligned}$$ for $t \in [t_k, t_{k+1}]$. A similar treatment applies to the EDFVF associated to the other multi-indices.

Extensions to the case of negative (Laurent polynomials), rational and real multi-indices bases {#sec:laurent}
===============================================================================================

In this section, we investigate the extension of the theory to treat the case of splitting in terms of the form $\mathbf{x^j}$ with $\mathbf{j}\in \mathbb{Z}^n$, or, more generally, $\mathbf{j} \in \mathbb{R}^n$.

We assume that the divergence of the vector field can be expressed in the form $$P(x_1, x_2, \ldots, x_n) = \sum_{\mathbf{j}\in \mathcal{J}}p_\mathbf{j} \mathbf{x}^\mathbf{j},
    \label{eq:divLaurent}$$ where the $\mathcal{J}\subset \mathbb{F}^n$ (where $\mathbb{F}= \mathbb{Z},\mathbb{R}$) is discrete finite set of multi-indices.

\[th:jine-1\] Assume that the divergence free vector field has components $$f_i (\mathbf{x}) = \sum_{\mathbf{l}\in \mathcal{L}} a^i_{\mathbf{l}} \mathbf{x}^{\mathbf{l}}, 
    \label{eq:gen_mon}$$ where $\mathcal{L} \subset \mathbb{F}^n$ (discrete set), $i=1,\ldots, n$. Then, the divergence of $\mathbf{f}$ is a generalized polynomial of the form . Moreover, either the indices $\mathbf{j}$ in have components $j_i \not=-1$, or, if some $j_i = -1$, then equation $i$ does not contribute to this multi-index.

Terms contributing to the divergence originate from $\partial f_i /\partial x_i$. Any $\mathbf{x^l} = x_1^{l_{1}}\cdots x_i^{l_{i}} \cdots x_i^{l_{n}}$ contributes $a^i_{\mathbf{l}}  l_{i} x_1^{l_{1}}\cdots x_i^{l_{i}-1} \cdots x_n^{l_{n}}= a^i_{\mathbf{l}}  l_i \mathbf{x^j}$, to , $\mathbf{j = l-e}_i$, and it is easily seen that linear combinations of these terms yield a “generalized” polynomial of the form . As far as the second part of the statement is concened, we see that $j_{i}=-1$ (yielding a term of the type $x_i^{-1}$ in the divergence) implies $l_i=0$ (this contribution is always identically zero). This rules out that a contribution to the term $x_i^{-1}$ can come from the $i$th equation but it does not rule out such a contribution to come from another equation.

As a matter of facts, a contribution to the term $x_1^{j_{1}} \cdots x_i^{-1}\cdots x_n^{j_{n}}$ from equation $i$ can occur if and only if a term of the type $x_1^{j_{1}} \cdots \log x_i \cdots x_n^{j_{n}}$ appears in $f_i$, which is ruled out from the hypotesis that $f_i$ is of a Laurent polynomial type.

The results in Theorem \[th:integrability\]-\[th:dcase\] can now be extended with minor modifications to treat the generalized polynomial case.

\[th:Laurent\] Let be given the differential equation , where $\mathbf{f}$ component-wise satisfies . Let $$0 = \sum_{\mathbf{j}\in\mathcal{J}} p_\mathbf{j} \mathbf{x}^\mathbf{j}
    \label{eq:divfree}$$ be the corresponding divergence-free condition. For each multi-index $\mathbf{j} \in \mathcal{J}$, there exists the associated elementary divergence-free vector field $$\dot x_i = (1-\delta_{i,l}) a^i_{ \mathbf{j+e}_i} x_i \mathbf{x^j}, \qquad i=1,\ldots, n, \qquad j_{l}=-1,
    \label{eq:xidotdelta}$$ where $\delta_{i,l}$ is the Kronecker delta and $l$ is any index such that $j_l=-1$. The elementary divergence free vector fields can be written in the form $$\dot x_i = a_i x_i \mathbf{x^j},$$ (disregarding the dependance on the index $\mathbf{j})$, they are volume preserving by construction and possess $n-1$ integrals of motions, hence they are integrable and their solution is explicitly given as in \[th:dcase\].

The proof is very similar to the monomial case. The only difference is the presence of the term $\delta_{i,l}$, taking care of excluding the contribution of the $l$th coefficient in $a_{ \mathbf{j} +\mathbf{e}_{i}}$ (because this contributes to a off-diagonal term, see above lemma) for the specific $\mathbf{j}$ with $j_l = -1$.

Numerical examples {#sec:num}
==================

We present some interesting applications to volume-preserving vector fields for a quadratic and a cubic Stokes flow. In our literature search, we didn’t find any relevant examples of Laurent volume-preserving fields, and we will test the methods on a artificial example, for the sake of illustration.

Quadratic Stokes flow
---------------------

We consider a quadratic volume-preserving system introduced in [@MR2433381] to study the distruction of adiabatic invariance under separatrix crossing, $$\arraycolsep=1.4pt
  \renewcommand{\arraystretch}{1.4}
  \begin{array}{rcl}
    \dot x_1 &=& -8x_1 x_2 + \varepsilon x_3,\\
    \dot x_2 &=& 11 x_1^2 +3x_2^2 + x_3^2 - 3,\\
    \dot x_3 &=& 2x_3x_2 - \varepsilon x_1.
\end{array}
\label{eq:QuadStokes}$$ The system is integrable for $\varepsilon=0$. When $\varepsilon \not=0$ the system is no longer integrable but solutions inside the unit sphere remain bounded in the sphere. We consider the non-integrable case and split the system into a diagonal and off-diagonal part, $$\mathbf{f}^\mathrm{diag}: \quad \begin{array}{lcl}
    \dot x_1 &=& -8x_1x_2 \\
    \dot x_2 &=& 3x_2^2 \\
    \dot x_3 &=& 2x_3x_2  \end{array}
    \hspace*{30pt}
    \mathbf{f}^\mathrm{offdiag}: \quad \begin{array}{lcl}
    \dot x_1 &=& \varepsilon x_3,\\
    \dot x_2 &=& 11 x_1^2 + x_3^2 - 3,\\
    \dot x_3 &=& - \varepsilon x_1. \end{array}$$ The diagonal part corresponds to the index $\mathbf{j} = (0,1,0)$, which is solved by a single flow as in Prop. \[th:dcase\], with $c_\mathbf{j} = [0,1,0] [-8,3,2]^T=3$ (scalar product) and $r_\mathbf{j} = \frac13 \times [-8,3,2]^T$.

![The quadratic Stokes flow with $\varepsilon=0.1$. Initial condition $\mathbf{x}_0=[ 0,0,0.96]^T$, integration time $T=500$. Top left: volume-preserving method, order 2, implemented with stepsize $h=0.01$. Top right: `ode45`, using step size control, with options `RelTol=1e-6` (default value `1e-3`). Bottom left: `ode45`, using step size control, default implementation, until $T=220$, the method becomes unstable at $T=2.2887$e+2. Bottom right: same as bottom left, but letting $T=250$.[]{data-label="fig:QuadStokes"}](fig2.pdf){width=".9\textwidth"}

![Histogram of the step sizes used by `ode45` to meet the prescribed tolerance for the quadratic volume-preserving flow . Dark gray: `RelTol=1e-6`. Light gray: `RelTol=1e-3`.[]{data-label="fig:hist45"}](hist2.pdf){width=".95\textwidth"}

In Figure \[fig:QuadStokes\] we show the numerical trajectories corresponding to $\varepsilon = 0.1$ and initial condition $\mathbf{x}_0=[ 0,0,0.96]^T$, for $T=500$. The second-order volume-preserving method (top left) is implemented as described in Algorithm 1, with a constant step size $h=0.01$, and gives a nice and bounded trajectory. The first-order method gives very similar results. For comparison, we use Matlab’s `ode45` method, that is explicit, fourth-order, and uses step size control. The `ode45` method is not volume preserving unless the solution is computed to machine accuracy. The standard implementation of `ode45` (with stepsize control) becomes unstable at $T=2.2887$e+2 (see illustration bottom right). In the bottom left picture we display the integration up to $T=220$. To obtain a result similar to the volume-preserving method, it is necessary to decrease the error tolerance. The top right plot illustrates the results by `ode45` with options set to `RelTol=1e-6` (default value `1e-3`).

Figure \[fig:hist45\] shows the histogram for the stepsize chosen by `ode45` to meet the required tolerance. For the smallest error tolerance, the average step size is around $h=0.01$, which justifies the choice for our explicit volume-preserving methods. Numerical tests revealed that our methods were stable up to a choice $h \approx 0.05$. We performed also long time integration, $T=100000$, with $h=0.05$, and the solution still stayed bounded for the volume-preserving method.

As far as computational time is concerned, our methods have the advantage of being explicit. The speedup, with respect to `ode45` implemented with the option `RelTol=1e-6` for stability, is approximatly $3.2$ for the first-order volume-preserving method and $3.15$ for the second-order method, over an average of 100 runs with the same initial condition as above, see Table \[tab:1\].

                       `ode45`   vol. pres., order 1   vol. pres., order 2
  ------------------ --------- --------------------- ---------------------
  cpu time \[sec\]      3.8687                1.2099                1.2260
  speed up                   1                3.1976                3.1557

  : CPU time (in seconds) for `ode45` with `RelTol=1e-6` and the explicit volume-preserving splitting methods. The values, for indication only, are computed as an average of 100 runs on an standard laptop (MacBook Pro) with initial condition $\mathbf{x}_0=[ 0,0,0.96]^T$ and $T=500$.[]{data-label="tab:1"}

![The cubic Stokes flow with $\mathbf{w}=\mathbf{0}$ (integrable case). Simulation of trajectories for $T=20000$ for various initial conditions (the circles on the line denote the starting value). Top left: volume-preserving method, order 2, implemented with stepsize $h=0.01$. Top right: volume-preserving method, order 1, implemented with stepsize $h=0.01$. Bottom left: `ode45`, using step size control, with `RelTol=1e-6`. Bottom right: `ode45`, default implementation (`RelTol=1e-3`). []{data-label="fig:CubicStokes1"}](fig4.pdf){width=".9\textwidth"}

Cubic Stokes flow
-----------------

We consider a cubic volume-preserving flow which has been used for the investigation of the kinematics of bounded steady Stokes flows, in particular, to investigate the streamlines inside a neutrally buoyant spherical drop immersed in a genera linear flow [@stone91csi]. This system is described by the divergence-free differential equation $$\dot{\mathbf{x}} = \frac12 [(5 r^2 -3) E \mathbf{x} - 2\mathbf{x} (\mathbf{x}^T E \mathbf{x})] + \frac12 \mathbf{w} \times \mathbf{x},$$ where $r^2 = \mathbf{x}^T\mathbf{x}$, $\mathbf{w}$ is the vorticity vector and $E$ is the symmetric traceless rate-of-strain tensor of the external motion, that, in dimensionless terms, takes the form $$E = \left( \begin{array}{ccc} 1/(1+\alpha) &0 & 0 \\
        0 & \alpha/(1+\alpha) &0  \\
        0 & 0 & -1\end{array} \right),$$ where $\alpha = E_{22}/E_{11}$, and $\times$ is the usual cross product.

In Cartesian coordinates, the dynamical system for investigating fluid particle motion takes the dimensionless form: $$\arraycolsep=1.4pt
  \renewcommand{\arraystretch}{1.4}
  \begin{array}{rcl}
 \dot x_1 &=& \frac12 \left[ (5r^2 -3) \frac{x_1}{1+\alpha} -2 x_1 \left( \frac{x_1^2}{1+\alpha} + \frac{\alpha x_2^2}{1+\alpha} - x_3^2 \right) \right] + \frac12 (w_2 x_3 - w_3 x_2) \\
 \dot x_2 &=& \frac12 \left[ (5r^2 -3) \frac{\alpha x_2}{1+\alpha} -2 x_2\left( \frac{x_1^2}{1+\alpha} + \frac{\alpha x_2^2}{1+\alpha} - x_3^2 \right) \right] + \frac12 (w_3 x_1 - w_1 x_3) \\
 \dot x_3 &=& \frac12 \left[ -(5r^2 -3)x_3 -2 x_3 \left( \frac{x_1^2}{1+\alpha} + \frac{\alpha x_2^2}{1+\alpha} - x_3^2 \right) \right] + \frac12 (w_1 x_2 - w_2 x_1). \\
\end{array}
\label{eq:CubicStokes}$$ The case $\mathbf{w=0}$ is integrable. The generic case $\mathbf{w \not=0}$ is not integrable, but also in this case, solutions with initial condition in the unit sphere are bounded to the unit sphere, see [@stone91csi] for more details.

In the following experiments, the cubic Stokes flow was computed numerically by volume-preserving methods and Matlab’s `ode45` method. Figure \[fig:CubicStokes1\] shows the comparison of the trajectories by the volume-preserving methods (order $1$ and $2$) and Matlab’s `ode45` method with options `RelTol=1e-6` and `RelTol=1e-3` for $\mathbf{w}=\mathbf{0}$. The volume-preserving methods preserve the qualitative behaviour for various initial conditions (denoted by small circles) in Figure \[fig:CubicStokes1\]. However, Matlab’s `ode45` method, with standard options control, has a dissipative behaviour. To obtain the same visual result as the volume-preserving methods, one has to reduce the tolerance on the error, for instance set `RelTol=1e-6`. See the picture at the bottom of Figure \[fig:CubicStokes1\] for more details.

In the numerical investigation presented below (see Figure \[fig:CubicStokes2\]), we study the effect of changing both the orientation and the magnititude of $\mathbf{w}$. Without loss of generality, we will take $\mathbf{w}=(w_1,0,w_3)$. The orientation of $\mathbf{w}$, measured from the $x_3$-axis in the $(x_1,x_3)$ plane, will be denoted by $\varTheta$. To describe the three-dimensional particle paths, we present Poincar[é]{} sections through the $(x_1,x_3)$ plane. We compute the particle paths by solving numerically using a volume-preserving method (order 2), Matlab’s `ode45` method with options `RelTol=1e-6` and standard setup (`RelTol=1e-3`) respectively. For more detail about the experiments design, see [@stone91csi] for reference.

To start with, we examine the case $\varTheta=0.275\pi$ and $\|\mathbf{w}\|=1.5$ illustrated in the first row of Figure \[fig:CubicStokes2\]. The left column is the Poincaré section obtained by the volume-preserving method (order 2), while the middle and right columns are obtained by Matlab’s `ode45` method with options `RelTol=1e-6` and standard implementation (`RelTol=1e-3`) respectively. The particle parth Poincaré section, computed by our volume-preserving method and Matlab’s `ode45` method with options `RelTol=1e-6` are visually identical to the figure showed in [@stone91csi] (see Figure 3 in [@stone91csi] for detail), while Matlab’s `ode45` method with standard options `RelTol=1e-3` gives a totally wrong section.

We next study the Poincaré section for a fixed orientation of the vorticity vector, but increasing magnitude of the vorticity. The second and third rows of Figure \[fig:CubicStokes2\] show simulations for an orientation $0.2\pi$ from the $x_3$-axis by increasing $\|\mathbf{w}\|$ from $2.5$ to $4.0$. As the magnitude of $\mathbf{w}$ increases, islands are created and destroyed until $\|\mathbf{w}\|=1.4$. For larger magnitudes of the vorticity vector, the streamlines structure is again modified and becomes much more complex [@stone91csi].

The last two rows in the Figure \[fig:CubicStokes2\] show the cases where the vorticity orientation increases from $0.02\pi$ to $0.4\pi$ with $\|\mathbf{w}\|=2.0$ fixed.

For all the experiments illustrated in Figure \[fig:CubicStokes2\], the simulations have been performed with $T=20000$ and initial condition $\mathbf{x}_0= [-0.1689,0,-0.0437]^T$. The volume-preserving method (order 2) is implemented with constant step size $h=0.01$ (the last row with step size $h=0.005$). Each row corresponds to different value of $\mathbf{w}$ as explained above.

From Figure \[fig:CubicStokes2\], we can immediately observe that the particle parths in the Poincaré sections by the volume-preserving method behave qualitatively similarly to the exact solution, while, for non volume-preserving methods (like Matlab’s `ode45`), stronger restrictions on the step size are required to obtain the desired dynamic result.

![The cubic Stokes flow with $\mathbf{w}\not=0$ (non-integrable case). Simulation of trajectories for $T=20000$ with initial condition $\mathbf{x}_0= [-0.1689,0,-0.0437]^T$. Left column: volume-preserving method, order 2, stepsize $h=0.01$ (with exception of last row, for which $h=0.005$). Middle column: `ode45` with option `RelTol=1e-6`. Right column: `ode45`, standard implementation (`RelTol=1e-3`). Each row corresponds to a different value of $\mathbf{w}=(w_x,0,w_y)$. From top to bottom: $\|\mathbf{w}\|=1.5, \varTheta=0.275\pi$; $\|\mathbf{w}\|=2.5, \varTheta=0.2\pi$; $\|\mathbf{w}\|=4.0, \varTheta=0.2\pi$; $\|\mathbf{w}\|=2.0, \varTheta=0.02\pi$; $\|\mathbf{w}\|=2.0, \varTheta=0.4\pi$.[]{data-label="fig:CubicStokes2"}](fig5.pdf){width="\textwidth"}

Laurent polynomial example
--------------------------

In order to illustrate our method for Laurent-type volume-preserving fields, we consider the following example: $$\arraycolsep=1.4pt
  \renewcommand{\arraystretch}{1.4}
  \begin{array}{rcl}
    \dot x_1 &=& 3x_1^{-2}x_2^2 + 2x_1^3x_2^{-3},\\
    \dot x_2 &=& 2x_1^{-3}x_2^3 +3x_1^2x_2^{-2}.
\end{array}
\label{eq:Laurent}$$ According to the approach presented in this paper, we split the above equation (\[eq:Laurent\]) into two vector fields $F_{1}$ and $F_{2}$, corresponding to the multi-index $\mathbf{j}_{1}=(-3,2)$ and $\mathbf{j}_{2}=(2,-3)$, $$F_{1} : \quad \begin{array}{lcl}
    \dot x_1 &=& 3x_1^{-2}x_2^2 \\
    \dot x_2 &=& 2x_1^{-3}x_2^3
        \end{array}
    \hspace*{30pt}
    F_{2} : \quad \begin{array}{lcl}
    \dot x_1 &=& 2x_1^3x_2^{-3}\\
    \dot x_2 &=& 3x_1^2x_2^{-2}.
    \end{array}$$

We can integrate the vector fields $F_{1}$ and $F_{2}$ explicitly according to Theorem \[th:Laurent\]. By choosing different initial values around $[-0.5689,0.0437]$, we compare the results obtained by our method to Matlab’s `ode45` method with options `RelTol=1e-8`, `RelTol=1e-8`, `RelTol=1e-10` and `RelTol=1e-6` and standard implentation `RelTol=1e-3` respectively.

![Laurent polynomials problem (\[eq:Laurent\]) with different initial values. Left: volume-preserving method, order 2, implemented with stepsize $h=0.001$. Middle: `ode45`, using stepsize control, with options `RelTol=1e-8`, `RelTol=1e-8`, `RelTol=1e-10` and `RelTol=1e-6`. Right: `ode45`, using step size control, with standard options `RelTol=1e-3`.[]{data-label="fig:L1"}](fig6.pdf){width=".9\textwidth"}

In Figure \[fig:L1\], we show the numerical results for the volume-preserving method and Matlab’s `ode45` method, for initial values $[-0.5689+k\Delta,0.0437+l\Delta]^T$, $k,l=0,1$, $\Delta = 0.02$ (a small square). The standard Matlab’s `ode45` implementation diverged for early values of $T$ (right column). Thus, to obtain the correct dynamical behaviour, stronger restrictions on the error had to be imposed (middle column): a relative tolerance `RelTol=1e-8` gave results similar to the volume-preserving method for two of the initial conditions, for another one, the error control needed to be sharpened (`RelTol=1e-10`), while for the last one, it was sufficient to use `RelTol=1e-6`.

![Laurent polynomials problem (\[eq:Laurent\]) with initial conditions $(0.0437,-0.5489)$, $(-0.5489,0.437)$ in the different Octant in the plane. Left: volume-preserving method, order 2, implemented with stepsize $h=0.001$. Middle: `ode45`, using stepsize control, with options `RelTol=1e-8`. Right: `ode45`, using stepsize control, with standard implementation `RelTol=1e-3`.[]{data-label="fig:LL1"}](fig7.pdf){width=".9\textwidth"}

![Histogram of the stepsizes used by `ode45` to meet the prescribed tolerance for in $[0,10]$. Gray: `RelTol=1e-8`. Black: `RelTol=1e-3`[]{data-label="fig:L2"}](laa2){width="95.00000%"}

Figure $\ref{fig:LL1}$ shows the trajectories obtained with initial values $[0.0437,-0.5489]^T$ and $[-0.5489,0.437]^T$ (opposite octants) using the volume-preserving method (order 2) with stepsize $h=0.001$ and Matlab’s `ode45` methods with options `RelTol=1e-8` and standard implementation. Matlab’s `ode45` method with standard implementation was again unstable, and the trajectories diverged before the end of the integration interval was reached. The trajectories of the system are very sensitive with respect to the initial condition (ill conditioned problem) and small changes of initial values may cause tremendous difference in the actual trajectories. This type of problems is well known to affect the numerical methods by imposing constraints on the step size (small step size must be used to remain sufficiently close to the true trajectory). Also our volume-preserving methods were affected by this, as well as Matlab’s `ode45`, but they required milder restrictions.

Figure \[fig:L2\] shows the histogram for the stepsize chosen by `ode45` with `RelTol=1e-8` and standard implementation (`RelTol=1e-3`) respectively, to meet the required tolerance in $[0,10]$. For the case `ode45` with options `RelTol=1e-8`, the average stepsize is the interval $(0,0.005)$. Looser restrictions gave diverging solutions. For comparison, the volume-preserving methods could use $h=0.25$ and still remain bounded. However, for longer time integrations, also the volume-preserving methods had to reduce the step size to avoid divergence.

Conclusion and further remarks {#sec:conclusion}
==============================

In this paper we have presented new, explicit, volume-preserving splitting methods for Laurent polynomial divergence-free vector fields. The methods rely on a decomposition of the divergence into a monomial basis. For each monomial basis element, a corresponding divergence-free elementary vector field is identified and integrated exactly (hence the solution is volume preserving). These elementary fields are composed by a splitting method, giving rise to a overall explicit volume-preserving method.

The method proposed are a significant contribution to the understanding of volume-preserving methods. To our knowledge, there is no explicit volume-preserving method that can deal with arbitrary polynomial vector fields of arbitrary degree, with exception of the case of linear and quadratic vector fields.

As far as stability is concerned, our methods were stable for sufficiently small step size, a feature that is common to explicit methods. An upper bound on the step size can be obtained by Prop \[th:dcase\]. Beware that this bound can be unrealistic for long time numerical integration and further step size restriction might occur, especially if the solution is unbounded, which is often the case of divergence-free vector fields and and high degrees of nonlinearity. These problems are harder to integrate and most numerical methods will require step size reduction for stability.

To conclude, let us sketch out a possible way to attach the problem of volume-preserving integration for generic functions: let $\{\mathbf{\phi}_i(\mathbf{x})\}$ a set of basis functions. Consider the decomposition $$p(\mathbf{x}) = \nabla \cdot \mathbf{f}(\mathbf{x}) = \sum_i p_i \mathbf{\phi}_i(\mathbf{x})$$ in the given basis. Then, for each $i$, because of the independence of the basis functions, one must have $p_i=0$. Thus, as long as we are able to recover all the terms in the given vector field contributing to the coefficient $p_i$, we will automatically have a volume-preserving split vector field. The difficulty is to find suitable basis functions, in the sense that the corresponding vector fields must be easy to integrate exactly. One such example is the monomial basis for polynomials, $\mathbf{x}^\mathbf{j}$, for which we have demonstrated the existence of exact formulas.

Acknowledgements {#acknowledgements .unnumbered}
================

The work of the first author has been supported by the Norwegian Research Council through the project “Geometric Numerical Integration and Applications”, NFR project no. 191178/V30.

[10]{}

P. Chartier and A. Murua. . , 27(2):381, 2007.

Alex J. Dragt and Dan T. Abell. Symplectic maps and computation of orbits in particle accelerators. In [*Integration algorithms and classical mechanics ([T]{}oronto, [ON]{}, 1993)*]{}, volume 10 of [*Fields Inst. Commun.*]{}, pages 59–85. Amer. Math. Soc., Providence, RI, 1996.

A. Iserles, G. R. W. Quispel, and P. S. P. Tse. B-series methods cannot be volume-preserving. , 47(2):351–378, 2007.

Feng Kang and Zai Jiu Shang. Volume-preserving algorithms for source-free dynamical systems. , 71(4):451–463, 1995.

R. I. McLachlan, H. Z. Munthe-Kaas, G. R. W. Quispel, and A. Zanna. Explicit volume-preserving splitting methods for linear and quadratic divergence-free vector fields. , 8(3):335–355, 2008.

R.I. McLachlan and G.R.W. Quispel. Explicit geometric integration of polynomial vector fields. , 44:515–538, 2004.

Robert I. Mclachlan and Clint Scovel. A survey of open problems in symplectic integration. , 10:151–180, 1998.

Anatoly Neishtadt, Dmitri Vainchtein, and Alexei Vasiliev. Adiabatic invariance in volume-preserving systems. In [*I[UTAM]{} [S]{}ymposium on [H]{}amiltonian [D]{}ynamics, [V]{}ortex [S]{}tructures, [T]{}urbulence*]{}, volume 6 of [*IUTAM Bookser.*]{}, pages 89–107. Springer, Dordrecht, 2008.

M.-Z. Qin and W.-J. Zhu. Volume-[P]{}reserving [S]{}chemes and [N]{}umerical [E]{}xperiments. , 26:33–42, 1993.

G. R. W. Quispel. Volume-preserving integrators. , 206(1-2):26–30, 1995.

Govindan Rangarajan. Symplectic completion of symplectic jets. , 37(9):4514–4542, 1996.

J. C. Scovel. Symplectic numerical integration of hamiltonian systems. In Tudor Ratiu, editor, [*The Geometry of Hamiltonian Systems*]{}, volume 22, pages 463–496. MSRI, Springer-Verlag, New York, 1991. Symplectic numerical integration of Hamiltonian systems.

Zai Jiu Shang. Construction of volume-preserving difference schemes for source-free systems via generating functions. , 12(3):265–272, 1994.

H. A. Stone, Ali Nadim, and Steven H. Strogatz. Chaotic streamlines inside drops immersed in steady [S]{}tokes flows. , 232:629–646, 1991.

[^1]: Email: <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>, <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>