Relationship between Unstable Point Symmetries and Higher-Order Approximate Symmetries of Differential Equations with a Small Parameter

: The framework of Baikov–Gazizov–Ibragimov approximate symmetries has proven useful for many examples where a small perturbation of an ordinary or partial differential equation (ODE, PDE) destroys its local exact symmetry group. For the perturbed model, some of the local symmetries of the unperturbed equation may (or may not) re-appear as approximate symmetries. Approximate symmetries are useful as a tool for systematic construction of approximate solutions. For algebraic and ﬁrst-order differential equations, to every point symmetry of the unperturbed equation, there corresponds an approximate point symmetry of the perturbed equation. For second and higher-order ODEs, this is not the case: a point symmetry of the original ODE may be unstable, that is, not have an analogue in the approximate point symmetry classiﬁcation of the perturbed ODE. We show that such unstable point symmetries correspond to higher-order approximate symmetries of the perturbed ODE and can be systematically computed. Multiple examples of computations of exact and approximate point and local symmetries are presented, with two detailed examples that include a fourth-order nonlinear Boussinesq equation reduction. Examples of the use of higher-order approximate symmetries and approximate integrating factors to obtain approximate solutions of higher-order ODEs are provided.


Introduction
A symmetry of a system of algebraic or differential equations is a transformation that maps solutions of the system to other solutions. Dating back to the works of Sophus Lie in the nineteenth century, symmetry ideas have seen significant development over the last century, relating to symmetry reduction and solutions of differential equations, integrating factors, conserved quantities and local conservation laws, Hamiltonian and Lagrangian structures, integrability, nonlocal extensions, invertible and non-invertible mappings between different classes of differential equations, and more (see, e.g., [1,2] and the references therein).
Let the term perturbed equations denote equations differing from some canonical or otherwise well-understood model by extra term(s) involving a small parameter. This small perturbation disturbs the local Lie symmetry properties of the unperturbed equations. Several approximate Lie symmetry methods have been developed to study symmetry properties of perturbed models, and relate and compare them to symmetry structure of the unperturbed equations. An approximate symmetry method (referred to here as the BGI method) was introduced by Baikov, Gazizov, and Ibragimov [3][4][5], where the approximate symmetry generator is expanded in a perturbation series. This approach preserves the Lie group structure, in particular, a commutator of two approximate symmetry generators yields an approximate symmetry generator [6]. Using the BGI framework, approximate symmetries, first integrals, and approximate solutions have been constructed for a number of models involving ordinary and partial differential equations (ODE and PDE) (e.g., [7][8][9].) A different approach to approximate symmetries, developed by Fushchich and Shtelen [10], combines a perturbation technique with the symmetry group method by expanding the dependent variables in a Taylor series in the small parameter and approximately replacing the original equations by a system of equations that are coefficients at different powers of the parameter. The classical Lie symmetry method is applied to obtain symmetries of the latter system. Using this method, approximate symmetries and approximate solutions have been found for some PDE models [11,12]. The BGI and Fushchich-Shtelen approaches are not equivalent. They have been compared and used to obtain approximate symmetries and approximate solutions for several PDE models [13][14][15][16]. In [17], approximate solutions of a singularly perturbed Boussinesq PDE were found using the approximate Fushchich-Shtelen symmetries and using a proposed method that is not based on the Lie symmetries. Burde [18] developed a new approach for approximate symmetries by constructing equations that could be reduced by exact transformations to an unperturbed equation and would, at the same time, coincide approximately with the perturbed equation.
Contact and higher-order exact symmetries can be used to construct solutions for ordinary differential equations (e.g., [19]). In [20], it was shown how integrating factors for linear and nonlinear ordinary differential equations can be determined. A perturbation method based on integrating factors was developed for a system of regularly perturbed first-order ODEs [21].
In this paper, we follow the BGI approximate symmetry framework for algebraic equations and ODEs. While BGI approximate symmetries have been found for many models, including ODEs and PDEs (e.g., [3][4][5][6][7]14,22]), it has not been clarified under what conditions point or local symmetries of exact equations can become unstable, disappearing from the classification of approximate symmetries of a perturbed system of the same differential order, and what form they, therefore, take.
In the sections below, we prove that for single algebraic and first-order differential equations, all point symmetries are stable. For second and higher-order ODEs, we show that a point or local symmetry of the unperturbed equation usually yields a higher-order (generally, of order n − 1) symmetry of the perturbed model.
The paper is organized as follows. In Section 2, we briefly review the necessary notation for the Lie group of transformations, infinitesimal transformations and determining equations for finding the exact symmetries, and provide an introduction to the theory of approximate transformations and approximate symmetries of the perturbed equations in the sense of Baikov, Gazizov, and Ibragimov. In Section 3, we consider exact and approximate point symmetries of algebraic and ordinary differential equations with a small parameter, and provide a relation between exact and approximate symmetries of the original and perturbed algebraic and first-order ordinary differential equations. We investigate the BGI approximate point symmetries of the perturbed higher-order ODEs. Point symmetries of the unperturbed equation may indeed disappear from the classification of approximate point symmetries of the perturbed model, and conditions for that are presented.
In Section 4, we consider point and higher-order local exact and approximate symmetries of second and higher-order ODEs in evolutionary form, and present a systematic way (Theorem 2) to find approximate symmetry components for approximate symmetries that correspond to every point and local symmetry of the unperturbed equation. Relations between exact and approximate symmetries are considered in detail for two examples, including a nonlinearly perturbed second-order ODE, and a fourth-order ODE arising as a traveling wave reduction of the Boussinesq partial differential equation modeling shallow water wave propagation.
Finally, in Section 5, we determine approximate integrating factors of perturbed firstorder ODEs using approximate point symmetries. We find the determining equations of approximate integrating factors, and show how these determining equations and higher-order approximate symmetries can be used to obtain approximate solutions of a perturbed Boussinesq ODE. A brief discussion is offered in Section 6.
In addition to providing a complete answer to the question of stability of point and local symmetries of unperturbed ODEs vs. their perturbed versions with a small parameter, the main value of this contribution lies in new detailed examples of computation and comparison of exact and approximate symmetry structures of multiple ODEs, and the use of point and higher-order approximate symmetries to calculate closed-form approximate solutions of such perturbed models.

