Group Foliation of Equations in Geophysical Fluid Dynamics

The method of group foliation can be used to construct solutions to a system of partial differential equations that, as opposed to Lie’s method of symmetry reduction, are not invariant under any symmetry of the equations. The classical approach is based on foliating the space of solutions into orbits of the given symmetry group action, resulting in rewriting the equations as a pair of systems, the so-called automorphic and resolvent systems, involving the differential invariants of the symmetry group, while a more modern approach utilizes a reduction process for an exterior differential system associated with the equations. In each method solutions to the reduced equations are then used to reconstruct solutions to the original equations. We present an application of the two techniques to the one-dimensional Korteweg-de Vries equation and the two-dimensional Flierl-Petviashvili (FP) equation. An exact analytical solution is found for the radial FP equation, although it does not appear to be of direct geophysical interest.

1. Introduction. The main theoretical objects of study in geophysical fluid dynamics are systems of partial differential equations describing the time evolution and thermodynamics of rotating stratified fluids [23,28]. Finding exact analytical solutions to these nonlinear partial differential equations (PDEs) can be exceptionally difficult, and increasingly physicists have instead turned to numerical solution methods. The rapid growth of computing power in the last few decades has made numerical methods more compelling, but perhaps more importantly, the computational tools now available allow many partial differential equations to be solved using only a few basic methods. Numerical results, however, can not provide the same confidence in the understanding of the underlying physics as analytical solutions do. Ideally, standard analytical methods could be applied towards identifying solutions to nonlinear PDEs with the same level of robustness as current computational techniques. Unfortunately, in sharp contrast with their linear counterparts, very few general analytical techniques exist for solving nonlinear PDEs.
One general approach that is applicable to many nonlinear PDEs is the method of symmetry reduction originally proposed by Sophus Lie. In Lie's method the point symmetry group of a system of PDEs is first computed algorithmically by solving the linear determining equations for symmetry generators. Any subgroup of the full symmetry group can then be used to reduce the system in order to construct solutions to the original equations that are invariant under the action of the subgroup [17]. This method is used implicitly by physicists when identifying radial solutions (rotation invariance) or traveling wave solutions (Galilean invariance), but also explicitly to reveal more complicated symmetry groups and corresponding invariant solutions (see, e.g., [27]). An unfortunate consequence of symmetry reduction is to eliminate any solution not invariant under the given group, which severally restricts the scope of available solutions and limits the applications to the physical sciences.
A more general, but less explored, technique to the analytical integration of PDEs is the method group foliation. The classical approach developed by Vessiot [30] can be applied to the same large class of linear and nonlinear PDEs as symmetry reduction, but it does not suffer from the limitations imposed by group invariance. Group foliation relies on a symmetry group's action on the space of solutions to the differential equation. Because the action of symmetry group transforms solutions to solutions, an individual solution can be transformed into other solutions along the same orbit by applying the group action. The method of group foliation provides a reduction of the original system of PDEs by constructing a set of equations identifying the individual orbits. Solutions found in the space of the reduced coordinates can then be returned to the full space by solving another system of PDEs. Importantly, in contrast with symmetry reduction, Vessiot's method, at least in principle, can be used to construct any solution of the original system from the reduced equations.
Although the method of group foliation was first put forth more than a century ago, renewed interest in the technique has only recently begun to develop. What might be described as the classical method of group foliation can be found in modern form in [19]. Here, a complete set of differential invariants of the differential equations' symmetry group are first computed and these are then chosen as the new coordinates with which to rewrite the equations. The resulting system is augmented by including differential relations between the invariants, known as syzygies. Applications of this method are still a rarity [16,2], but owing to recent advances in the complete algorithmic description of a symmetry group's differential invariants and their syzygies [18,8], the technique has become more readily applicable.
More recent work [4,13,14,25], recasts a system of PDEs as an exterior differential system (EDS) which is reduced under the action of the symmetry group of the equations; a process we call the EDS method of group foliation. The EDS encodes the differential equation as a set of differential forms known as contact forms that are defined on a subspace determined by the differential equations. This set of contact forms is reduced by including only the semi-basic forms, which are differential forms in the EDS annihilated by the infinitesimal generators of the symmetry group action. Solutions to the original system are then constructed from those of the reduced system via integrating a system of differential equations of generalized Lie type for the group parameters. Applications of this technique for finite dimensional symmetry groups can be found in [4], while the infinite dimensional case of pseudogroups is treated in [25].
Here we provide an introduction to both the classical and the EDS methods of group foliation in the context of geophysical fluid dynamics (GFD). Our goal is to introduce these techniques to the GFD community, by illustrating their application to familiar equations, and to highlight some of the challenges that remain.
2. Geophysical fluid equations. The geophysical fluid dynamical equations that we consider here are the one-dimensional Korteweg-de Vries (KdV) equation and various forms of the Flierl-Petviashivili (FP) equation. We are interested in the Flierl-Petviashivili equation in particular as an extension to quasi-geostrophic theory. Satellite altimetry observations of oceanic mesoscale features were initially compared to attributes predicted for linear Rossby waves [9], but higher resolution data now shows that the observed features are dominated by a high degree of nonlinearity and remain coherent structures for long durations [10]. In particular, the ratio of the observed eddy heights to the equivalent depth of the first baroclinic mode strongly violate the small amplitude assumptions made in quasi-geostrophic theory. The Flierl-Petviashvili equation extends quasi-geostrophy by allowing for larger amplitude features by introducing an additional height nonlinearity.
The FP equation (1) arises in several different contexts including the Great Red Spot on Jupiter [24], as well as modifications to quasi-geostrophy with an exterior mean shear flow [15], finite amplitude height changes [3,7], and thermobaricity of the equation of state [11]. In a general non-dimensionalized form, the equation can be written with only one constant β −1 as where φ(x, y, t) is a two-dimensional time-dependent stream function, J(a, b) = a x b y − a y b x is the 2-dimensional Jacobian, and β −1 = O(100) for oceanic values. This differs from the quasi-geostrophic potential vorticity equation [23,28] by the addition of the φ 2 term. No general analytical solutions to equation (1) are known, although the equation has been modeled numerically with β −1 = 1 [29]. The FP equation (1) can be reduced into a number of physically interesting forms. By assuming that φ takes the form of a scaled traveling wave such that φ(x, y, t) = λ 2 ψ (λ(x − ct), λy) where λ = 1 + c −1 , solutions to the reduced twodimensional equation, also satisfy the FP equation (1). With the additional assumption of radial symmetry, equation (2) is further reduced to ∂ ∂r Solutions to the radially symmetric equation (3) have been found numerically [15] and by approximation [24,5], but as of yet, no exact analytical solutions seem to have appeared in the literature. Finally, when equation (2) is restricted to only one dimension [15,26] by positing ψ(x, y) = u(x), the equation reduces to the one-dimensional traveling wave form of the Korteweg-de Vries equation, We will apply the method of group foliation to the three different equations (2), (3) and (4), beginning by solving the KdV equation (4) using both the classical and EDS methods of group foliation. The KdV equation reduces from a third order nonlinear ordinary differential equation to a first order ordinary differential equation on the quotient manifold. For the radial FP equation (3), reduction based on the scaling symmetry turns out to result in the obvious integration, but the integrated equations, for a particular value of the constant of integration, admits a symmetry that can used to construct an exact analytical solution. Finally, we consider the more complicated 2-D FP equation (2) in polar coordinates and reduce it by the EDS method under the rotational symmetry. We conclude with some remarks on the challenges of applying the EDS method to the full symmetry group of the FP equation (1).

