1. Introduction
Around 20 years ago, Stephen Anco and George Bluman [Reference Anco and Bluman2, Reference Anco and Bluman3] introduced a comprehensive practical method for determining conservation laws of partial differential equations (PDEs) in Kovalevskaya form. The method is based on finding adjoint symmetries and applying Helmholtz conditions.Footnote ^{1} A key part of the calculation is the inversion of the total divergence operator $\mathrm{Div}$ to obtain the components of the conservation law. Usually, this can be done by using a homotopy operator, but the following three problems may occur with the standard homotopy formula (which is given by Olver in [Reference Olver8]).

1. The homotopy formula uses definite integrals. This works well if the divergence is a differential polynomial; by contrast, rational polynomials commonly have a singularity at one limit. Hickman [Reference Hickman7] and Poole & Hereman [Reference Poole and Hereman9] suggest avoiding these by working in terms of indefinite integrals, an approach that we use throughout this paper. Alternatively, one can move the singularity by modifying the dependent variable (see Anco & Bluman [Reference Anco and Bluman3] and Poole & Hereman [Reference Poole and Hereman9]).

2. Scaling is fundamental to the homotopy approach to inversion. For instance, varying the scaling parameter in the standard homotopy formula moves contributions from variables along a ray to the origin. However, scaling does not change rational polynomials that are homogeneous of degree zero, so the standard inversion process does not work for such terms. Deconinck & Nivala [Reference Deconinck and Nivala5] discussed this problem in some detail (for one independent variable only) and suggested possible workarounds, but commented, ‘We are unaware of a homotopy method that algorithmically avoids all problems like the ones demonstrated $\ldots$ ’. Poole & Hereman [Reference Poole and Hereman9] proposed an approach that works well for problems with one independent variable, but noted the difficulties of extending this to multiple independent variables (in a way that can be programmed).