Lie Groups of Exact and Approximate Point and Local Symmetries
We denote a general system of N algebraic or differential equations by and its first-order perturbation in terms of a small parameter by x = (x 1 , x 2 , . . . , x n ), n ≥ 1, and v = (v 1 , v 2 , . . . , v m ), m ≥ 1, are respectively independent and dependent variables, and ∂ q v denotes all qth-order derivatives of all components of v.

Exact and Approximate Transformation Groups
A one-parameter Lie group of transformations with the group parameter a, and the corresponding infinitesimal generator where summation in repeated indices is assumed, is a point symmetry of the system (1) when, for each σ = 1, 2, . . . , N, that is, (5) holds on solutions of (1) (e.g., [1,2]). The evolutionary (characteristic) form of the Lie group of transformations (3) is the one-parameter family of transformations with the evolutionary infinitesimal generator where Higher-order local transformations generalize (6) by allowing the infinitesimal components ζ[v] to depend more generally on derivatives of v, including higher-order derivatives. Example 1. The classical example (whose various perturbed versions we will use below) is that of a second-order ODE with a maximal Lie group of point symmetries, (Here and below, we use primes to denote ordinary derivatives.) The computation of the prolongation of X 0 (4) to the second order and the solution of determining Equation (5) yields the general point symmetry components (see, e.g., [1]) where C i are arbitrary constants. The resulting eight-parameter Lie group of point symmetries of (8) is spanned by the generators For a general model (2) involving a small parameter , exact point and local symmetry generators have the form Solving the determining Equation (5), one finds exact symmetries of (2) holding for an arbitrary . It is commonly the case that, due to the perturbation term, some (or even all) point and/or local symmetries of the unperturbed Equation (1) disappear from the local symmetry classification of the perturbed model (2).

Example 2.
Consider an ODE y = (y ) −1 , which is a perturbed version of (8). It is straightforward to show that the only point symmetries of (8) that are also point symmetries of (12), holding for an arbitrary , are the translations Approximate symmetries have been developed a tool to seek additional symmetry structure of perturbed models. For Equation (2) with a small parameter , Baikov-Gazizov-Ibragimov (BGI) approximate point symmetries [22,23] are defined by approximate symmetry generators and similarly, local approximate BGI transformations correspond to generators in evolutionary form given byX The approximate invariance condition yields determining equations for the components ξ i j , η µ j of first-order BGI approximate symmetries. The solution of the determining Equation (16) can be subdivided in the following steps:

1.
Compute an exact point/local symmetry generator X 0 of the unperturbed Equation (1) using determining Equation (5) for exact local or point symmetries.

2.
Find the corresponding first-order deformation (the part X 1 of the generator (14)) using the equation where G is obtained from the coefficients of in , σ = 1, . . . , N.
The following theorem holds.  (16) that any symmetry generator X 0 of an unperturbed system (1) yields an approximate symmetry X = X 0 of the perturbed system (2). We call such approximate generators trivial approximate symmetries of (2).

Remark 2.
Similarly to the above-described procedure, one can consider higher-order expansions of both the perturbed Equation (2) and symmetry generators in terms of the small parameter .