Korteweg-de Vries equation.
3.1. Nonlinear wave solution. The Korteweg-de Vries equation is a classical example of a completely integrable system and has been thoroughly analyzed [1]. The physically interesting (real and bounded) cnoidal solution of (4) can be written in general form as The constant α scales the length and amplitude, x 0 translates the origin, and k varies the amplitude and the modulus of the elliptic function cn.
In terms of the original wave function φ(x, t), these solutions represent a train of traveling waves with amplitude dependent speed c = α 2 − 1 −1 . The well-known sech 2 solitary wave (soliton) solution is recovered from (5) in the limit k 2 → 1, for which the distance between the periodic waves is infinite.
3.2. Point symmetries. Central to both methods of group foliation is the existence of nontrivial symmetries for the differential equation of study. A basic introduction to analyzing point symmetries of differential equations can be found in e.g. [17], but we will briefly introduce here the concepts necessary for the basic computations.
A differential equation can be thought of as the zero set of some function of independent and dependent variables, and the derivatives of the dependent variables. In the case of one independent and dependent variable (x, u), the equation ∆(x, u, u x , u xx , ...) = 0 is simply an ordinary differential equation. An example of this is provided by the KdV equation (4) with ∆ given by We denote the total space of independent and dependent variables by E = {(x, u)} and the associated third order jet space by J 3 = {(x, u, u x , u xx , u xxx )}. In general, the order n of the jet space J n reflects the highest order derivative appearing in the equations.
By definition, the group of point symmetries of a system of differential equations consist of a local Lie group of transformations acting the total space E that map any solution of the equations to a new solution. More concretely, let u = f (x) be a solution to the the KdV equation (4) and G be a group of transformations acting on E = (x, u). Then G is a symmetry group of the equation provided that for each g ∈ G it is true that the transformed functionũ = g · f (x) is also a solution to the KdV equation (4).
Central to the algorithm for computing the point symmetries of a differential equation is the correspondence between finite transformations and their infinitesimal generators. For example, a one-parameter family of finite transformations, such as the rotations (x, u) → (x cos θ − u sin θ, x sin θ + u cos θ), θ ∈ R, define a curve through each point in the total space E, in this case a series of concentric circles. The corresponding infinitesimal transformation is the vector field everywhere tangent to these curves, namely w = −u∂ x + x∂ u for the group of rotations. Conversely, a vector field v = ξ(x, u)∂ x + η(x, u)∂ u defined on the total space E = {(x, u)} gives rise to a family of finite transformations obtained by exponentiating the generator, (x,ũ) = exp[ǫv](x, u), i.e., by finding the flow of v by solving a system of ODEs.
The next step is to determine the action of an infinitesimal transformation v on the derivatives of the dependent variable, a process known as the prolongation of v to J n . The prolongation formula can be derived by first finding the action of finite transformations on the derivative variables by an application of the chain rule and then computing the associated infinitesimal generator [17]. In this paper we will only treat differential equations involving one unknown function, so it suffices to consider a vector field on the total space The prolongation of v to the n-th order jet space J n can then be written as with the coefficients φ J given by Here J = (j 1 , j 2 , . . . , j n ) stands for a multi-index of integers and D J = D j1 D j2 · · · D jn denotes the repeated application of the coordinate total derivative operators The point symmetries arising from a local Lie group action can now identified by solving the determining equations for the infinitesimal generators of the group action. In the case of the KdV equation, these are found by prolonging a generic vector field v = ξ(x, u)∂ x +η(x, u)∂u on the total space E = {(x, u)} to the jet space J 3 = {(x, u, u x , u xx , u xxx )} and then demanding that the infinitesimal determining equations are satisfied. One thus finds two independent point symmetries with the infinitesimal generators corresponding to a translation and a scaling. By exponentiating the generators with (x,ũ) = exp[ǫv](x, u) we recover the corresponding finite transformations As a symmetry of the differential equation, the scaling transformation in (10) has the property that if u(x) is a solution to the KdV equation, then so isũ ǫ (x) = 1 − e −2ǫ (1 − u(e −ǫ x)). In terms of the solution (5) the group action simply rescales α. Indeed, we find that which changes the value of α to αe −ǫ , where ǫ ∈ R. It is in this sense that, given fixed values of x 0 and k, the action of the scaling symmetry traces out an orbit of solutions belonging to the same equivalence class. Simply by finding a single solution at any point along the orbit, say when α = 1, we can recover the entire equivalence class of solutions by applying the scaling transformations. In this sense the entire space of solutions can be partitioned into orbits, the equivalence classes, and the quotient manifold, denoted Sol/G, is the space of these equivalence classes.
The central idea of the method of group foliation is to exploit this partitioning afforded to us by the structure of the equation's symmetry group, and, roughly speaking, rewrite the equation on a certain quotient manifold determined by its point symmetries. Specifically, for the KdV equation, we will use v x to translate x 0 and v s to scale α. We should expect then that we will reduce the KdV equation from a third order ODE to a first order ODE, with the single constant of integration representing the modulus k of equation (5). 3.3. Classical method of group foliation. In the first step, the classical method of group foliation requires finding the automorphic system, which describes the orbits of the symmetry group action on the solution space Sol. The automorphic system is charaterized by the property that all of its solutions can be obtained from a fixed solution by symmetry transformations. As is expounded in [19], the automorphic system is described in terms of differential invariants of the symmetry group, which are split into new independent and dependent variables. For the KdV equation (4) this is the only equation, but in general, the automorphic system must be augmented with the resolving system, which encodes the integrability conditions of the former by way of syzygies among the differential invariants.
Differential invariants are functions of the jet space coordinates that are unaffected by the action of the group. Thus differential invariants ζ(x, u, u x , ...) of the whole group generated by the two vector fields v x and v s must vanish when acted upon by the prolongations of either vector field. The prolongations take the following form pr (2) pr (2) and together the prolonged group action maps coordinates by Invariance under the first vector field v x implies that differential invariants are independent of x. Consequently, the lowest order invariant can be found by looking for the x independent characteristics of the system pr (2) v s ζ = 0. The lowest order invariant is therefore found by integrating Appealing to our physical sense, using the differential invariant ζ poses a problem because 1−u may take on negative values in which case ζ will become imaginary. To resolve this, we simply use the squared value and define the lowest order differential invariant as The next order invariant is similarly found and is given by To verify that y and w are indeed differential invariants, simply apply the group action (14) to see that the group parameters λ and ǫ cancel out, as expected. When writing equation (5) in terms of these invariants, neither x 0 nor α will play a role.
As the next step in the algorithm we choose y as the new independent and w as the new dependent variable. We still need to compute a third order invariant to rewrite the KdV equation (4) in terms of the new variables. As is well known (see. e.g., [17], proposition 2.53), the derivative yields a new differential invariant. For simplicity, we will use w 1 = u xxx (1 − u) −5/2 as our third order differential invariant so that in terms of w and y, w 1 = dw dy 2wy 1/2 + 3y 3/2 − 2wy 1/2 .
Now one can easily see that the KdV equation becomes Writing this using y and w only we see that Equation (21) can be written in a standard form by letting a = 2w + 3y and z = 5y, which yields where a and z denote the differential invariants Equation (22) is an Abel equation of the second kind [31], a first order ODE, exactly as anticipated. Any solution to this equation yields a family of solutions to the KdV equation (4) when the original coordinates are reintroduced. Before examining solutions to (22), we will consider the alternative method of group foliation based on exterior differential systems.
3.4. EDS method of group foliation. The exterior differential systems (EDS) method of group foliation [4,13,14,25] requires rewriting the equation in terms of a differential ideal involving differential forms on the jet space J n . An exterior differential system (also called a differential ideal ) on J n is a subset I ⊂ Ω * (J n ) of the exterior algebra of all differential forms on J n satisfying the following properties: (i) if ω,ω ∈ I, then ω + fω ∈ I for any smooth function f on J n , 2) if ω ∈ I and η is any form in Ω * (J n ), then the wedge product ω ∧ η ∈ I, and 3) if ω ∈ I, then the exterior derivative dω ∈ I.
A differential equation ∆ : J n → R determines a submanifold S ∆ in the jet space J n , and so the contact ideal on J n , when restricted to S ∆ , encodes the differential equation and, as a consequence, its solutions. Recall that contact forms are characterized by the property that their pull-backs under the prolongation (i.e., under computing repeated derivatives) of any function f : X → U vanishes, where we have written the total space E = X × U as the product of the spaces X and U of the independent and dependent variables. For the KdV equation (4) it suffices to consider the contact ideal on the second order jet space J 2 with coordinates (x, u, u x , u xx ). Setting the contact forms ω 1 = du − u x dx and ω 2 = du x − u xx dx to zero requires that The KdV equation (4) is therefore described by the exterior differential system I generated by the forms ω 1 , ω 2 and the form ω 3 = du xx − (1 − u)u x dx; in short The EDS method of group foliation requires us to find a reduced exterior differential systemĪ on the quotient manifold associated with the symmetry group action. In practiceĪ is constructed by identifying the forms in I, the so-called semi-basic forms, that are annihilated by the prolonged infinitesimal generators of the group action and restricting these to the quotient manifold. As our prolonged infinitesimal generators are Γ (2) = span{pr (2) v x , pr (2) v s }, the semi-basic one-forms are determined by An arbitrary one-form in I takes the form β = β 1 ω 1 + β 2 ω 2 + β 3 ω 3 , where β 1 , β 2 , β 3 are smooth functions on J 2 . Now β is semi-basic provided that the two conditions are satisfied. We write equations (26) in matrix form as By performing a basic matrix row reduction to solve for the unknowns, we find one independent semi-basic form, which is given by In terms of coordinate differentials, we have where the dx terms cancel by virtue of (26). The next step is to choose a cross-section and pull back β. In order to write our equation in the same form as equation (22), we define our cross-section by where a, z denote the differential invariants (23). The du term vanishes on pullback, giving This form generates the idealĪ on the quotient manifold corresponding to the reduced equation, which is Equation (32) is identical to equation (22), which we found by Vessiot's method that directly employs differential invariants of the group action. While in the EDS computation it may initially appear as if we had avoided the entire complication of using differential invariants, these do indirectly play a role in our choice of the crosssection (30), as can be seen via a simple application of the moving frames algorithm [12]. Compare the cross-section of equation (30) with the group action (14) and use the u coordinate to solve for the group parameter, which gives ǫ = (1 − u) −1/2 . When this value of ǫ is substituted into the equations and the resulting expressions are solved for a and z, the specific choice of differential invariants is recovered, namely those in (23).