3. The standard homotopy operator applies to a starshaped domain and integrates along rays to the origin, changing all Cartesian coordinates at once. This gives an inefficient inversion, in the sense that the number of terms is generally very much greater than necessary; the homotopy formula creates superfluous elements of $\mathrm{ker}(\mathrm{Div})$ . For polynomial divergences, Poole & Hereman [Reference Poole and Hereman9] removed curls by parametrizing all terms given by the homotopy formula and optimizing the resulting linear system. This approach is very effective, because (except where there are cancellations) the homotopy formula tends to include every possible term that can appear in an inversion. However, it is unclear whether this approach can be generalized to nonpolynomial divergences. Moreover, the removal of curl terms takes extra processing time and does not allow for the possibility that new terms might appear in the most concise inversion of the divergence. If inversion could be done with respect to one independent variable at a time, this might prevent the occurrence of superfluous curls from the outset.
Example 1.1. To illustrate the inefficiency of homotopy operators on starshaped domains, consider the following divergence in $\mathbb{R}^3$ :
The homotopy operator based on a starshaped domain gives ${\mathcal{C}}=\mathrm{div}(x\phi,y\phi,z\phi )$ , where
By comparison, for a given divergence ${\mathcal{C}}(x,y,z)$ that has no singularities on the coordinate axes, using a homotopy formula that integrates one variable at a time gives ${\mathcal{C}}=\mathrm{div}(F,G,H)$ , where
This recovers the concise form $(F,G,H)=(x^2y\cos z,0,0)$ . However, the use of the lower limit makes the formula overelaborate (and unsuitable for treating singularities). Indefinite integration is far more straightforward:
So the homotopy formula for starshaped domains gives $14$ more terms than the simple form above; the superfluous terms amount to
The current paper extends the efficient onevariableatatime approach to total derivatives. Indefinite integration is used, as advocated by Hickman [Reference Hickman7] for standard homotopy operators; it avoids the complications resulting from singularities. From the computational viewpoint, the biggest advantage of integration with respect to just one independent variable is that the major computer algebra systems have efficient procedures for computing antiderivatives.
The keys to inverting a total divergence one variable at a time are ‘partial Euler operators’ (see Section 3). These enable the inversion of a total derivative $D_x$ to be written as an indefinite line integral. Section 4 introduces a new iterative method for inverting a given total divergence; typically, this does not produce superfluous terms and very few iterations are needed. Furthermore, it can cope with components that are unchanged by the relevant scaling operator.
The methods in this paper are systematic, leading to procedures that are intended to be suitable for implementation in computer algebra systems.
2. Standard differential and homotopy operators
Here is a brief summary of the standard operators that are relevant to total divergences; for further details, see Olver [Reference Olver8]. The independent variables $\textbf{x} = (x^1,\ldots, x^p)$ are local Cartesian coordinates, and the dependent variables $\textbf{u}=(u^1,\ldots,u^q)$ may be real or complexvalued. The Einstein summation convention is used to explain the main ideas and state general results. In examples, commonly used notation is used where this aids clarity. All functions are assumed to be locally smooth, to allow the key ideas to be presented simply.
Derivatives of each $u^\alpha$ are written as $u^\alpha _{\textbf{J}}$ , where ${\textbf{J}}=(j^1,\ldots,j^p)$ is a multiindex; each $j^i$ denotes the number of derivatives with respect to $x^i$ , so $u^\alpha _{\textbf{0}}= u^\alpha$ . The variables $x^i$ and $u^\alpha _{\textbf{J}}$ can be regarded as jet space coordinates. The total derivative with respect to $x^i$ ,
treats each $u^\alpha _{\textbf{J}}$ as a function of $\textbf{x}$ . To keep the notation concise, write
note that $u^\alpha _{\textbf{J}}=D_{\textbf{J}}(u^\alpha )$ . Let $[{\textbf{u}}]$ represent $\textbf{u}$ and finitely many of its derivatives; more generally, square brackets around an expression denote the expression and as many of its total derivatives as are needed.
A total divergence is an expression of the form
(If all $F^i$ depend on $\textbf{x}$ only, $\mathcal{C}$ is an ordinary divergence.) A conservation law of a given system of partial differential equations (PDEs), ${\mathcal{A}}_\ell ({\textbf{x}},[{\textbf{u}}])=0,\ \ell =1,\ldots,L$ , is a total divergence that is zero on all solutions of the system; each $F^i$ is a finite sum of terms. By using elementary algebraic operations (in particular, expanding logarithms of products and products of sums), the number of linearly independent terms may be maximized. When the number of linearly independent terms is maximal for each $i$ , we call the result the fully expanded form of $\textbf{F}$ .
When $p\gt 1$ , the $p$ tuple of components, $\textbf{F}$ , is determined by $\mathcal{C}$ up to a transformation of the form
(If $p=3$ , such a transformation adds a total curl to $\textbf{F}$ .) The total number of terms in $\textbf{F}$ is the sum of the number of terms in all of the fully expanded components $F^i$ . If this cannot be lowered by any transformation (2.1), we call $\textbf{F}$ minimal. Commonly, there is more than one minimal $\textbf{F}$ , any of which puts the inversion of $\mathrm{Div}$ in as concise a form as possible. If $p=1$ , the sole component $F^1$ (also denoted $F$ ) is determined up to an arbitrary constant, so the number of nonconstant terms is fixed.
The formal adjoint of a differential operator (with total derivatives), $\mathcal{D}$ , is the unique differential operator ${\mathcal{D}}^{\dagger}$ such that
is a total divergence for all functions $f({\textbf{x}},[{\textbf{u}}])$ and $g({\textbf{x}},[{\textbf{u}}])$ . In particular,
Thus, the (standard) Euler–Lagrange operator corresponding to variations in $u^\alpha$ is
Total divergences satisfy a useful identity: a function ${\mathcal{C}}({\textbf{x}},[{\textbf{u}}])$ is a total divergence if and only if
Given a Lagrangian function $L({\textbf{x}},[{\textbf{u}}])$ , the Euler–Lagrange equations are ${\textbf{E}}_{u^\alpha }(L)=0$ . Given a set of Euler–Lagrange equations that are polynomial in the variables $({\textbf{x}},[{\textbf{u}}])$ , the function $\overline{L}$ given by the homotopy formula
differs from $L$ by a total divergence. (The same applies to many, but not all, nonpolynomial Euler–Lagrange equations.)
When $p=1$ , the equation $P(x,[{\textbf{u}}])=D_xF$ is invertible (at least, for polynomial $P$ ) by the following standard homotopy formula:
The operator acting on $P$ in the braces above is the higher Euler operator of order $i$ for $p=1$ . When $p\geq 2$ , the standard homotopy formula is similar, but somewhat more complex (see Olver [Reference Olver8] for details); it is based on higher Euler operators and integration along rays in a totally starshaped domain. The following example illustrates that even for quite simple divergences, this formula commonly yields inversions with many superfluous terms.
Example 2.1. The Benjamin–Bona–Mahony (BBM) equation, $u_{t}uu_xu_{xxt}=0$ , has a conservation law
The standard homotopy formula gives
a total of 17 terms. By contrast, careful integration by inspection yields
which is minimal, having only five terms in the components.
The homotopy formulae above can be applied or adapted to some, but not all, classes of nonpolynomial Lagrangians and divergences.
3. Partial Euler operators and partial scalings
This section introduces some ideas and results that underpin integration with respect to one independent variable at a time. The independent variable over which one integrates is distinguished; this is denoted by $x$ . For instance, if $x=x^1$ , replace the derivative index $\textbf{J}$ by $({\textbf{I}},j)$ , where $j=j^1$ and ${\textbf{I}}=(j^2,\ldots,j^{p})$ . So the dependent variables and their derivatives are denoted
In examples, however, we write each $u^\alpha _{\textbf{I}}$ more simply (as $u$ , $v_y$ , $u_{yt}$ , and so on), using, j for $D_x^{\,j}$ .
3.1 Partial Euler operators
The partial Euler operator with respect to $x$ and $u^\alpha _{\textbf{I}}$ is obtained by varying each $u^\alpha _{\textbf{I}}$ independently, treating $x$ as the sole independent variable:
Consequently, the standard Euler operator with respect to $u^\alpha$ amounts to
Similarly, the partial Euler operator with respect to $x$ and $u^\alpha _{\textbf{I},k}$ is
Note that
The following identities are easily verified; here, $f({\textbf{x}},[{\textbf{u}}])$ is an arbitrary function.
where the last term in (3.7) is zero if $j^i=0$ .
3.2 Inversion of $\boldsymbol{D}_\boldsymbol{x}$
The identities (3.5), (3.6) and (3.7) are central to the inversion of total divergences, including the following inversion of $P=D_xF$ as an indefinite line integral.
Lemma 3.1. If $P(\boldsymbol{x},[\boldsymbol{u}])=D_xF$ , then, up to an (irrelevant) arbitrary function of all independent variables other than $x$ ,
Proof. By the identity (3.6),
Moreover,
Substituting $P$ for $D_xF$ completes the proof.
Example 3.1. Locally, away from its singularities, the function
belongs to $\mathrm{im}(D_x)$ , but cannot be inverted using the standard homotopy formula. Substituting
into (3.8) yields the inversion:
3.3 Integration by parts
From here on, we will restrict attention to total divergences $\mathcal{C}$ whose fully expanded form has no terms that depend on $\textbf{x}$ only. Such terms can be inverted easily by evaluating an indefinite integral, as explained in the Introduction. Henceforth, all indefinite integrals denote antiderivatives with the minimal number of terms in their fully expanded form. Any arbitrary constants and functions that would increase the number of terms are set to zero. This restriction facilitates the search for minimal inversions.
The indefinite line integral formula (3.8) is closely related to integration by parts. To see this, we introduce a positive ranking on the variables $u^\alpha _{{\textbf{J}}}\,$ ; this is a total order $\preceq$ that is subject to two conditions:
The leading part of a differential function is the sum of terms in the function that depend on the highestranked $u^\alpha _{\textbf{J}}$ , and the rank of the function is the rank of its leading part (see Rust et al. [Reference Rust, Reid and Wittkopf12] for details and references). Let $f({\textbf{x}},[{\textbf{u}}])$ denote the leading part of the fully expanded form of $F$ and let $\mathrm{U}_{,k}$ denote the highestranked $u^\alpha _{\textbf{I},k}$ ; then the highestranked part of $P=D_xF$ is $\mathrm{U}_{,k+1}\partial f/\partial \mathrm{U}_{,k}$ . Then (3.8) includes the contribution
Integration by parts gives the same result. Subtracting $f$ from $F$ and iterating shows that evaluating the line integral (3.8) is equivalent to integrating by parts from the highestranked terms downwards.
Integration by parts is useful for splitting a differential expression $P({\textbf{x}},[{\textbf{u}}])$ , with $P({\textbf{x}},[0])=0$ , into $D_xF$ and a remainder, $R$ , whose $x$ derivatives are of the lowestpossible order. The splitting is achieved by the following procedure.
Procedure A. Integration by parts

Step 0. Choose a positive ranking in which $u^\alpha _{{\textbf{I}},0}\prec D_xu^\beta$ for all $\alpha,{\textbf{I}}$ and $\beta$ . (We call such rankings $x$ dominant.) Initialize by setting $F\;:\!=\;0$ and $R\;:\!=\;0$ .

Step 1. Identify the highestranked $u^\alpha _{\textbf{I},k}$ in $P$ ; denote this by $\mathrm{U}_{,k}$ . If $k=0$ , add $P$ to $R$ and stop. Otherwise, determine the leading part, $g$ , of $P$ .