Stable and Unstable Symmetries in the BGI Framework
The converse to Theorem 1 does not hold: not every point or local symmetry of a given model (1) yields an approximate BGI symmetry of its perturbed version (2). In [22], an exact point (or local) symmetry of the unperturbed Equation (1) with the generator (4) (or (7)) is called stable if there exists a point (local) approximate generator (14) or respectively (15) that is a BGI approximate symmetry of the perturbed Equation (2). If all symmetries of the Equation (1) are stable, the perturbed Equation (2) are said to inherit the symmetries of the unperturbed equations. Its unperturbed version y = 0 has eight exact point symmetries given by (10). The approximate BGI symmetry generator of (12) can be sought in the form where X 0 is an exact symmetry generator of the unperturbed ODE. The determining Equation (16) for approximate symmetries yield where ξ 0 , η 0 are exact symmetry components (9) computed in Example 1. The determining Equation (17) splits into a system of PDEs for ξ 1 , η 1 with the solutions and the additional conditions: 3ξ 0 which is a different perturbed version of (8). The determining Equation (16) for approximate symmetries of (21) yields where ξ 0 , η 0 are the unperturbed symmetry components (9). Clearly, Equation (22) splits into a system of PDEs in ξ 1 , η 1 with no change on ξ 0 , η 0 . Consequently, the perturbed ODE (21) admits 16 approximate symmetries given by All exact symmetries (10) of the unperturbed ODE (8) are inherited by the approximate symmetries (23), and thus are stable by definition.

Algebraic and First-Order Differential Equations
The relation between exact and approximate point symmetries of algebraic equations and fist-order ODEs is quite simple. In summary, to every exact Lie point symmetry of an unperturbed equation, there correspond: • an infinite set of exact Lie point symmetries of the perturbed equation; and • an infinite set of approximate BGI point symmetries of the perturbed equation.
It follows that all point symmetries of algebraic systems and first-order ODEs are stable in the BGI approximate symmetry sense.
The above statement is the result of both algebraic equations and ODEs having infinite sets of point symmetries, in both the classical Lie and approximate BGI framework. In particular, for algebraic equations, let define a family of surfaces (curves) in R n , with F 0 being a scalar or vector function of m components, 1 ≤ m < n. Let also denote a perturbation of (24). Suppose the point symmetry generator preserves the solution set of (24) in the sense that For the linear system (27), the dimension of the solution space is d = n − rank (DF 0 /Dx) ≥ 1, that is, the point symmetry generator (26) is parameterized by d arbitrary functions. In the same manner, assuming that rank (DF 0 /Dx) = rank (DF/Dx), for the point symmetry generator Y = η i (x; )∂/∂x i of the perturbed system (25), the symmetry determining equations YF j (x) = 0, j = 1, . . . , m yields infinite point symmetries involving d arbitrary functions, with lim →0 η i = ξ 0i . For BGI approximate point symmetries admitted by the family of perturbed surfaces (25), the generator has the form The determining Equation (16) yield the conditions on the first-order infinitesimals ξ 1i For each point symmetry (26) of the unperturbed Equation (24) given by a set of the infinitesimal components, the determining Equation (29) have multiple nontrivial solutions parameterized by d = n − rank (DF 0 /Dx) ≥ 1 functions; thus, to every point symmetry of (24), there corresponds a d-dimensional set of approximate BGI point symmetries (28), again satisfying lim →0 X = X 0 . For a first-order ODE let denote an exact point symmetry generator admitted by (30). Exact point symmetry components (ξ 0 (x, y), η 0 (x, y)) of the ODE (30) satisfy the determining Equation (5) η 0 that is, for a fixed arbitrary function ξ 0 (x, y), a linear non-homogeneous PDE on η 0 (x, y), which has infinitely many solutions, corresponding to infinite point symmetries of the first-order ODE (30). In particular, for an arbitrary For a perturbed version of the ODE (30) the exact symmetry generator has the form Since the right-hand side of the ODE (33) is just another function of x, y, by the same reason as above, the perturbed ODE (33) has an infinite set of exact point symmetries with generator (34), arising as solutions of the determining equation Again, for an arbitrary ξ(x, y; ) analytic in , one can find η(x, y; ) analytic in . Consequently, when = 0, each symmetry (34) of the perturbed ODE (33) reduces to the exact point symmetry (31) of the unperturbed ODE (30). For the perturbed ODE (33), one can also seek a BGI approximate symmetry generator in the form Applying the approximate symmetry determining Equation (16), one has which is a linear nonhomogeneous PDE in two unknowns ξ 1 and η 1 . Hence, for an arbitrary fixed ξ 0 (x, y), η 0 (x, y), and ξ 1 (x, y), an infinite set of solutions η 1 (x, y) can be found, corresponding to an infinite set of BGI approximate symmetries (36) of the perturbed ODE (33), corresponding to any exact symmetry (31) of the unperturbed ODE (30).