Reconstruction. Equation (32) possesses the parametric solution
see [31]. The constant C is related to the constant k of equation (5). All that remains is the reconstruction problem, that is, recovering solutions of the original equation (4) from (33). The classical group foliation and the EDS methods provide two different approaches. Starting with the EDS method, we look for solutions of the form where µ denotes the group action (14) and γ(t) is a symmetry transformation to be chosen so that the resulting s(t) vanishes on the ideal I defined in (24), see [4]. In the standard coordinates of J 2 , equation (34) becomes Here the unknown functions ξ(t), η(t) are related to the group parameters λ and ǫ in (14) by but they are now assumed to explicitly depend on the variable t. Next we simply compute the pullbacks of the generators (24) of the EDS I under s(t) and demand that the resulting forms vanish. This produces the equations The two equations can be solved simultaneously to yielḋ Requiring that the pullback of ω 3 under s(t) vanishes imposes no further constraints on the function ξ(t) or η(t). Reconstruction using Vessiot's method requires that we solve the expressions for the invariants (23) for u. Given that our solution (33) involves a parameter t, we will also construct x and u as functions of t. The expression for the z invariant in (23) produces the equationẋ We obtain the second equation by differentiating the z invariant and then using the expression for the a invariant in the resulting equation, which yieldṡ We note that the above equations are equivalent to the equations (36) for ξ(t), η(t), which we will integrate next. The group parameter α is recovered as the constant of integration from the equation for η(t) in (36). If we use t = 1 5 2s + 3 s + 1 to reparameterize equation (33), we find that η(s) = α 2 s. It then follows from equation (36) that dx ds which produces the translation group parameter x 0 upon integration in terms of an elliptic integral. To fully recover the solution (5), redefine the constant C = δ 3 1+δ and change variables by s = δ − q to see that The roots of the quadratic are r ± = 3 2 δ ± δ 2 δ−3 δ+1 . These are real if δ > 3, in which case r + > r − > 0, or if δ < −1, in which case 0 > r − > r + . For brevity, we only consider the first case, in which the integral can be written as This can be expressed in terms of a standard elliptic integral [6] as Solving (41) for q and writing the result in terms of u gives exactly the solution (5).