Step 2. Determine the sum $h\mathrm{U}_{,k}$ of all terms in the fully expanded form of $g$ that are of the form $\gamma \mathrm{U}_{,k}$ , where $\gamma$ is ranked no higher than $\mathrm{U}_{,k1}$ , and let
\begin{equation*} H=\int h\,\mathrm {d}\mathrm {U}_{,k1}\,. \end{equation*} 
Step 3. Update $F, R$ and $P$ , as follows:
\begin{equation*} F\;:\!=\;F+H,\qquad R\;:\!=\;R+gh\mathrm {U}_{,k},\qquad P\;:\!=\;Pg+h\mathrm {U}_{,k}\!\!D_xH. \end{equation*}If $P\neq 0$ , return to Step 1. Otherwise output $F$ and $R$ , then stop.
The reason for choosing an $x$ dominant ranking is to ensure that the derivative order with respect to $x$ outweighs all other ranking criteria. Consequently, the minimally ranked remainder cannot contain $x$ derivatives of unnecessarily high order.
Example 3.2. To produce a concise inversion of a conservation law of the Harry Dym equation (see Example 4.1 below), it is necessary to split
Procedure A gives the splitting
Example 3.3. The ranking criterion in Step 2 of Procedure A ensures that there are no infinite loops. It is not enough that terms are linear in the highestranked $x$ derivative, as shown by the following splitting of
For the positive $x$ dominant ranking defined by $v\prec u\prec v_y$ , Procedure A yields
Both terms in the remainder are linear in their highestranked $x$ derivatives, which are $v_{,2}$ and $v_{y,1}$ , respectively. However, further integration by parts would return $P$ to a form with a higherranked remainder.
3.4 Partial scalings
To investigate partial Euler operators further, it is helpful to use a variant of the homotopy approach. The partial scaling (by a positive real parameter, $\lambda$ ) of a function $f({\textbf{x}},[{\textbf{u}}])$ with respect to $x$ and $u^\alpha _{\textbf{I}}$ is the mapping
Again, each $u^\alpha _{\textbf{I}}$ is treated as a distinct dependent variable. Note the identity
Definition 3.1. The partial scaling $\sigma _{u^\alpha _\boldsymbol{I}}^{x}$ is a good scaling for a given differential function $f(\boldsymbol{x},[\boldsymbol{u}])$ if
for all $\lambda$ in some neighborhood of $1$ .
By definition, the partial scaling $\sigma _{u^\alpha _{\textbf{I}}}^{x}$ fails to be a good scaling for $f$ if and only if there are terms that are independent of $\lambda$ in the fully expanded form of $\sigma _{u^\alpha _{\textbf{I}}}^{x}(f;\;\lambda )$ . The simplest cause of this is that the fully expanded form of $f$ has terms that are independent of $u^\alpha _{\textbf{I}}$ and its $x$ derivatives. However, this is not the only cause, as the following example illustrates.
Example 3.4. The scalings $\sigma ^y_{u}$ and $\sigma ^y_v$ are not good scalings for
because (in fully expanded form),
The terms in braces are independent of $\lambda$ ; in (3.14) (resp. (3.15)), some of these depend on $u$ (resp. $v$ ) and/or its $y$ derivatives. Part of the scaled logarithmic term is independent of $\lambda$ , though part survives differentiation. Note that $\sigma ^y_{u}$ is a good scaling for the term $u_{x}/u^2$ ; the singularity at $u=0$ is not an obstacle.
Lemma 3.2. The partial scaling $\sigma _{u^\alpha _\boldsymbol{I}}^{x}$ is a good scaling for $f(\boldsymbol{x},[\boldsymbol{u}])$ if and only if
Proof. If $\sigma _{u^\alpha _\boldsymbol{I}}^{x}$ is a good scaling, (3.16) is a consequence of $f=\sigma _{u^\alpha _{\textbf{I}}}^{x}(f;\;1)$ and local smoothness. Conversely, suppose that (3.16) holds and let $\mu$ be a positive real parameter that is independent of $\lambda$ . Then for $\mu$ sufficiently close to $1$ ,
Therefore, $\sigma _{u^\alpha _{\textbf{I}}}^{x}$ is a good scaling.
The use of the limit in (3.16) is needed to deal with any values of $u^\alpha _{\textbf{I},k}$ for which the integral is an indeterminate form. For other values, simple substitution of $\lambda =1$ gives the limit.
Given a partial scaling $\sigma _{u^\alpha _{\textbf{I}}}^{x}$ and a differential function $f({\textbf{x}},[{\textbf{u}}])$ , let
For $\mu$ sufficiently close to $1$ (using $\widehat{f}$ as shorthand for $\pi _{u^\alpha _{\textbf{I}}}^{x}(f)$ ),
the first equality comes from the proof of Lemma 3.2. Therefore, $\sigma _{u^\alpha _{\textbf{I}}}^{x}$ is a good scaling for $\pi _{u^\alpha _{\textbf{I}}}^{x}(f)$ . Moreover, there are no terms in the fully expanded form of the remainder, $f\pi _{u^\alpha _{\textbf{I}}}^{x}(f)$ , for which $\sigma _{u^\alpha _{\textbf{I}}}^{x}$ is a good scaling, because
So $\pi _{u^\alpha _{\textbf{I}}}^{x}$ is the projection that maps a given function onto the component which has $\sigma _{u^\alpha _{\textbf{I}}}^{x}$ as a good scaling.
Definition 3.2. The partial scaling $\sigma _{u^\alpha _\boldsymbol{I}}^{x}$ is a poor scaling for a given differential function $f(\boldsymbol{x},[\boldsymbol{u}])$ if $f\pi _{u^\alpha _\boldsymbol{I}}^{x}(f)$ depends on any $u^\alpha _{\boldsymbol{I},k}\,$ .
For instance, both $\sigma ^y_{u}$ and $\sigma ^y_v$ are poor scalings for (3.13), as explained in Example 3.4. Section 4.3 addresses the problem of inverting divergences such as (3.13) that have poor scalings. First, we develop the inversion process for general divergences. The following results are fundamental.
Theorem 3.1. Let $f(\boldsymbol{x},[\boldsymbol{u}])$ be a differential function.

1. If $f=D_x F$ , then $f\in \mathrm{ker}(\boldsymbol{E}_{u^\alpha _\boldsymbol{I}}^{x})$ for all $\alpha$ and $\boldsymbol{I}$ ; moreover,
(3.18) \begin{equation} \pi _{u^\alpha _\boldsymbol{I}}^{x}(F)=\lim _{\lambda \rightarrow 1}\int \sum _{j\geq 0} u^\alpha _{\boldsymbol{I},j}\,\sigma _{u^\alpha _\boldsymbol{I}}^{x}\!\left (\boldsymbol{E}^x_{u^\alpha _{\boldsymbol{I},j+1}}(f);\;\lambda \right )\,\mathrm{d} \lambda. \end{equation} 
2. If $f\in \mathrm{ker}(\boldsymbol{E}_{u^\alpha _\boldsymbol{I}}^{x})$ , then $\pi _{u^\alpha _\boldsymbol{I}}^{x}(f)\in \mathrm{im}(D_x)$ .