Second and Higher-Order ODEs
The situation with stability of point symmetries of second and higher-order ODEs in the BGI framework is significantly different: determining equations on BGI approximate symmetry components may (or may not) contain conditions on the exact symmetry components, which can lead to unstable point symmetries. Consider the unperturbed higher-order ODE and its perturbed version Perturbed ODEs generally have fewer exact point symmetries than their unperturbed versions; Example 2 for the ODE y = (y ) −1 illustrates this trend. The exact and BGI approximate symmetry generators for (38), (39) are given by (31) and (36). To find the BGI approximate symmetries of the perturbed ODE (39), we apply the approximate invariance condition (16). In the zeroth order in , they are the same as (5) for exact point symmetries.
At the first order in , one has equivalent to Note that η 0 (n) is linear in y (n) , and satisfies the equation Hence the general form for the determining equation for approximate symmetries of the perturbed ODE (39) is The additional determining Equation (42) is a PDE on the BGI first-order perturbation point symmetry components (ξ 1 , η 1 ) which are functions of x, y. After replacing y (n) by f 0 (x, y, y , . . . , y (n−1) ), and splitting with respect to different combinations of y , . . . , y (n−1) on which (ξ 1 , η 1 ), one obtains a set of linear homogeneous PDEs. These involve the unknown BGI perturbation components ξ 1 , η 1 as well as the exact point symmetry components (ξ 0 , η 0 ) obtained from (5) in the previous step.
In particular, depending on the form of f 0 and f 1 in (39), these split determining equations may contain additional conditions on the exact point symmetry components ξ 0 , η 0 . When that is the case, some exact point symmetries of the unperturbed ODE (38) may disappear from the approximate symmetry classification of the perturbed ODE (39), thus becoming unstable (see Example 3). If the additional determining Equation (42) contains no restrictions on the exact point symmetry components ξ 0 , η 0 , all symmetries of the unperturbed ODE remain stable (Example 4).

Exact and Approximate Local Symmetries of Higher-Order ODEs
While for algebraic equations and first-order ODEs, every point symmetry of the unperturbed equation is stable and reappears in the BGI approximate symmetry classification, we have seen that for second-and higher-order ODEs, this is not the case: point symmetries of a second or higher-order ODEs may or may not be stable.
By analogy with ODE systems, for higher-order ODEs, it is natural to expect that the correct framework is provided by local (including higher-order) symmetries. Indeed, below we show that to every point or local symmetry of an unperturbed ODE of second or higher order, there corresponds a local BGI approximate symmetry of the perturbed ODE.

Exact Local Symmetries of the Unperturbed ODE
The infinitesimal generator of a local symmetry (6) admitted by an unperturbed ODE (38) has the formX with the infinitesimal component The nth prolongation of (43) is given bŷ The determining equations for the exact symmetries of the general ODE (38) arises from the invariance conditionX or in detail, If s = n − 1, Equation (46) is a linear homogeneous PDE for ζ 0 with independent variables x, y, y , . . . , y (n−1) . This PDE can be written in solved form for the highest derivative of ζ 0 by x, where all derivatives with respect to x appearing in the right-hand side of (47) are of lower order than those appearing on the left-hand side. It follows that when s = n − 1, the PDE (46) is solvable for ζ 0 for any "initial condition", and hence any given ODE of order n admits an infinite number of local symmetries of order n − 1, parameterized by solutions of the PDE (47).

Approximate Local Symmetries of the Perturbed ODEs
The higher-order approximate symmetry generator for the ODE (39) with a small parameter is given byX where ζ 0 [y] is given by (44), and The prolongation of this generator has the form with ζ 1 (j) = D j ζ 1 , j = 1, 2, . . . , n. To find the approximate symmetries of the perturbed ODE (39), we apply the determining equations for approximate symmetrieŝ First, one computes an exact local symmetry generator (43) of the unperturbed ODE (38). Then, the first-order deformationX 1 can be found from the equation where G is the coefficient of in The additional determining Equation (51) becomes D n ζ 1 When = n − 1, Equation (52) yields a linear nonhomogeneous PDE in ζ 1 which has a Cauchy-Kovalevskaya form with respect to the independent variable x, and its solutions that can be obtained by the method of characteristics. No additional restrictions on the unperturbed symmetry component ζ 0 arise. Hence, the following theorem holds. Theorem 2. For each exact point or local symmetry (43) of an unperturbed ODE (38), there is an approximate local symmetry (48) of the perturbed ODE (39), with the symmetry component ζ 1 being of differential order at most n − 1.
We now consider two examples in detail.