4.1.
Integration. The approximate and numerical solutions of the radial FP equation appearing in [15], [24], [5], were found by first integrating equation (3) to and then setting the constant of integration to zero, c = 0. An important difference between equations (3) and (43) for the purposes of group foliation is that for c = −1/2, equation (43) admits no point symmetries. If c = 1/2, equation (3) admits the same scaling symmetry as was found for the KdV equation, namely, with the finite action given by In this special case the method of group foliation can be applied to reduce the equation to a first order ODE.
unfortunately, does not appear to be of physical interest as the FP equation is derived under the assumption ψ ≪ 1, while equation (46) would appear to require that as r → ∞, ψ rr → 0 and ψ → 1. Nonetheless, we present this new solution in the appendix as it is the only known exact analytical solution of (3).

2-D FP equation in polar coordinates.
Applying the method of group foliation to an ordinary differential equation as in the preceding sections is essentially equivalent to using Lie's symmetry based method for integrating these equations (see, e.g., [17], section 2.5). While there is no direct extension of Lie's algorithm to equations with more than one independent variable, the two group foliation methods can be applied to systems of PDEs with the same basic reduction and reconstruction steps as are described above.
As an example, we consider the once integrated form 1 r ∂ ∂r r ∂ψ ∂r of the two-dimensional FP equation (2) written in polar coordinates The total space for equation (47) is E = (r, θ, ψ) and the associated jet space is given by J 2 = {(r, θ, ψ, ψ r , ψ θ , ψ rr , ψ rθ , ψ θθ )}. To apply the EDS approach of section 3, we encode the equation as the differential ideal on J 2 . Equation (47) admits translational symmetries in x and y, as well as a rotational symmetry which, in polar coordinates, is simply represented by a translation in the θ direction and is thus generated by the vector field v θ = ∂ θ . Here we will apply the EDS foliation algorithm with the rotational symmetry only.
Taking an arbitrary one-form in the ideal, β = β 1 ω +β 2 ω r +β 3 ω FP , and requiring that it vanish under the interior product with pr (2) v θ , we find the semi-basic forms which, when written in coordinate basis, become where ∆ = −r 2 ψ rr + 1 r ψ r − ψ + 1 2 ψ 2 . By choosing our cross-section with coordinates coordinates q, s, A, B, C and F as σ(q, s, A, B, C, F ) = (r = q, θ = 0, ψ = s, ψ r = A, ψ θ = B, ψ rr = C, ψ rθ = F ) , (53) we find that where Annihilating σ * ω 1 sb and σ * ω 2 sb and regarding (q, s) as the new independent variables, we obtain the system of equations After some algebraic manipulations, this system of equations together with (55) yield the system only involving the unknown functions A and H = B 2 /2. Thus, by applying the EDS method, we are able to reduce the elliptic equation (47) to a system of two quasi-linear equations (57), and, again, any solution of (47), including non-invariant ones, can, at least in principle, be reconstructed from the solutions of the system (57).
6. Conclusions. The classical and EDS methods of group foliation, when applied to the ODEs (3) and (4), provide useful techniques for reducing the order of the equations. In the case of the KdV equation (4), the general solution to the third order equation is recovered from the solutions of the Abel equation of the second kind, a first order ODE. The reduction of the radial FP equation (3) using the scaling symmetry simply results in the integration apparent by inspection. However, for a particular value of the constant of integration, the integrated equation (46) admits the same symmetry as its progenitor. This additional symmetry allows the equation to be integrated to yield an exact analytical solution to the radial FP equation. The application of the two methods of group foliation to PDEs, such as the 2-D FP equation (2), follows the same basic steps as in the case of ODEs. Various forms of the reduced equations can now be derived by choosing different symmetry subgroups defining the foliation or different cross-sections to the group orbits, allowing multiple approaches to solving the equations by the use of the same technique. As an example, the reduction of the 2-D FP equation written in polar coordinates (47) under the rotational symmetry is considered here. This results in rewriting the elliptic equation (47) as a system of two quasi-linear equations (57) that may be easier to analyze.
In the examples considered in this paper, the EDS method generally proves to be a simpler and more flexible approach than Vessiot's classical method. In the application of the EDS method to the KdV and the radial FP equation, the reduced equations on the quotient manifold were obtained without the need of having to compute the differential invariants of the symmetry group action, in contrast with the classical method that, in its original form, requires explicit expressions for the invariants and for their syzygies. Furthermore, choosing a different cross-section (and therefore a new set of differential invariants), the EDS method allows one to rewrite the reduced equations in new coordinates in an efficient manner. The reconstruction step in the EDS method is also straightforward to set up, and amounts to solving a system of differential equations of Lie type for the symmetry group parameters.
As the equation and its symmetry group increase in complexity, the advantages of the EDS method over Vessiot's approach become less pronounced. With the help of the moving frames method for computing differential invariants and their syzygies [18,8], the classical method of group foliation can be similarly used to write the reduced equations on the quotient manifold, identified as a cross section to the symmetry group orbits, without having to employ the explicit expressions for the differential invariants. For the symmetry subgroup of the 2-D FP equation (2) consisting of scalings and translations, this produces a system of three first order coupled differential equations. Here the main difficulty with applying the EDS method is in finding a managable basis for the semi-basic forms. However, as is implicit in the moving frames techniques applied in [25], if invariant contact forms are introduced in place of the coordinate differentials, the ensuing computations are vastly simplied and yield a basis for semi-basic forms that are comparably complex to the expressions one encounters in the application of Vessiot's method. Ultimately a combination of both the classical and EDS methods may be required in order to treat the more intricate examples.