3. If $g=\boldsymbol{E}_{u^\alpha _\boldsymbol{I}}^{x}(f)$ , then, up to terms in $\mathrm{im}(D_x)$ ,
(3.19) \begin{equation} \pi _{u^\alpha _\boldsymbol{I}}^{x}(f)=\lim _{\lambda \rightarrow 1}\int u^\alpha _\boldsymbol{I}\,\sigma _{u^\alpha _\boldsymbol{I}}^{x}\!\left (g;\;\lambda \right )\,\mathrm{d} \lambda. \end{equation}
Proof. All three statements are proved by expanding $\pi _{u^\alpha _{\textbf{I}}}^{x}(f)$ :
Here, $h({\textbf{x}},[{\textbf{u}}])$ is obtained by integrating by parts, using the identity (3.11).
If $f=D_x F$ , the identity (3.5) amounts to $f\in \mathrm{ker}({\textbf{E}}_{u^\alpha _{\textbf{I}}}^{x})$ . Replace $f$ by $F$ in (3.20) and use the identity (3.6) to obtain (3.18). Statements $2$ and $3$ come directly from (3.21).
Note that (3.18) is a homotopy formula for (at least partially) inverting $D_xF$ , giving a third way to do this. The line integral formula (3.8) carries out the full inversion in one step, but may take longer to compute.
4. The inversion method for Div
This section introduces a procedure to invert $\mathrm{Div}$ , with a ranking heuristic (informed by experience) that is intended to keep the calculation short and efficient. To motivate the procedure, it is helpful to examine a simple example.
Example 4.1. Wolf et al. [Reference Wolf, Brand and Mohammadzadeh14] introduced computer algebra algorithms that can handle general (nonpolynomial) conservation laws and used these to derive various rational conservation laws of the Harry Dym equation. In the (unique) $x$ dominant positive ranking, the equation is ${\mathcal{A}}=0$ , with
The highestorder conservation law derived in Wolf et al. [Reference Wolf, Brand and Mohammadzadeh14] is ${\mathcal{C}}={\mathcal{Q}}{\mathcal{A}}$ , where
Note that $\sigma ^x_{u}$ is a good scaling for $\mathcal{C}$ . The first step in inverting ${\mathcal{C}}=D_xF+D_tG$ is to apply the partial Euler operator $\boldsymbol{E}_{u}^{x}$ , to annihilate the term $D_xF$ . There are only two independent variables, so the identity (3.2) shortens the calculation to
Applying (3.19), then using Procedure A to integrate by parts (see Example 3.2) gives
Therefore,
where $\boldsymbol{E}_{u}^{x}(\widetilde{{\mathcal{C}}})=0$ . As $\sigma ^x_{u}$ is a good scaling for
the second part of Theorem 3.1 states that $\widetilde{{\mathcal{C}}}\in \mathrm{im}(D_x)$ ; consequently,
Either the line integral (3.8) or Procedure A completes the inversion, giving $\widetilde{{\mathcal{C}}}=D_xF$ , where
The fully expanded form of $(F,G)$ is minimal, having $11$ terms rather than the $12$ terms in Wolf et al. [Reference Wolf, Brand and Mohammadzadeh14]. (Note: there is an equivalent conservation law, not in the form ${\mathcal{Q}}{\mathcal{A}}$ , that has only $10$ terms.)
4.1 A single iteration
The basic method for inverting a given total divergence one independent variable at a time works similarly to the example above. Suppose that after $n$ iterations the inversion process has yielded components $F^i_n$ and that an expression of the form ${\mathcal{C}}=D_if_n^i$ remains to be inverted. For the next iteration, let ${\textbf{E}}_{u}^{x}$ be the partial Euler operator that is applied to $\mathcal{C}$ . Here, $u$ is one of the variables $u^\alpha _{\textbf{I}}$ , which is chosen to ensure that for each $i$ such that $x^i\neq x$ ,
This requires care, in view of the identity (3.7). However, it is achievable by using the variables $u^\alpha _{\textbf{I}}$ in the order given by a ranking that is discussed in Section 4.2. This ranking is entirely determined by userdefined rankings of the variables $x^j$ and $u^\alpha$ .
Taking (3.5) into account leads to the identity
which, with together with Theorem 3.1, is the basis of the inversion method. The method works without modification provided that:

• there are no poor scalings for any terms in $\mathcal{C}$ ;

• the fully expanded form of $\mathcal{C}$ has no terms that are linear in $[{\textbf{u}}]$ .
We begin by restricting attention to divergences for which these conditions hold, so that
where every term in ${\textbf{E}}_{u}^{x}(\pi _u^x f_n^i)$ depends on $[{\textbf{u}}]$ . The modifications needed if either condition does not hold are given in Sections 4.3 and 4.4.
The iteration of the inversion process runs as follows. Calculate ${\textbf{E}}_{u}^{x}({\mathcal{C}})$ , which is a divergence $D_iP^i$ with no $D_x$ term, by (4.3); it involves at most $p1$ (but commonly, very few) nonzero functions $P^i$ . Invert this divergence, treating $x$ as a parameter. If it is possible to invert in more than one way, always invert into the $P^i$ for which $x^i$ is ranked as low as possible; the reason for this is given in the next paragraph. If ${\textbf{E}}_{u}^{x}({\mathcal{C}})$ has nonlinear terms that involve derivatives with respect to more than one $D_i$ (excluding $D_x$ ), this is accomplished by iterating the inversion process with as few independent variables as are needed. Otherwise, $P^i$ can be determined more quickly by using integration by parts (Procedure A, with $x$ replaced by the appropriate $x^i$ ), and/or the method for linear terms (see Procedure B in Section 4.4). Note that this shortcut can be used whenever there are only two independent variables.
At this stage, check that the fully expanded form of each $P^i$ has no terms that are ranked lower than $u$ . If any term is ranked lower than $u$ , stop the calculation and try a different ranking of the variables $x^j$ and/or $u^\alpha$ . This is essential because, to satisfy (4.1) and avoid infinite loops, the variables $u^\alpha _{\textbf{I}}$ that are chosen to be $u$ in successive iterations must progress upwards through the ranking. Where there is a choice of inversion, the rank of each term in $P^i$ is maximized by using the $x^i$ of minimum order; this avoids unnecessary reranking.
Having found and checked $P^i$ , use (3.21) to obtain
for arbitrary functions $h^i({\textbf{x}},[{\textbf{u}}])$ . Apply Procedure A to the function in braces and choose $h^i$ to make the righthand side of (4.4) equal the remainder from this procedure. This yields the representation of $\pi _u^x(f_n^i)$ that has the lowestorder derivatives (with respect to $x$ ) consistent with the inversion of $D_iP^i$ ; call this representation $f^i$ . Commonly, such a lowestorder representation is needed to obtain a minimal inversion.
By Theorem 3.1, there exists $\phi$ such that
because (by construction) the expression in parentheses belongs to $\mathrm{ker}({\textbf{E}}_{u}^{x})$ . Use the line integral formula (3.8) or Procedure A to obtain $\phi$ , then set $f^i\;:\!=\;\phi$ for $x^i=x$ . Now update: set
4.2 Ranking and using the variables
Having described a single iteration, we now turn to the question of how to choose $x$ and $u$ effectively. The starting point is to construct a derivativedominant ranking of the variables $u^\alpha _{\textbf{J}}$ . This is a positive ranking that is determined by:

• a ranking of the independent variables, $x^1\prec x^2\prec \cdots \prec x^p$ ;

• a ranking of the dependent variables, $u^1\prec u^2\prec \cdots \prec u^q$ .
(Later in this section, we give a heuristic for ranking the dependent and independent variables effectively.) The derivativedominant ranking (denoted ${\textbf{u}}_p$ ) is constructed iteratively, as follows.
In practice, very few $u^\alpha _{\textbf{J}}$ are needed to carry out many inversions of interest, but it is essential that these are used in the order given by their ranking, subject to a constraint on ${\textbf{I}}$ that is explained below.
Given an independent variable, $x$ , we call $u^\alpha _{\textbf{I}}$ relevant if the updated $\mathcal{C}$ depends on $u^\alpha _{\textbf{I},k}$ for some $k\geq 0$ . The first set of iterations uses $x=x^1$ . For the initial iteration, $u$ is the lowestranked relevant $u^\alpha$ . In the following iteration, $u$ is the nextlowestranked relevant $u^\alpha$ and so on, up to and including $u^q$ . (From (3.7), the condition (4.1) holds whenever $u=u^\alpha,\ \alpha = 1,\ldots, q$ .) After these iterations, the updated $\mathcal{C}$ is independent of $\textbf{u}$ and its unmixed $x$ derivatives.
If the updated $\mathcal{C}$ has any remaining $x$ derivatives, these are mixed. Thus, as $\mathcal{C}$ has no linear terms, a necessary condition for the inversion to be minimal is that every $f_n^i$ is independent of $\textbf{u}$ and its unmixed $x$ derivatives. Consequently, (4.1) holds for $u=u^\alpha _{\textbf{I}}$ whenever ${\textbf{I}}=1$ , because
Therefore, the process can be continued using each relevant $u=u^\alpha _{\textbf{I}}$ with ${\textbf{I}}=1$ in the ranked order. Iterating, the same argument is used with ${\textbf{I}}=2,3,\ldots$ , until $\mathcal{C}$ is independent of $x$ derivatives. Now set $x=x^2$ and iterate, treating $x^1$ as a parameter. In principle, this can be continued up to $x=x^p$ ; in practice, only a very few iterations are usually needed to complete the inversion. The best rankings invert many terms at each iteration. On the basis of some experience with conservation laws, the following heuristic for ranking the variables $x^j$ and $u^\alpha$ is recommended.
Ranking heuristic. Apply criteria for ranking independent variables, using the following order of precedence.

1. Any independent variables that occur in the arguments of arbitrary functions of $\textbf{x}$ should be ranked as high as possible, if they multiply terms that are nonlinear in $[{\textbf{u}}]$ . For instance, if nonlinear terms in a divergence depend on an arbitrary function, $f(t)$ , set $x^p=t$ .

2. If independent variables occur explicitly in nonarbitrary functions, they should be ranked as high as possible (subject to 1 above), with priority going to variables with the most complicated functional dependence. For instance, if $\mathcal{C}$ is linear in $x^i$ and quadratic in $x^j$ , then $x^i\prec x^j$ (so $i\lt j$ in our ordering).

3. If an unmixed derivative of any $u^\alpha$ with respect to $x^i$ is the argument of a function other than a rational polynomial, rank $x^i$ as low as possible.

4. Set $x^i\prec x^j$ if the highestorder unmixed derivative (of any $u^\alpha$ ) with respect to $x^i$ is of higher order than the highestorder unmixed derivative with respect to $x^j$ .

5. Set $x^i\prec x^j$ if there are more occurrences of unmixed $x^i$ derivatives (in the fully expanded divergence) than there are of unmixed $x^j$ derivatives.

6. Apply criteria $3$ , $4$ , and $5$ in order of precedence, replacing unmixed by ‘minimally mixed’ derivatives. Minimally mixed means that there are as few derivatives as possible with respect to any other variable(s).
The derivative indices $\textbf{J}$ in a derivativedominant ranking are ordered according to the ranking of the independent variables. This can be used to rank the dependent variables; if there is more than one dependent variable in $\mathcal{C}$ , use the following criteria in order.

1. Let $u^\alpha \prec u^\beta$ if $\mathcal{C}$ is linear in $[u^\alpha ]$ and nonlinear in $[u^\beta ]$ .

2. Let $u^\alpha \prec u^\beta$ if the lowestranked derivative of $u^\alpha$ that occurs in $\mathcal{C}$ is ranked lower than the lowestranked derivative of $u^\beta$ in $\mathcal{C}$ . (In conservation laws, the lowestranked derivative of $u^\alpha$ is commonly the undifferentiated $u^\alpha$ , which corresponds to ${\textbf{J}}=0$ .)

3. Let $u^\alpha \prec u^\beta$ if the lowestranked derivative of $u^\alpha$ in the fully expanded form of $\mathcal{C}$ occurs in more terms than the corresponding derivative of $u^\beta$ does.
These two sets of criteria are not exhaustive (allowing ties, which must be broken), but the aim that underlies them is to carry as few terms as possible into successive iterations of the procedure. Partial Euler operators with respect to unmixed derivatives are used in the earliest iterations; commonly, these are sufficient to complete the inversion.
4.3 How to deal with poor scalings
To remove an earlier restriction on the inversion process, we now address the problem of poor scalings, namely, that $u^\alpha _{\textbf{I}}$ and its $x$ derivatives (denoted $[u^\alpha _{\textbf{I}}]_x$ ) may occur in terms that belong to ${\mathcal{C}}\pi _{u^\alpha _{\textbf{I}}}^{x}({\mathcal{C}})$ . Such terms (when fully expanded) are products of homogeneous rational polynomials in $[u^\alpha _{\textbf{I}}]_x$ of degree zero and logarithms of a single element of $[u^\alpha _{\textbf{I}}]_x$ . We refer to these collectively as zerodegree terms.
To overcome this difficulty, we modify the approach used by Anco & Bluman [Reference Anco and Bluman3] to treat singularities. In our context, $[u^\alpha _{\textbf{I}}]_x$ is replaced by $[u^\alpha _{\textbf{I}}+U^\alpha _{\textbf{I}}]_x$ , where $U^\alpha _{\textbf{I}}$ is regarded as a new dependent variable that is ranked higher than $u^\alpha _{\textbf{I}}$ . This approach works equally well for logarithms, ensuring that $\pi _{u^\alpha _{\textbf{I}}}^{x}$ is a good scaling for all terms that depend on $[u^\alpha _{\textbf{I}}]_x$ , so that its kernel consists only of terms that are independent of these variables. At the end of the calculation, all members of $[U^\alpha _{\textbf{I}}]_x$ are set to zero. Note that there is no need to replace $u^\alpha _{\textbf{I}}$ in terms that are not zerodegree in $[u^\alpha _{\textbf{I}}]_x$ , as total differentiation preserves the degree of homogeneity.Footnote ^{2}
Example 4.2. To illustrate the inversion process for divergences that have zerodegree terms, we complete Example 3.4 by inverting
this has no linear terms. The ranking heuristic gives $y\prec x$ and $u\prec v$ . The rest of the calculation goes as follows.