The First Detailed Example
For the second-order ODE (12) with a small parameter, y = (y ) −1 , we apply Theorem 2 to find approximate symmetries of order n − 1 = 1 corresponding to unstable point symmetries of (12) (see Example 3). This ODE is a perturbed version of y = 0. In total, it admits 12 approximate point symmetries; this set does not include the following unstable point symmetries of y = 0: be the symmetry generator of the ODE y = 0 in evolutionary form. Therefore, ζ 0 has the form The eight point symmetries (10) of y = 0 have evolutionary formŝ be a local approximate symmetry generator admitted by the perturbed ODE (12) where ζ 0 is given by (55). The determining Equation (52) on ζ 1 requires that By a change of variables t = y − xy , ζ 1 (x, y, y ) = u(x, t), the homogeneous part of (58), becomes u xx = 0, with the general solution u(x, t) = R 1 (t) + xR 2 (t), where R 1 , R 2 are arbitrary functions of their arguments. Hence, the homogeneous PDE (59) has a general solution ζ 1 c = R 1 (y − xy ) + xR 2 (y − xy ). Now, let be a particular solution for the nonhomogeneous PDE (58). Substituting this particular solution into the Equation (58) yields the following system of PDEs R xx + 2Q xy + P yy = 2α 3 x + 4α 6 , 2R xy + Q yy = 0, R yy = 0.
The corresponding approximate symmetry of the perturbed ODE (12) is given bŷ which is exactly the evolutionary form of the approximate point symmetry X 11 in (19b).

Remark 3.
We note that, in the current example, one would obtain an infinite set of first-order approximate symmetries corresponding to each unstable point symmetry (53) of y = 0, if a more general form (60) of ζ 1 (x, y, y ) was used instead of the simplified ansatz (61). This, however, does not make such first-order approximate symmetries trivial; they can be used, for example, for construction of approximate solutions of the perturbed ODE (12) through mappings or approximate reduction of order (see Section 5 below).

The Second Detailed Example
In the following example, we compute exact point and local symmetries of the fourthorder Boussinesq differential equation [24,25] and discuss their stability. Consider a linear ODE y (4) + y = 0 (63) and its perturbed version, the Boussinesq ODE y (4) + y − 2yy + 2y 2 = 0.
The latter ODE can be obtained as a time-independent or a traveling wave reduction of the Boussinesq partial differential equation which was introduced in 1871 to describe the propagation of long waves in shallow water [26]. In this example, some point and local symmetries of the unperturbed ODE (63) are shown to correspond to third-order local approximate BGI symmetries of the perturbed ODE (64), as guaranteed by Theorem 2. The calculated approximate symmetries are used in the next section to illustrate the construction of an approximate solution of the perturbed Boussinesq Equation (64).

Exact Point Symmetries of (63); Approximate Point Symmetries of (64)
First, we seek exact point symmetries for (63) and approximate point symmetries for (64). Let be an exact point symmetry generator of the ODE (63). After the prolongation of X 0 to the fourth-order and applying the determining Equation (5), one finds involving six arbitrary constants. Consequently, the ODE (63) has a six-dimensional Lie algebra of point symmetry generators, spanned by Next, we proceed to find approximate point symmetries of the Boussinesq ODE (64). Let denote the approximate BGI symmetry generator admitted by (64), where X 0 is an exact symmetry generator (66) of the unperturbed ODE (63). The determining equation for approximate symmetries (42) yields The above system has the solution also involving six arbitrary constants. Specifically, the perturbed ODE (64) admits six trivial approximate symmetries X j = X 0 j , j = 1, 2, . . . , 6, corresponding to the free constants a 1 . . . , a 6 , where X 0 j are the exact point symmetries (68) of the unperturbed ODE (63) and two nontrivial approximate point symmetries It follows that the only two stable point symmetries of (63) are X 0 1 and X 0 2 , while the remaining ones X 0 j , j = 3, . . . , 6 in (68) are unstable.

Exact
Second-Order Local Symmetries of (63); Approximate Second-Order Local Symmetries of (64) We now extend the above analysis, seeking exact local symmetries admitted by (63) up to second-order, in the form Applying the determining Equation (46), one finds The above equation splits into system of PDEs. Solving this system gives ϕ 0 = k 3 y + k 2 + k 3 x + k 4 y + k 5 sin x + k 6 cos x + k 7 y + k 8 (y sin x + y cos x) +k 9 y 2 + y 2 + k 10 y 2 − y 2 cos x + 2y y sin x + k 11 y 2 + y 2 sin x + 2y y cos x +k 12 y 2y − x + 2y − xy 2 + k 13 2 sin x − x cos x y − x sin x + cos x y + 2y sin x +k 14 x sin x + 3 cos x y + 2 sin x − x cos x y + y cos x + k 15 (y sin x − y cos x), involving 15 arbitrary constants k j . Hence, the ODE (63) admits local symmetries These generators were numbered to match the point symmetry classification (68) of the unperturbed ODE (63). In particular, the generators V 1 , . . . , V 6 in (76) are evolutionary forms of the point symmetries (68). Now, we will find the approximate local symmetries for the perturbed ODE (64). Let be the local approximate symmetry generator admitted by the perturbed ODE (64) where ϕ 0 is given by Equation (75). Using the determining Equation (52), one obtains ϕ 1 = Q 1 (y) + y Q 2 (y) + a 3 x + a 4 y + a 5 sin x + a 6 cos x + a 7 y 2y − x + 2y − xy 2 +a 8 y 2 + y 2 + a 9 y 2 + y 2 sin x + 2y y cos x + a 10 y 2 − y 2 cos x + 2y y sin x +a 11 (y sin x − y cos x) + a 12 2 sin x − x cos x y − x sin x + cos x y + 2y sin x +a 13 (y sin x + y cos x) + a 14 x sin x + 3 cos x y + 2 sin x − x cos x y + y cos x −k 2 xy + k 3 2xy − 1 2 This set includes the evolutionary forms of the approximate point generators X 1 and X 2 of (72) in their evolutionary forms V 1 and V 2 . Moreover, V 3 is a second-order approximate symmetry of the perturbed ODE (64) corresponding to the unstable point symmetry X 3 0 in (68), and V 7 is an evolutionary form of the approximate point symmetry X 7 in (72).

Higher-Order Approximate Symmetries Corresponding to Unstable Point and Local
Symmetries of (63) be the evolutionary form of the exact point or local symmetry generator of the unperturbed ODE (63). Here, ζ 0 = ζ 0 (x, y, y ) for point symmetries (68), and ζ 0 = φ 0 (x, y, y , y ) for second-order local symmetries (76) of the unperturbed ODE (63). Following Theorem 2, for each unstable local symmetry V j 0 , j = 4, 5, 6, 8, . . . , 15 in (76) of the ODE (63), there is a corresponding higher-order approximate symmetry for the perturbed ODE (64) of the form First, we consider the unstable point symmetry X 4 0 in (68) (V 0 4 in (76)). Its evolutionary components is ζ 0 = y. The corresponding ζ 1 is any solution of the linear nonhomogeneous PDE A simple particular solution is given by One consequently obtainŝ as a third-order approximate symmetry for the perturbed ODE (64) corresponds to the unstable point symmetry X 4 0 , V 4 0 . In the same way, one can find a third-order approximate symmetry corresponding to each unstable point symmetry of (63) in (68) or unstable local symmetry in (76). Let V = ϕ 0 (x, y, y , y ) + φ 1 (x, y, y , y , y ) ∂ ∂y (81) be an approximate symmetry generator for the perturbed ODE (64) where ϕ 0 is given by the Equation (75) . From the determining Equation (52), one can find thatφ 1 corresponds to each local symmetry of (76). For example, consider the unstable local symmetry V 0 9 = y 2 + y 2 ∂/∂y. By substituting ϕ 0 = y 2 + y 2 into the determining Equation (52), one obtains The above equation has a particular solution given bŷ Hence, is a third-order local approximate symmetry of the Boussinesq ODE (64) corresponding to the exact local symmetry V 0 9 of the unperturbed Equation (63), which used to be unstable in the class of second-order local symmetries.