1. ${\mathcal{C}}\pi _u^y({\mathcal{C}})$ has just one zerodegree term (in $[u]_y$ ), namely $(2u_y/u)\ln u$ . Replace this term by $(2(u_y+U_y)/(u+U))\ln u+U$ .

2. Calculate $\boldsymbol{E}_u^y({\mathcal{C}})=2u_x+v_{xy}2u_xu^{3}=D_x\{2u+v_y+u^{2}\}=D_x\{\boldsymbol{E}_u^y(u^2+uv_yu^{1})\}$ . No term in $2u+v_y+u^{2}$ is ranked lower than $u$ , as $u$ is the lowestranked variable.

3. Then $\pi _u^y({\mathcal{C}}\!\!D_x\{u^2+uv_yu^{1}\})=D_y\{uv_x+(\ln u+U)^2\}$ .

4. Now set $U=0$ to yield the remainder ${\mathcal{C}}_1={\mathcal{C}}\!\!D_x\{u^2+uv_yu^{1}\}\!\!D_y\{uv_x+(\ln u)^2\}$ at the close of the first iteration. This amounts to ${\mathcal{C}}_1=2v_xv_{yy}+v_{yy}/v_y$ .

5. The second iteration starts with ${\mathcal{C}}_1\pi _v^y({\mathcal{C}}_1)=v_{yy}/v_y$ ; as this term is zerodegree in $[v]_y$ , replace it by $(v_{yy}+V_{yy})/(v_y+V_y)$ .

6. Calculate $\boldsymbol{E}_v^y({\mathcal{C}}_1)=2v_{xyy}=D_x\{2v_{yy}\}=D_x\{E_v^y(vv_{yy})\}=D_x\{E_v^y(v_y^2)\}$ . Note that $2v_{yy}$ is not ranked lower than $v$ .

7. Then $\pi _v^y({\mathcal{C}}_1\!\!D_x\{v_y^2\})=D_y\{2v_xv_y+\ln v_y+V_y\}$ .

8. Now set $V=0$ to yield the remainder ${\mathcal{C}}_2={\mathcal{C}}_1\!\!D_x\{v_y^2\}\!\!D_y\{2v_xv_y+\ln v_y\}=0$ at the close of the second iteration. The inversion process stops, having yielded the output
\begin{equation*} F=u^2+uv_yu^{1}+v_y^2,\qquad G=uv_x+(\ln u)^22v_xv_y+\ln v_y. \end{equation*}
Note that the homotopy formula (3.18) for inverting $D_xF$ can be adjusted in the same way, whenever $\sigma _{u^\alpha _\boldsymbol{I}}^{x}$ is a poor scaling for $\boldsymbol{E}^x_{u^\alpha _{\boldsymbol{I},j+1}}(D_xF)$ .
4.4 Linear divergences
The inversion process runs into a difficulty when a given divergence has terms that are linear in $[{\textbf{u}}]$ , with mixed derivatives. Then it is possible to invert in more than one way, some of which may not produce a minimal result. To address this, it is helpful to invert using a different process. Suppose that $\mathcal{C}$ is a linear divergence (in fully expanded form). Instead of using the derivativedominant ranking, integrate by parts, working down the total order ${\textbf{J}}$ of the derivatives $u^\alpha _{\textbf{J}}$ . For a given total order, we will invert the mixed derivatives first, though this is not essential.
Starting with the highestorder derivatives, one could integrate each term $f({\textbf{x}})u^\alpha _{\textbf{J}}$ in $\mathcal{C}$ by parts with respect to any $x^i$ such that $j^i\geq 1$ , yielding the remainder $\!\!D_i(f)u^\alpha _{{\textbf{J}}{\boldsymbol{1}}_i}$ . If $D_{\textbf{J}}$ is a mixed derivative, we seek to choose $D_i$ in a way that keeps the result concise. Here are some simple criteria that are commonly effective, listed in order of precedence.

1. $f$ is independent of $x^i$ .

2. $\mathcal{C}$ includes the term $D_i(f)u^\alpha _{{\textbf{J}}{\boldsymbol{1}}_i}$ .

3. $f$ is linear in $x^i$ .
These criteria can be used as an initial pass to invert $\mathcal{C}$ at least partially, leaving a remainder to be inverted that may have far fewer terms than $\mathcal{C}$ does.
Integrating the remainder by parts is straightforward if each $u^\alpha _{\textbf{J}}$ is an unmixed derivative. If $\textbf{J}$ denotes a mixed derivative, integrate with respect to each $x^i$ such that $j^i\geq 1$ in turn, multiplying each result by a parameter (with the parameters summing to $1$ ). Iterate until either the remainder has a factor that is zero for some choice of parameters or there is no remainder. The final stage is to choose the parameters so as to minimize the number of terms in the final expression. Although this produces a minimal inversion, it comes at the cost of extra computational time spent doing all possible inversions followed by the parameter optimization.
Example 4.3. Consider the linear divergence
where $f$ is an arbitrary function. The first of the simple criteria above yields
The second criterion is not helpful at this stage, but the third criterion gives
The remainder has no mixed derivatives; integrating it by parts produces the minimal inversion
Example 4.4. To illustrate the parametric approach, consider
The simple criteria are irrelevant at present, so instead introduce parameters $\lambda _l$ and consider all possible inversions of the mixed derivative terms. Stepbystep, one obtains the following.
As the remainder has a factor $(1\lambda _2)(1\lambda _1)$ , the inversion is complete if either parameter is set to $1$ . The minimal (twoterm) inversion has $\lambda _1=1$ , which gives ${\mathcal{C}}=D_xF+D_tG$ , where
For $\lambda _1\neq 1, \lambda _2=1$ , the inversion has three terms if $\lambda _1=0$ , or five terms otherwise.
Note the importance of factorizing the remainder to stop the calculation once a possible parameter choice occurs. If we had continued the calculation without setting either $\lambda _i$ to $1$ , it would have stopped at order zero, not one, giving an eleventerm inversion for general $\lambda _i$ .
To summarize, one can invert a divergence that is linear in $[{\textbf{u}}]$ by applying the following procedure.
Procedure B. Iinversion of a linear total divergence

Step 0. Identify the maximum derivative order, $N={\textbf{J}}$ , of the variables $u^\alpha _{\textbf{J}}$ that occur in $\mathcal{C}$ . Set $F^i\;:\!=\;0,\ i=1\ldots,p$ .

Step 1. While $\mathcal{C}$ has at least one term of order $N$ , do the following. Select any such term, $f({\textbf{x}})u^\alpha _{\textbf{J}}$ , and determine a variable $x^i$ over which to integrate. If desired, parametrize for mixed derivatives, as described above. Set $F^i\;:\!=\;F^i+fu^\alpha _{{\textbf{J}}{\boldsymbol{1}}_i}$ for the chosen $i$ (with the appropriate modification if parameters are used), and update the remainder by setting ${\mathcal{C}}\;:\!=\;{\mathcal{C}}f({\textbf{x}})u^\alpha _{\textbf{J}}\!\!D_i(f)u^\alpha _{{\textbf{J}}{\boldsymbol{1}}_i}$ . Once $\mathcal{C}$ has no terms of order $N$ , continue to Step 2.