Reduction of Order and Approximately Invariant Solutions of Perturbed Differential Equations
In this section, we discuss approximate reduction techniques, including approximate integrating factors and approximate first integrals of perturbed differential equations, and the use of the higher-order approximate symmetries to find approximate solutions of perturbed ODEs.
This follows from taking ξ(x, y; ) = ξ 0 (x, y) + ξ 1 (x, y) + o( ) and η(x, y; ) = η 0 (x, y) + η 1 (x, y) + o( ), substituting these values into (86) and taking the Taylor expansion about = 0. Example 5. The first-order ODE y = y + xy (88) admits the approximate symmetry generator X = (1 + )y ∂/∂y. An approximate integrating factor for (88) has the form µ(x, y; ) = (1 − )/y. Using this integrating factor, one finds Hence, is a family of approximate solution curves for the perturbed ODE (88). Note that the first two terms of the Taylor expansion in of (90) agree with the first two terms of the Taylor expansion in of the exact solution of the ODE (88):
Equation (97) can be written in the form Hence, the Boussinesq ODE (97) reduces to the second-order ODE The general solution of (98) is unknown. An approximate solution can be constructed in the assumption of C 1 , C 2 = O( ). Let C 1 = c 1 , C 2 = c 2 ; then, the ODE (98) becomes Using the determining Equation (96), one can easily find that µ = y + (y − c 1 ) is an approximate integrating factor for the ODE (99). Multiplying this integrating factor by (99) yields and consequently an approximate first integral: Hence, the perturbed Boussinesq ODE (97) is reduced to the first-order ODE where c 1 , c 2 , c 3 are arbitrary constants. A series ansatz y(x; ) = y 0 (x) + y 1 (x) + o( ) into the ODE (100) leads to the system of ODEs (y 0 ) 2 + y 2 0 = 2c 2 3 , with the solutions Finally, a general approximate solution for the Boussinesq ODE (97) involving four arbitrary constants is obtained: y(x; ) = c 3 (sin x + cos x)
The differential functions ω k (x, y, y , . . . , y (k) ; ) = ω 0 k (x, y, y , . . . , y (k) ) + ω 1 k (x, y, y , . . . , y (k) ) + o( ), k = 1, . . . , n are approximate differential invariants for the ODE (102) ifX (k) ω k (x, y, y , . . . , y (k) ; ) = o( ). Note that ω 0 k are exact differential invariants for the unperturbed ODE (38). They arise as constant of integrations of the characteristic equations Then, the approximate differential invariant components ω 1 k are determined from the condition where H is a differential expression in terms of ω 1 k arising from the coefficient of in Example 7. The first example of using approximate differential invariants to reduce ODEs is rather elementary and is used here for illustration purposes. Consider the second-order ODE This ODE admits an approximate contact symmetry given bŷ We determine the invariants ω(x, y, y ; ) = ω 0 (x, y, y ) + ω 1 (x, y, y ) + o( ) satisfyingX (1) ω = o( ). Clearly, one invariant is x. Other invariants are determined by first finding ω 0 satisfyinĝ which has a general solution ω 0 (x, y, y ) = F(xy − y) based on the fundamental invariant xy − y. Let ω 0 (x, y, y ) = xy − y. Then, one finds that the first-order correction satisfies the inhomogeneous linear PDE The simplest particular solution is given by ω 1 (x, y, y ) = −(x 3 y 2 + y 3 )/3. Consequently, is an approximate invariant for the ODE (105); here, C 1 = const is a constant of integration. Thus, the ODE (105) approximately reduces to a first-order ODE. By substituting y(x; ) = y 0 (x) + y 1 (x) + o( ) into the ODE (106) and setting to zero coefficients at 0 and 1 , one finds a system of ODEs Its solution yields an approximate solution of the perturbed ODE (105) We note that the ODE (105) is solvable by separation of variables, which makes it easy to compare its general solution with the approximate solution (107). The general solution is The first two terms of its Taylor series with respect to indeed coincide with the approximate solution (107).
Substituting a series expansion y(x; ) = y 0 (x) + y 1 (x) + o( ) into the Equation (112) and equating the coefficients of 0 , 1 , we find y 0 = C 1 x + C 2 sin x + C 3 cos x + C 4 and y 1 = −h(x, y 0 , y 0 , y 0 , y 0 ). Hence, the approximate solution of the Boussinesq ODE (97) is given by The unperturbed ODE (63) with the initial conditions y(0) = 1, y (0) = 1, y (0) = −1, y (0) = −1 has a particular solution y(x) = sin x + cos x. (115) Using this particular solution and the corresponding different set of initial conditions in (114), one finds C 1 = 0, C 2 = 1, C 3 = 1, and C 4 = 0. This particular approximate solution (114) of the perturbed ODE (97) has the form In order to test the accuracy of the approximate solution (117), we convert the perturbed fourth-order ODE (97) into a system of four first-order ODEs, and compute numerical solutions of the resulting system with the initial conditions (116) using the Matlab native ODE solver ode45. The solver employs an adaptive Dormand-Prince algorithm [27] based on the use of a fourth-and a fifth-order Runge-Kutta (RK) method pair. At every discrete independent variable step i → i + 1, the algorithm chooses the optimal Runge-Kutta coefficients to minimize the error of the fifth-order RK solution and also find the optimal variable step h i for efficient computation.
In particular, on each step, the difference between the fourth-and the fifth-order RK solution values is given by e i+1 = u [5] i+1 − u [4] i+1 , where each u k is a four-component vector providing a numerical approximation of the exact solution u = [y(x k ), y (x k ), y (x k ), y (x k )]. The one-step difference (118) is controlled by user-defined relative and absolute tolerances RelTol, AbsTol according to e i+1 ≤ max{RelTol · |u i |, AbsTol}. (119) If the ODE (97) is solved numerically for x ∈ [0, L] using N numerical steps, the conservative estimate of the global numerical error at x = L, for the small parameter value , is given by The difference between the numerical and the approximate solution at a numerical grid node x i is given by d( For a sample numerical-approximate solution computation, we use tolerance values For example, for = 0.1, this choice yields N = 381 steps in x, with variable step sizes h ranging from 0.00616 to 0.031948. Figure 1 shows numerical and approximate solution curves of y(x) as functions of x ∈ [0, L], L = 5, for = 0.02 and = 0.1. It is observed that, for = 0.02, the difference stays small for all x in the interval, while for a larger = 0.1, the numerical and approximate solutions begin to differ significantly after x 1. (We note that, for = 0, the approximate solution (117) becomes exact, and the difference (121) is negligible). To provide further details about the error and difference behaviour for the numerical and approximate solutions, Figure 2 shows the conservative estimate (120) of the total numerical error at x = L, the difference between the numerical and approximate solutions d(l; ) (121) at x = L as a function of , and also the typical behaviour of the difference (121) as a function of x for the specific small parameter value = 0.05.
The above analysis indicates that for sufficiently small values of the parameter , the approximate solution (117) of the perturbed Boussinesq Equation (97) indeed provides a precise approximation of the exact solution, with the error growing as the interval x ∈ [0, L] lengthens and/or the parameter is increased.

Discussion
In this paper, local symmetries of algebraic and ordinary differential equations involving a small parameter were considered in comparison to the symmetry structure of their unperturbed versions (small parameter equal to zero). Exact symmetries of the unperturbed equations, and exact and approximate symmetries (in the BGI framework [3][4][5]) of the perturbed models were investigated. The main goal of the paper was to address the question of stability of symmetries, when some given equation is perturbed by the addition of small O( ) terms.
It was observed by the original authors of the BGI method that while new and useful approximate symmetries can be sometimes found for perturbed models, some point symmetries of the unperturbed model may not appear in any form in the approximate point symmetry classification of a perturbed model, being thereby unstable. The aims of this paper were to find out the conditions under which a local symmetry becomes unstable, the form it can assume in the approximate point symmetry classification of a perturbed equation, and applications of approximate symmetries (in particular, higher-order ones) to compute approximate solutions of the given ODE with a small parameter.
It is straightforward to check that, for algebraic equations and first-order ODEs, every point symmetry of the unperturbed equation is stable (Section 3.1): a corresponding approximate point symmetry of the perturbed equation always exists; moreover, approximate point symmetry generators of perturbed algebraic equations and ODEs are more general than the exact symmetry generators of perturbed algebraic equations, and the approximate symmetry components arise as first-order Taylor terms in the expansion of exact symmetry components of the perturbed equation in the small parameter.
For second and higher-order ODEs and PDEs, the situation is more complex (Section 3.2): some original symmetries of the unperturbed model (38) can be unstable, in the sense of not being inherited as nontrivial approximate point symmetries of a perturbed ODE (39) (Example 3). At the same time, for some ODEs, all point symmetries of the unperturbed model might be stable (Example 4). This occurs because in the approximate point symmetry computation of a perturbed ODE, additional conditions on the O( 0 ) approximate symmetry components may or may not arise.
The situation is clarified in Section 4, where symmetries (point or local, exact, and approximate) are written in the evolutionary form. Theorem 2 is proven, showing that to every point or local symmetry of an exact ODE (38) of any order, there corresponds an approximate symmetry of the perturbed ODE (39), being possibly a higher-order symmetry of order at most n − 1. Two examples are considered in detail: a nonlinearly perturbed second-order ODE (12) (Section 4.3), and a fourth-order Boussinesq reduction ODE (64) (Section 4.4).
One of the most important applications of the approximate symmetry framework is the construction of closed-form approximate solutions to nonlinear ODE models with a small parameter. In Section 5, two approaches to obtain such solutions are developed. The first approach is based on approximate integrating factors using approximate point symmetries (Section 5.1).
Equations satisfied by approximate integrating factor components are derived (Theorem 3) and applied to obtain a four-parameter approximate solution family (101) of the fourth-order Boussinesq ODE (97). Another technique, approximate reduction of order under contact and higher-order symmetries, is presented in Section 5.3 and illustrated on two examples: an ODE (105) with a small parameter for which the exact general solution is known (Example 7), and again the fourth-order Boussinesq ODE (97) (Example 8). In the latter, the approximate solution is validated via a comparison to numerical solutions of the Boussinesq Equation (97).
In future work, it will be important to extend the understanding of relationships between symmetry structures of unperturbed and perturbed models in the cases of systems of ODEs, scalar PDEs, and systems of PDEs. Moreover, it is of high importance to investigate approaches to the computation of approximate symmetry properties of singularly perturbed models, including both ODE models (e.g., [28]) and PDE models, such as almost-inviscid Navier-Stokes fluids and shallow water equations [29] as well as weakly nonlinear models in elastodynamics [30].