Step 2. If $\mathcal{C}$ is nonzero and cannot be set to zero by any choice of parameters, set $N\;:\!=\;N1$ and return to Step 1; otherwise, set $\mathcal{C}$ to zero and carry out parameter optimization (if needed), give the output $F^i,\ i=1,\ldots,p$ , then stop.
4.5 Summary: a procedure for inverting Div
Having addressed potential modifications, we are now in a position to summarize the inversion process for any total divergence $\mathcal{C}$ whose fully expanded form has no terms depending on $\textbf{x}$ only.
Procedure C. Inversion of ${\mathcal{C}}=D_iF^i$

Step 0. Let ${\mathcal{C}}_\ell$ be the linear part of $\mathcal{C}$ . If ${\mathcal{C}}_\ell =0$ , set ${\mathcal{C}}_0\;:\!=\;{\mathcal{C}}$ and $F_0^i\;:\!=\;0,\ i=1,\ldots,p$ . Otherwise, invert ${\mathcal{C}}_\ell$ using the technique described in Procedure B above, to obtain functions $F_0^i$ that satisfy ${\mathcal{C}}_\ell =D_iF_0^i$ . Set ${\mathcal{C}}_0\;:\!=\;{\mathcal{C}}{\mathcal{C}}_\ell$ . Choose a derivativedominant ranking for ${\mathcal{C}}_0$ (either using the ranking heuristic or otherwise). In the notation used earlier, set $x\;:\!=\;x^1$ and $u$ to be the lowestranked relevant $u^\alpha$ ; typically, $u\;:\!=\;u^1$ . Set $n\;:\!=\;0$ ; here $n+1$ is the iteration number.

Step 1. Calculate ${\mathcal{C}}_n\pi _u^x({\mathcal{C}}_n)$ ; if this includes terms depending on $[u]_x$ , replace $[u]_x$ in these terms only by $[u+U]_x$ .

Step 2. Apply the process detailed in Section 4.1 (with ${\mathcal{C}}_n$ replacing $\mathcal{C}$ ). Provided that the ranking check is passed, this yields components $f^i$ , which may depend on $[U]_x$ . If the ranking check is failed, choose a different ranking of $x^j$ and $u^\alpha$ for the remainder of the inversion process and return to Step 1, starting with the lowestranked $x$ and (relevant) $u$ and working upwards at each iteration.

Step 3. Replace all elements of $[U]_x$ in $f^i$ by zero.

Step 4. Update: set $F^i_{n+1}\;:\!=\;F^i_n+f^i,\ {\mathcal{C}}_{n+1}\;:\!=\;{\mathcal{C}}_n\!\!D_if^i$ and $n\;:\!=\;n+1$ . If ${\mathcal{C}}_{n+1}= 0$ , output $F^i=F^i_{n+1}$ and stop. Otherwise, update $u$ and $x$ as detailed in Section 4.2 and return to Step 1.
Example 4.5. As an example with linear terms, the Khokhlov–Zabolotskaya equation,
has a conservation law (see Poole & Hereman [11]) that involves an arbitrary function, $f(t)$ :
The linear part ${\mathcal{C}}_\ell$ of $\mathcal{C}$ is inverted by Procedure B, as shown in Example 4.3. From (4.6), Step 0 in Procedure C gives
The remainder ${\mathcal{C}}{\mathcal{C}}_\ell$ is
With the derivativedominant ranking $x\prec y\prec t$ , only one iteration is needed to complete the inversion, because $\boldsymbol{E}_{u}^{x}({\mathcal{C}}_0)=0$ . Then
and as ${\mathcal{C}}_0=\pi _u^x({\mathcal{C}}_0)$ , the calculation stops after updating, giving the output
This corresponds to the minimal result in Poole & Hereman [11].
Example 4.6. In every example so far, all iterations have used $u=u^\alpha$ , but not $u^\alpha _\boldsymbol{I}$ for $I\neq 0$ . The BBM equation from Example 2.1 illustrates the way that the process continues unhindered once the current $\mathcal{C}$ is independent of all elements of $[u^\alpha ]_x$ . The conservation law (2.6) is
it has no linear part or zerodegree terms. The ranking heuristic gives $x\prec t$ . The first iteration using $\boldsymbol{E}_{u}^{x}$ is similar to what we have seen so far, giving the following updates:
At this stage, ${\mathcal{C}}_1$ is independent of $[u]_x$ and equals $\pi ^x_{u_t}({\mathcal{C}}_1)$ . In the second iteration, the ranking requires us to apply $\boldsymbol{E}^x_{u_t}$ . This annihilates ${\mathcal{C}}_1$ ; consequently, ${\mathcal{C}}_1$ is a total derivative with respect to $x$ . Inverting this gives the final (minimal) result,
which was obtained by inspection in Example 2.1.
The Appendix lists inversions of various other divergences; in each case, the inversion produced by Procedure C and the ranking heuristic is minimal and takes very few iterations to complete.
4.6 Splitting a divergence using discrete symmetries
A given divergence may have discrete symmetries between various terms in its fully expanded form. If the divergence has very many terms that are connected by a particular discrete symmetry group, it can be worth splitting these into disjoint divergences that are mapped to one another by the group elements. Then it is only necessary to invert one of these divergences, using the symmetries to create inversions of the others without the need for much extra computation. However, to use Procedure C, it is necessary to check that all split terms are grouped into divergences; this check is done by using the criterion (2.3).
Polynomial divergences can first be split by degree, yielding divergences that are homogeneous in $[{\textbf{u}}]$ . Such splitting does not add significantly to the computation time, nor does it need to be checked using (2.3), which holds automatically. Splitting by degree can make it easy to identify terms that are linked by discrete symmetries, as illustrated by the following example.
Example 4.7. In Cartesian coordinates $(x,y)$ , the steady nondimensionalized von Kármán equations for a plate subject to a prescribed axisymmetric load function, $p(x^2+y^2)$ , are ${\mathcal{A}}_\ell =0$ , $\ell =1,2$ , where
Here $u$ is the displacement and $v$ is the Airy stress. This is a system of Euler–Lagrange equations. By Noether’s Theorem, the oneparameter Lie group of rotational symmetries yields the following conservation law (see Djondjorov and Vassilev [Reference Djondjorov, Vassilev, Bertram, Constanda and Struthers6]):
This conservation law has linear, quadratic and cubic terms, so $\mathcal{C}$ can be inverted by summing the inversions of each of the following divergences:
The quadratic terms have an obvious discrete symmetry, $\Gamma _1\;:\;(x,y,u,v)\mapsto (\!\!x,y,v,u)$ , which gives a splitting into two parts, each of which is a divergence:
where
Consequently, we can invert ${\mathcal{C}}_q$ by inverting $\overline{{\mathcal{C}}}_q$ and applying the symmetry $\Gamma _1$ .
Note that $\overline{{\mathcal{C}}}_q$ has a discrete symmetry, $\Gamma _2\;:\;(x,y,v)\mapsto (y,x,v)$ , which gives
Checking (2.3) shows that this is not a valid splitting into divergences, because
If necessary, split divergences can be inverted using different rankings. However, in this simple example, a single ranking works for all (nonlinear) parts. One tiebreak is needed: let $x\prec y$ . The ranking heuristic gives $v\prec u$ for the cubic terms; the variables $[u]$ are not relevant in the inversion of $\overline{{\mathcal{C}}}_q$ . Procedure C gives the following minimal inversions,
where
Applying $\Gamma _1$ to $\overline{{\mathcal{C}}}_q$ gives the following minimal inversion for $\mathcal{C}$ :
Note that the tiebreak $x\prec y$ causes the inversion to break the symmetry $\Gamma _2$ that is apparent in $\overline{{\mathcal{C}}}_q$ . As it turns out, this symmetry can be restored without increasing the overall number of terms, by adding components of a trivial conservation law (which are of the form (2.1)). It is an open question whether such preservation of symmetry and minimality is achievable in general.
In Djondjorov & Vassilev [Reference Djondjorov, Vassilev, Bertram, Constanda and Struthers6], the inversion of $\mathcal{C}$ for the special case $p=0$ has 62 terms, which are grouped according to their physical meaning. By contrast, the minimal inversion above has just 46 (resp. 48) terms when $p$ is zero (resp. nonzero), a considerable saving. Moreover, by exploiting the symmetry $\Gamma _1$ , only 28 (resp. 30) of these terms are determined using Procedure C. However, in seeking an efficient inversion, we have ignored the physics. It would be interesting to understand the extent to which the use of a minimal inversion obscures the underlying physics.
5. Concluding remarks
Partial Euler operators and partial scalings make it possible to invert divergences with respect to one independent variable at a time, a byproduct being that some contributions to other components are determined at each iteration step. Although each iteration involves a fair amount of computation, very few iterations are needed for many systems of interest.
Given the potential complexity of functions, it is unlikely that every divergence can be inverted, even in principle, by Procedure C. The question of how to prove or disprove this is open. In practice, products of mixed derivatives present the greatest challenge to concise inversion, although the option of reranking partway through the procedure enables a divideandconquer approach to be taken.
The focus of this work has been on inverting the total divergence operator Div. However, this immediately applies to expressions that can be recast as total divergences. For instance, for $p=3$ , the total curl ${\textbf{F}}=\mathrm{Curl}(\textbf{G})$ can be inverted by writing
where $\epsilon ^{ijk}$ is the Levi–Civita symbol, then inverting one component at a time and using the results at each stage to simplify the remainder of the calculation. Once $H^{ij}$ is known, the identity $G_l=\frac{1}{2}\epsilon _{ijl}H^{ij}$ recovers $\textbf{G}$ . Typically, a minimal inversion is achieved by using a different ranking for each $F^i$ , in accordance with the ranking heuristic. Here is a simple illustration of the general approach.
Example 5.1. In Cartesian coordinates, invert $(F^x,F^y,F^z)=\mathrm{Curl}(G_x,G_y,G_z)$ , where
Begin by using Procedure C with $y\prec z\prec x$ . In two iterations, this gives
With $(x,y,z)$ replacing the indices $(1,2,3)$ in (5.1), let
Therefore
One could invert this by a further iteration with the ranking $z\prec y$ , though it is inverted more quickly by the line integral formula, which gives $H^{yz}=u_yu_z$ . Finally,
Acknowledgments
I am grateful to Willy Hereman for discussions on homotopy operators and many helpful comments on a draft of this paper, and to the Centre International de Rencontres Mathématiques in Luminy for support and hospitality during the conference Symmetry and Computation. I would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the program Geometry, Compatibility and Structure Preservation in Computational Differential Equations; this provided the opportunity to develop much of the theory described above. I also thank the referees for their constructive comments.
Conflicts of interest
None.
A Appendix
Procedure C has been tested on the following conservation laws. In every case, the procedure coupled with the ranking heuristic yields a minimal inversion.
Kadomtsev–Petviashvili (KP) equation
There are two forms of the KP equation, depending on which of $\epsilon =\pm 1$ is chosen. In either case, it has a conservation law,
With the ranking $x\prec y\prec t$ , the procedure requires three iterations to obtain this inversion.
Potential Burgers equation (see Wolf et al. [Reference Wolf, Brand and Mohammadzadeh14])
The potential Burgers equation has a conservation law for each $f(x,t)$ such that $f_t+f_{xx}=0$ :
The inversion requires one iteration with $x\prec t$ , exchanging $f_t$ and $f_{xx}$ twice.
Zakharov–Kuznetzov equation (see Poole & Hereman [11])
The Zakharov–Kuznetzov equation is ${\mathcal{A}}=0$ , where
It has a conservation law
The fully expanded form has just eleven terms, and the inversion requires two iterations.
Shortpulse equation (see Brunelli [Reference Brunelli4])
The shortpulse equation, ${\mathcal{A}}=0$ , with
has the following conservation law involving a square root:
With $x\prec t$ , this can be inverted in one iteration.
Nonlinear Schrödinger equation
Splitting the field into its real and imaginary parts gives the system ${\mathcal{A}}_1=0,\ {\mathcal{A}}_2=0$ , with
One of the conservation laws is
With the ranking $u\prec v$ and $x\prec t$ , the procedure requires two iterations.
Itô equations (see Wolf [Reference Wolf13])
The equations are ${\mathcal{A}}_1=0,\ {\mathcal{A}}_2=0$ , with
This system has a rational conservation law,
With the ranking $u\prec v$ and $x\prec t$ , the procedure requires two iterations.
Navier–Stokes equations
The (constantdensity) twodimensional Navier–Stokes equations are ${\mathcal{A}}_\ell =0,\ \ell =1,2,3$ , where
This system has a family of conservation laws involving two arbitrary functions, $f(t)$ and $g(t)$ , namely
Procedure C consists of a linear inversion and two further iterations, using the ranking $x\prec y\prec t$ . The threedimensional Navier–Stokes equations have a similar conservation law, which requires a linear inversion and three further iterations.
Procedure C and the ranking heuristic have also been tested on some divergences not arising from conservation laws, with high order or complexity. Again, the output in each case is a minimal inversion.
Highorder derivatives
The divergence is
With the ranking $x\prec y$ , the procedure requires three iterations. The other ranking, $y\prec x$ , also produces a minimal inversion after three iterations; it is higher order in $x$ but lower order in $y$ :
Highorder derivatives and explicit dependence
The divergence is
With the ranking $x\prec y\prec t$ , the procedure requires one iteration. This illustrates the value of the second criterion for ranking independent variables; if $t$ is ranked lower than $x$ and $y$ , the procedure fails at the first check.
Exponential dependence
The divergence is
With the ranking $x\prec y$ , the procedure requires one iteration.
Trigonometric dependence
The divergence is
With the ranking $x\prec y$ , the procedure requires one iteration.