2.1. Local symmetries of a PDE

A one-parameter (a) Lie group of transformations

(2.2) \begin{equation}\begin{array}{rcl} (x^*)^i & = & f^i (x,u;\,a) =x^i+a\xi_{0}^i(x,u)+ O\left(a^2\right) ,\quad i = 1,2, \ldots ,n, \\[1ex] u^* &= & g (x,u;\,a)= u +a \eta_0(x,u)+ O\left(a^2\right), \end{array}\end{equation}

in the space of the problem variables (x,u), with the corresponding infinitesimal generator

(2.3) \begin{equation}X^0 = \xi_{0}^i(x,u)\dfrac{\partial }{{\partial x^i }}+ \eta_0(x,u) \dfrac{\partial }{{\partial u }},\end{equation}

defines a point symmetry of the PDE (2.1) when the solution set of (2.1) is invariant under the action of the mapping (2.2); that is, each solution u(x) is mapped to a solution $u^*(x^*)$ of the same PDE. The local condition of this invariance is given by the determining equation

(2.4) \begin{equation}X^{0(k)} F_0[u]\big|_{F_0[u]=0}=0\end{equation}

in terms of the prolonged generator $X^{0(k)}$ . The determining equation splits into an overdetermined PDE system on the unknown symmetry components $\xi_{0}^i, \eta_0$ (see, e.g., [Reference Bluman, Cheviakov and Anco9, Reference Olver33]).

For a point symmetry (2.2) of the PDE (2.1), the evolutionary (characteristic) form providing the same mapping between solutions has the form

(2.5) \begin{equation}\begin{array}{rcl}(x^*)^i &= &x^i, \quad i = 1, 2, \ldots ,n,\\[1ex] u^* &= &u +a \zeta_0[u]+ O\left(a^2\right),\end{array}\end{equation}

with the evolutionary component $\zeta_0[u]={\eta_0}(x,u)-u_i \xi_{0}^i(x,u)$ , and the generator is given by

(2.6) \begin{equation}\hat{X}^0=\zeta_0[u]\dfrac{\partial}{\partial u}\,.\end{equation}

Local (point, contact and higher-order) transformations and the related point, contact and higher-order local symmetries of the PDE (2.1) generalize (2.5), (2.6) by allowing the evolutionary infinitesimal components $\zeta_0=\zeta_0[u]$ to be general differential functions of u, depending on first- and/or higher-order derivatives of u (see, e.g., [Reference Bluman, Cheviakov and Anco9, Reference Olver33] and references therein). The invariance condition (2.4) is replaced by

(2.7) \begin{equation}\hat{X}^{0(k)} F_0[u]\big|_{F_0[u]=0}=0.\end{equation}

The prolongation of $\hat{X}^0$ is defined by [Reference Olver33]

(2.8) \begin{equation} \hat{X}^{0(k)}=\zeta_0 \dfrac{\partial}{\partial u}+{\zeta_0}_{i}^{(1)} \dfrac{\partial }{{\partial u_i }}+\ldots + {\zeta_0}_{i_1 i_2 \ldots i_k }^{(k)} \dfrac{\partial}{{\partial u_{i_1 i_2 \ldots i_k } }}, \end{equation}

where the higher-order components are computed using

\begin{equation*} {\zeta_0}_{i}^{(1)} = D_i\zeta_0, \qquad {\zeta_0}_{i_1 i_2 \ldots i_p }^{(p)} =D_{i_p}{\zeta_0}_{i_1 i_2 \ldots i_{p-1} }^{(p-1)}, \end{equation*}

for $i,i_j=1,\ldots,n$ , $p=2,3,\ldots, k$ .

A point or local symmetry of a PDE is trivial if its components vanish on the solution set of the PDE (2.1). Trivial symmetries provide identity transformations (2.2), (2.5). Symmetry generators are commonly simplified modulo trivial symmetries.

The following elementary example will serve as a basis of further examples involving PDEs with a small parameter, and their exact and approximate symmetries.

Example 2.1. Consider a nonlinear wave-type equation [Reference Oron and Rosenau34]

(2.9) \begin{equation}u_{tt}=u_xu_{xx},\quad u=u(x,t).\end{equation}

The exact symmetry generator for the PDE (2.9) is given by

(2.10) \begin{equation}X^0 = \xi_{0}^1 (x,t,u)\dfrac{\partial }{{\partial x }} +\xi_{0}^2 (x,t,u)\dfrac{\partial }{{\partial t }}+ \eta_0 (x,t,u) \dfrac{\partial }{{\partial u }}.\end{equation}

The determining equations (2.4) yield the solution

(2.11) \begin{equation} \xi_{0}^1 = C_4+C_6 x,\quad \xi_{0}^2 = C_3+\left(C_6-\dfrac{C_5}{2}\right) t,\quad \eta_0 = C_1+C_2t+(C_5+C_6) u.\end{equation}

Consequently, the PDE (2.9) admits a six-dimensional Lie algebra of point symmetries spanned by

(2.12) \begin{equation}\begin{array}{rcl} X_1^0 = \dfrac{\partial}{\partial u},\quad X_2^0 = t\dfrac{\partial}{\partial u},\quad X_3^0 = \dfrac{\partial}{\partial t},\quad X_4^0 = \dfrac{\partial}{\partial x},\\[2ex] X_5^0 = u\dfrac{\partial}{\partial u}-\dfrac{t}{2}\dfrac{\partial}{\partial t},\quad X_6^0 =t\dfrac{\partial}{\partial t}+x\dfrac{\partial}{\partial x}+u\dfrac{\partial}{\partial u}.\end{array}\end{equation}

2.2. Exact and approximate symmetries of a PDE with a small parameter

A general first-order perturbation of a PDE (2.1) is a partial differential equation

(2.13) \begin{equation} F[u;\,\epsilon]=F_0[u]+\epsilon F_1[u]=o(\epsilon)\end{equation}

involving a small parameter $\epsilon$ . We assume that the perturbation $F_1[u]$ is regular, in the sense that the Kovalevskaya forms of the unperturbed PDE (2.1) and its perturbation (2.13) have the same leading derivatives.

2.2.1. Exact local symmetries

Exact point and local symmetry generators of (2.13) have the forms

(2.14) \begin{equation} Y=\alpha^{i}(x,u;\,\epsilon)\dfrac{\partial}{\partial x^i}+ \beta(x,u;\,\epsilon)\dfrac{\partial}{\partial u},\qquad \hat{Y}=\zeta[u;\,\epsilon]\dfrac{\partial}{\partial u}.\end{equation}

To find the exact symmetries of (2.13) holding for an arbitrary $\epsilon$ , one solves the determining equations (2.4) or (2.7).

Since (2.13) is a PDE family that includes an arbitrary element $\epsilon$ , the dimension of Lie algebra of point or local symmetries holding for a general $\epsilon$ cannot exceed that for some fixed $\epsilon$ , including $\epsilon=0$ . Therefore, the family (2.13) of perturbed PDEs will admit the same or smaller number of local symmetries than its unperturbed version (2.1).

Example 2.2. We compute exact point symmetries of a perturbed version of the PDE (2.9),

(2.15) \begin{equation} u_{tt}+\epsilon u u_t=u_xu_{xx}\end{equation}

holding for an arbitrary $\epsilon$ . The leading derivative $u_{tt}$ can be chosen for both (2.9) and (2.15). We obtain a Lie algebra of point symmetries spanned by

(2.16) \begin{equation} Y_1=X_3^0= \dfrac{\partial}{\partial t},\quad Y_2=X_4^0=\dfrac{\partial}{\partial x},\quad Y_3=\dfrac{4}{3}X_5^0-\dfrac{1}{3}X_6^0 = -t\dfrac{\partial}{\partial t}-\dfrac{x}{3}\dfrac{\partial}{\partial x} + u \dfrac{\partial}{\partial u},\end{equation}

a three-dimensional subalgebra of the six-dimensional Lie algebra of point symmetries (2.12) of the unperturbed wave equation (2.9).

2.2.2. Baikov–Gazizov–Ibragimov approximate point and local symmetries

Approximate symmetries can be useful for finding additional symmetry-like structures for the perturbed equation (2.13) that are not its exact symmetries but rather preserve (2.13) approximately, up to $o(\epsilon)$ [Reference Ibragimov27]. A one-parameter Lie group of Baikov–Gazizov–Ibragimov (BGI) approximate point transformations with the group parameter a, acting on the (x,u)-space, is given by

(2.17) \begin{equation}\begin{array}{rcll} {(x^*)}^{i} & = &f^i(x,u;\,a,\epsilon) = f_0^i(x,u;\,a)+\epsilon f_1^i(x,u;\,a) + o(\epsilon),& i=1,...,n, \\[1ex] {(u^*)} & = &g(x,u;\,a,\epsilon) = g_0(x,u;\,a)+\epsilon g_1(x,u;\,a) + o(\epsilon),&\end{array}\end{equation}

where, as in (2.2), $f_j^i$ , $g_j$ are sufficiently smooth functions. The generator of the approximate group (2.17), up to $o(\epsilon)$ , is

(2.18) \begin{equation}X = X^0+\epsilon X^1 = \left(\xi_0^i(x,u)+\epsilon \xi_1^i(x,u)\right)\dfrac{\partial}{\partial x^i}+\left(\eta_0(x,u)+\epsilon \eta_1(x,u)\right)\dfrac{\partial}{\partial u}.\end{equation}

Definition 2.1. The approximate group (2.17) of point transformations defines a BGI approximate point symmetry of the PDE (2.13) if it satisfies the approximate invariance condition of (2.13) under the action of (2.18):

(2.19) \begin{equation} (X^{0\,{(k)}}+\epsilon X^{1\,{(k)}})(F_0[u] +\epsilon F_1[u])\bigg |_{F_0[u] +\epsilon F_1[u]=o(\epsilon)}=o(\epsilon).\end{equation}

In (2.19), O(1) and $O(\epsilon)$ terms must vanish independently. It is easy to see that the O(1) term yields the determining equation (2.4) for the invariance of the unperturbed equation (2.1) under a point transformation $X^0$ (2.3). Hence, the following result holds.

Theorem 2.1. If the PDE (2.13) is approximately invariant under an approximate group of BGI point transformations with the generator (2.18) such that its O(1) part $X^0\not\equiv 0$ , then the infinitesimal operator $X^0$ (2.3) is a generator of an exact point symmetry group of the unperturbed PDE (2.1).

The converse of the above result does not always hold. Indeed, as it will be seen in examples below, if $X^0$ (2.3) generates an exact point symmetry group of the unperturbed PDE (2.1), there may be no corresponding BGI transformation (2.17) that approximately preserves the perturbed PDE (2.13). The following definition is important.

Definition 2.2. Suppose that the vector field $X^0$ given by (2.3) is a generator of an exact point symmetry group of the unperturbed PDE (2.1). If the perturbed PDE (2.13) admits an approximate generator X (2.18) with its O(1) part given by $X^0$ , then $X^0$ corresponds to a stable point symmetry of the unperturbed PDE (2.1) (in the BGI sense). Otherwise, it corresponds to an unstable point symmetry of (2.1).

Below in this paper, Definition 2.2 will be used not only for BGI approximate point symmetries, but more generally, for BGI and FS approximate point and local symmetries.

It is straightforward to show that $O(\epsilon)$ term of the determining equation (2.19) leads to the PDEs

(2.20) \begin{equation} X^{1\,{(k)}}F_0[u]\bigg|_{F_0[u]=0}=H[u],\end{equation}

where H is obtained from the coefficients of $\epsilon$ in the expression

(2.21) \begin{equation} -X^{0\,{(k)}}(F_0[u]+\epsilon F_1[u])\bigg|_{F_0[u]+\epsilon F_1[u]=o(\epsilon)}.\end{equation}

It follows that in order to calculate all approximate BGI symmetries (2.17), (2.18) of the perturbed PDE (2.13), one can take the following steps.

1. 1. Compute all generators $X^0$ (2.3) of exact point symmetry groups of the unperturbed PDE (2.1).

2. 2. Use each $X^0$ and the perturbation term $F_1[u]$ to compute H using (2.21).

3. 3. Solve (2.20) to find the components of $X^1$ .

An alternative procedure for the calculation of BGI symmetries involves writing down exact symmetry determining equations for $F[u;\,\epsilon]=0$ (cf. (2.13)), substituting $\zeta[u]=\zeta_0[u]+\epsilon \zeta_1[u]$ , and collecting O(1) and $O(\epsilon)$ coefficients of each split determining equation.

Remark 2.1. The first-order condition (2.20) may (or may not) contain additional conditions on the components $\xi_{0}^i, \eta_0$ of the unperturbed symmetry generator $X^0$ (2.3). This leads to the symmetry generated by $X^0$ being unstable (or respectively, stable). If all symmetries of the equations (2.1) are stable, the perturbed equations (2.13) are said to inherit the symmetry structure of the unperturbed equations [Reference Ibragimov27].

Remark 2.2. Similarly to exact local transformations with generators of the form (2.3), one can define more general local approximate BGI transformations with generators in evolutionary form given by

(2.22) \begin{equation}\hat{X} = \hat{X}^0+\epsilon \hat{X}^1 = \left(\zeta_0[u]+\epsilon \zeta_1[u]\right)\dfrac{\partial}{\partial u}\, .\end{equation}

Approximate local (including point, contact and higher-order) BGI symmetries of the perturbed PDE (2.13) can be found using the same procedure as described above for point approximate BGI symmetries. In particular, the analog of the first-order condition (2.20) takes the form

(2.23) \begin{equation} \left(\zeta_1 \dfrac{\partial}{\partial u}+{\zeta_1}_{i}^{(1)} \dfrac{\partial }{{\partial u_i }}+\ldots + {\zeta_1}_{i_1 i_2 \ldots i_k }^{(k)} \dfrac{\partial}{{\partial u_{i_1 i_2 \ldots i_k } }}\right)F_0[u]\bigg|_{F_0[u]=0}=H[u], \end{equation}

Theorem 2.1, the stability Definition 2.2 and Remark 2.1 concerning stability conditions of approximate symmetries directly carry over to the case of general local BGI symmetries.

2.2.3. Fushchich–Shtelen approximate point and local symmetries

Unlike the BGI approach where the symmetry generator is expanded in a power series in terms of the small parameter, the Fushchich–Shtelen method [Reference Fushchich and Shtelen21] applies the perturbation technique to the solution u(x) and the given PDE. In particular, the solution is written as

(2.24) \begin{equation} u(x) = v(x)+\epsilon w(x) + o(\epsilon).\end{equation}

Substituting (6.3) into the PDE (2.13) with a small parameter, expanding the result and setting to zero the O(1) and $O(\epsilon)$ terms independently, one obtains a coupled system of equations on v(x) and w(x) without the small parameter, given by

(2.25a) \begin{align} &G_1[v,w]\equiv F_0[v] = 0,\\[-5pt]\nonumber \end{align}

(2.25b) \begin{align} &G_2[v,w]\equiv {F_0}_v w+{F_0}_{v_i}w_i+{F_0}_{v_{ij}}w_{ij}+...+{F_0}_{v_{i_1i_2...i_k}}w_{i_1i_2...i_k}+F_1[v] = 0.\end{align}

It is clear that the first PDE (2.25a) is independent of w, and the second PDE (2.25b) is linear in w, with the linear operator being the Frèchet derivative of the $F_0[u]$ . We refer to equations (2.25) as the Fushchich–Shtelen system for the PDE (2.13). The PDE system (2.25) approximates the given PDE (2.13), in the sense that each exact solution pair (v(x), w(x)) of (2.25) yields an approximate solution (2.24) of the given PDE (2.13) up to the order $o(\epsilon)$ . Note that the Fushchich–Shtelen system (2.25) is in the extended Kovalevskaya form with respect to the similar leading derivatives of v and w as the leading derivatives of u in the original PDE (2.13).

Definition 2.3. A Lie group of point transformations with the group parameter a

(2.26) \begin{equation}\begin{array}{rcll} {(x^*)}^{i} & = &f^i(x,v,w;\,a) = x^i+a\, \lambda^i(x,v,w) + O\left(a^2\right),& i=1,...,n, \\[1ex] {(v^*)} & = &g(x,v,w;\,a) = v+ a\, \phi_1(x,v,w) + O\left(a^2\right),\\[1ex] {(w^*)} & = &h(x,v,w;\,a) = w+ a\, \phi_2(x,v,w) + O\left(a^2\right)\end{array}\end{equation}

with the generator

(2.27) \begin{equation} Z= \lambda^i(x,v,w)\dfrac{\partial}{\partial x^i}+ \phi_1(x,v,w) \dfrac{\partial}{\partial v}+ \phi_2(x,v,w) \dfrac{\partial}{\partial w}\end{equation}

defines a FS approximate point symmetry of the PDE (2.13) if it is an exact Lie point symmetry group of the Fushchich–Shtelen system (2.25).

In a similar manner, a generalized local (point or higher-order) transformation group in the evolutionary form

(2.28) \begin{equation}\begin{array}{ll} {(x^*)}^{i} &= x^i, \qquad i=1,...,n, \\[1ex] {(v^*)} &= v+ a \psi_1[v,w] + O\left(a^2\right),\\[1ex] {(w^*)} &= w+ a \psi_2[v,w] + O\left(a^2\right)\end{array}\end{equation}

with the generator

(2.29) \begin{equation} \hat{Z}=\psi_1[v,w]\dfrac{\partial}{\partial v}+\psi_2[v,w] \dfrac{\partial}{\partial w}\end{equation}

defines a local (point or higher-order) FS approximate symmetry of the PDE (2.13) if it is a local symmetry of the Fushchich–Shtelen system (2.25).

It is important to know whether the FS approximate symmetry structure of a PDE (2.13) with a small parameter is in some sense inherited from exact local symmetries of the unperturbed PDE (2.1). Similarly to the BGI case, one can define stable and unstable symmetries in the Fushchich–Shtelen framework.

Definition 2.4. Suppose $\hat{X}_0=\zeta_0[u]\,{\partial}/{\partial u}$ (2.22) is a generator of an exact local symmetry group of the unperturbed PDE (2.1). If the perturbed PDE (2.13) admits an approximate FS symmetry with generator (2.29) where the v-component $\psi_1[v,w]\equiv\zeta_0[v]$ , then $\hat{X}_0$ corresponds to a stable point symmetry of the unperturbed PDE (2.1) (in the FS sense). Otherwise, it corresponds to an unstable point symmetry of (2.1).

Similarly to the case for BGI approximate symmetries, a FS approximate symmetry of a PDE (2.13) given by (2.29) may be unstable because the second symmetry determining equation for the Fushchich–Shtelen system (2.25)

\begin{equation*} \hat{Z}^{(k)}G_2[v,w]\bigg|_{G_1[v,w]=G_2[v,w]=0} = 0\end{equation*}

could contain additional conditions on the v-component $\psi_1$ .

Remark 2.3. An important feature of the Fushchich–Shtelen approximate symmetry framework is the possibility of existence of approximate FS symmetries where the O(1) component $\psi_1[v,w]$ of the generator (2.29) depends on $O(\epsilon)$ solution component w. An example is provided by the linear PDE $u_{tt}+\epsilon u_t = u_{xx}$ which admits a point FS symmetry (2.29) with

(2.30) \begin{equation}\begin{array}{ll}\psi_1[v,w] = tv+2w +2xt v_x +(x^2+t^2)v_t,\\[1ex]\psi_2[v,w] = \dfrac{1}{2}x^2 v -tw +2xt w_x +(x^2+t^2)w_t\,.\end{array}\end{equation}

Such FS symmetries do not originate from stable local symmetries of the unperturbed PDE $u_{tt} = u_{xx}$ and cannot arise in the BGI framework.

3.1. Trivial approximate symmetries

3.1.1. Trivial BGI approximate symmetries

Consider a local BGI approximate transformation with the evolutionary generator (2.22):

\begin{equation*}\hat{X} = \hat{X}^0+\epsilon \hat{X}^1 = \left(\zeta_0[u]+\epsilon \zeta_1[u]\right)\dfrac{\partial}{\partial u}\, .\end{equation*}

The determining equations (2.19) for the generator (2.22) to define a local symmetry of the PDE (2.13) with a small parameter split into the O(1) part (2.4) and $O(\epsilon)$ part (2.20) with H defined by (2.21). Suppose that the O(1) part of the generator vanishes: $\hat{X}_0=0$ . In that case, the O(1) part (2.4) of the approximate symmetry determining equations is satisfied identically, and (2.21) yields $H=0$ . Consequently, the $O(\epsilon)$ part (2.20) of the determining equations (2.19) becomes

(3.1) \begin{equation} \hat{X}^{1\,{(k)}}F_0[u]\bigg|_{F_0[u]=0}=0,\end{equation}

which means that such $\hat{X}^{1}$ must be a local symmetry generator of the unperturbed equation (2.1). The opposite is also true: if $\hat{X}^{0}$ is a local symmetry generator of the unperturbed equation (2.1), then

(3.2) \begin{equation}\hat{X} = \epsilon \hat{X}^0\end{equation}

is a BGI approximate symmetry generator of the perturbed PDE (2.13). In the light of the above, we call a BGI approximate symmetry that has a generator with vanishing O(1) part

(3.3) \begin{equation}\hat{X} = \epsilon \hat{X}^1 = \epsilon \zeta_1[u]\dfrac{\partial}{\partial u}\end{equation}

a trivial BGI approximate symmetry. This triviality relates not to the trivial action of such symmetries but rather to the fact that every local symmetry $\hat{X}^0$ of the unperturbed equation (2.1) is guaranteed to yield a BGI approximate symmetry of the perturbed equation (2.13) having the form (3.2). The local action of a trivial BGI approximate symmetry in the evolutionary form defined by (3.3) is given by

(3.4) \begin{equation}\begin{array}{rcl}(x^*)^i &= &x^i, \quad i = 1, 2, \ldots ,n,\\[1ex] u^* &= &u + a\epsilon\zeta_1[u]+ O\left(a^2\right),\end{array}\end{equation}

with the first Taylor term of the transformation having the order $\sim a\epsilon = o(a,\epsilon)$ .

3.1.2. Trivial FS approximate symmetries

In a parallel fashion, one can define a trivial FS approximate symmetry of the perturbed PDE (2.13) as one for which the local generator (2.29) has a special form with the vanishing transformation component of the O(1) part of the solution $\psi_1=0$ , and $\psi_2[v,w]=\psi_2[v]$ :

(3.5) \begin{equation} \hat{Z}=0+\psi_2[v] \dfrac{\partial}{\partial w}\,.\end{equation}

(The respective trivial point FS symmetry generator (2.27) has $\lambda^i=\lambda^i(x,v)$ , $\phi_1=0$ , $\phi_2=\phi_2(x,v)$ .) For FS local symmetries with the generator of the form (3.5), it is straightforward to show that $\psi_2[v]$ is an evolutionary component of the local symmetry of the unperturbed equation (2.1) generated by

(3.6) \begin{equation} \hat{X}^0= \psi_2[u] \dfrac{\partial}{\partial u}.\end{equation}

Indeed, the action of (3.5) on the first PDE (2.25a) of Fushchich–Shtelen system is trivial, and the action on the linear equation (2.25b) is equivalent to the local symmetry determining equation (2.7) of the unperturbed PDE (2.1). The converse also holds: every local symmetry generator (3.6) of the unperturbed equation yields a trivial FS approximate symmetry generator (3.5).

3.2. Types of approximate symmetries arising in classifications

In the computation of BGI approximate symmetries of a PDE (2.13) with a small parameter, the following three types of symmetries can arise.

1. 1. BGI approximate symmetries with generators (2.22) having $\hat{X}^0 \ne 0$ , $\hat{X}^1=0$ correspond to exact local symmetries of the perturbed equation (2.13) (see Section 2.2.1).

2. 2. BGI approximate symmetries with generators having $\hat{X}^0 = 0$ , $\hat{X}^1 \ne 0$ correspond to trivial BGI approximate symmetries (Section 3.1.1).

3. 3. Genuine BGI approximate symmetries have generators with both $\hat{X}^0\ne 0$ and $\hat{X}^1 \ne 0$ .

For FS approximate symmetries, the following three types can arise.

1. 1. Symmetries with the same action on O(1) solution part v and $O(\epsilon)$ solution part w correspond to exact local symmetries of the perturbed equation (2.13). For example, an exact scaling symmetry with the generator $u\,\partial/\partial u$ admitted by the perturbed equation (2.13) is equivalent to a FS scaling symmetry with the generator $v\,\partial/\partial v + w\,\partial/\partial w$ .

2. 2. Trivial FS approximate symmetries, as defined in Section 3.1.2.

3. 3. Genuine FS approximate symmetries.

Genuine BGI and FS approximate symmetries are the main focus of the approximate symmetry study.

3.3. A computational example: Exact and approximate point symmetry classification for a second-order nonlinear PDE with a small parameter

In this section, we compare exact point symmetries (2.12) of the (1+1)-dimensional wave-type equation (2.9)

\begin{equation*}u_{tt}=u_xu_{xx}\end{equation*}

with BGI and FS approximate point symmetry classifications for the family of perturbed equations

(3.7) \begin{equation} u_{tt}+\epsilon F_1(u, u_t)=u_xu_{xx},\end{equation}

where $F_1(u, u_t)$ is an arbitrary function. The computation is limited to point symmetries for simplicity of presentation. The cases for the arbitrary function $F_1(u, u_t)$ are given modulo equivalence transformations

(3.8) \begin{equation} t=C_1^{-1}\,\tilde{t},\quad x=C_1^{-2/3}\,\tilde{x},\quad u=\tilde{u}+\dfrac{1}{2}\,\epsilon \,\left(1-C_1^{-2}\right)\, C_2 \,\tilde{t}^2,\quad F_1 = C_1^2\tilde{F}_1-\left(C_1^2-1\right)C_2,\end{equation}

where $C_1$ and $C_2$ are arbitrary constants.

Table 1. Stability of point symmetries (2.12) of the wave-type equation in terms of BGI and FS approximate point symmetries of the perturbed PDE (3.7), depending on the form of the arbitrary function $F_1$

The PDE (2.9) admits six exact point symmetries given by (2.12). The BGI approximate point symmetries are computed and classified following the procedure described in Section 2.2.2, and the FS approximate local symmetries are obtained following Section 2.2.3. In particular, the Fushchich–Shtelen system (2.25) for (3.7) is given by

(3.9) \begin{equation} v_{tt}-v_xv_{xx} = 0,\qquad w_{tt}+F_1(v,v_t)-v_xw_{xx}-w_xv_{xx}=0, \end{equation}

where $u(x,t)=v(x,t)+\epsilon w(x,t)$ . For the Fushchich–Shtelen system (3.9), point equivalence transformations are given by (3.8) for x, t and $F_1$ , with dependent variable transformations

(3.10) \begin{equation}v=\tilde{v}, \qquad w=\tilde{w}+\dfrac{1}{2}\,\left(1-C_1^{-2}\right)\, C_2 \,\tilde{t}^2.\end{equation}

The resulting classification is presented in Table 1. In the table, $Q_i$ denote arbitrary functions of their arguments and $a_j$ arbitrary constants. The table is organized as follows. The first column lists evolutionary forms of the six point symmetry generators $X^0_k$ (2.12), $k=1,\ldots, 6$ , of the unperturbed PDE (2.9). The second column lists the forms of the arbitrary function $F_1(u, u_t)$ for which the corresponding $X^0_k$ is stable in the BGI sense, and the corresponding BGI approximate point symmetry of the perturbed wave equation (3.7). The third column contains the same information for the FS approximate point symmetries of the perturbed equation (3.7). Trivial approximate symmetries are not listed (see Section 3.1).

In Table 1, the generators of BGI and FS symmetries were sought in the evolutionary forms, depending on at most first derivatives of the fields. As a result, approximate BGI and FS point symmetries were obtained, as well as some local first-order approximate FS symmetries. We also note that for the point symmetry defined by $\hat{X}^0_5$ to be stable in the FS sense, the arbitrary function $F_1=F_1(v, v_t)$ must satisfy three coupled nonlinear PDEs. Because of the complexity of the latter, a general form of $F_1$ could not be obtained explicitly; however, it was verified that the form of $F_1$ as listed satisfies those equations. The same takes place for the point symmetry defined by $\hat{X}^0_6$ .

Table 1 illustrates differences between BGI and FS frameworks as there are stable symmetries in one framework and unstable in the other framework. Some specific examples are listed below. The following observations can be made.

• • When $F_1(u,u_t)=u_t$ or $F_1={\rm const}$ , all point symmetries (2.12) of the unperturbed PDE (2.9) are stable as BGI and FS symmetries.

• • For all forms $F_1(u,u_t)$ , the PDE (3.7) and the FS system (3.9) are invariant under t- and x-translations. Consequently, the exact symmetries $\hat{X}^0_3$ and $\hat{X}^0_4$ reappear as BGI and FS approximate symmetries without change.

• • As an example, as seen in the first row corresponding to $\hat{X}^0_1$ , for $F_1=uu_t$ , this symmetry is stable as a point BGI approximate symmetry and as a local (but not point) first-order FS approximate symmetry.

• • It can be separately proven that for every form of $F_1$ where a symmetry $\hat{X}^0_i$ is stable as a BGI point approximate symmetry, it is also stable as a local (point or first-order) FS approximate symmetry.

• • The converse is not true; for example, the scaling symmetry $\hat{X}^0_6$ is stable as a local FS approximate symmetry for $F_1= v^{a_3}Q_6(v_t)+a_2v_t + a_1$ , but is only stable as a BGI approximate symmetry when $a_3=-1$ .

We also note that genuine BGI approximate symmetries are given by $\hat{X}_1$ , $\hat{X}_2$ , $\hat{X}_5$ and $\hat{X}_6$ ; similarly, genuine FS approximate symmetries correspond to the local generators $\hat{Z}_1$ , $\hat{Z}_2$ , $\hat{Z}_5$ and $\hat{Z}_6$ .

3.4. Instability of local symmetries of unperturbed PDEs in terms of higher-order approximate symmetries: An example

For an ordinary differential equation (ODE), all local symmetries are stable in the BGI sense: each local symmetry of a given ODE corresponds to an approximate BGI local, often higher-order, symmetry of its perturbed version [Reference Tarayrah and Cheviakov40]. For a PDE, in general, this is not the case. For the BGI framework, differential functions ${\zeta_1}_{i_1 i_2 \ldots i_p }^{(p)} $ in the determining equation (2.23) contain derivatives of u of orders higher than those in the differential function $\zeta^1$ . It follows that the left-hand side of equation (2.23) splits into a system of linear PDEs in $\zeta_1$ . On the other side, the function H may contain derivatives of u with respect to other variables different than those in the left-hand side of equation (2.23). This can lead to some constraints on the unperturbed symmetry component $\zeta_0$ ; in that case, an exact local symmetry of the unperturbed PDE (2.1) may not correspond to a local approximate BGI symmetry of the perturbed PDE (2.13). A similar argument holds for FS approximate symmetries. The main reason, as it can be seen in the example below, is the existence of multiple kinds of derivatives in PDEs and thus more restrictive conditions that arise for $\zeta_0$ when the determining equations are being split with respect to higher-order derivatives.

As an illustration, consider the PDE (3.7) with $F_1(u,u_t)=uu_t$ :

(3.11) \begin{equation} u_{tt}+\epsilon uu_t=u_xu_{xx}\end{equation}

and the related Fushchich–Shtelen system (3.9)

(3.12) \begin{equation}v_{tt}-v_xv_{xx} = 0,\qquad w_{tt}+vv_t-v_xw_{xx}-w_xv_{xx}=0.\end{equation}

From Table 1, one can see that $\hat{X}^0_2=t \,\partial/\partial u$ is unstable as a point symmetry in both BGI and FS frameworks; that is, the point symmetry $\hat{X}^0_2$ admitted by the PDE (3.11) with $\epsilon=0$ corresponds to no approximate point or first-order symmetry arising from BGI or FS approaches. First, we examine whether or not it is possible to construct a local, possibly higher-order, BGI approximate symmetry of (3.11) that would correspond to $\hat{X}^0_2$ . The generator of such a symmetry would have the form

(3.13) \begin{equation}\hat{X}_2 =C_2\hat{X}_2^0+\epsilon \hat{X}_2^1=\left( C_2 t +\epsilon \zeta_1[u]\right)\dfrac{\partial}{\partial u},\end{equation}

where $C_2={\rm const}\ne 0$ . The determining equation (2.23) for BGI local symmetries reads

(3.14) \begin{equation} \left(D_t^2\zeta_1-u_xD_x^2\zeta_1-u_{xx}D_x\zeta_1\right)\bigg|_{u_{tt}=u_xu_{xx}}={H},\end{equation}

where one readily finds

(3.15) \begin{equation}H=C_2(tu_t+u).\end{equation}

One can show by a direct computation that whatever the dependence of $\zeta_1$ on partial derivatives of u is chosen to be, higher-order derivatives of u that arise in (3.14) lead to constraints on $C_2$ that result in $C_2=0$ , which means that no nontrivial BGI point symmetry (3.13) corresponding to $\hat{X}^0_2$ exists.

Second, we seek a local, possibly higher-order, approximate FS symmetry of the PDE (3.11) corresponding to $\hat{X}^0_2$ . Such a symmetry would arise as an exact local symmetry of Fushchich–Shtelen system (3.12). The corresponding evolutionary generator (2.29) has the form

(3.16) \begin{equation} Z= \psi_1[v,w] \dfrac{\partial}{\partial v}+\psi_2[v,w] \dfrac{\partial}{\partial w}.\end{equation}

As noted in Section 2.2.3, the determining equation for the first equation of the system (3.12) is satisfied when $\psi_1=C_2t$ as in $\hat{X}^0_2$ . Now the determining equation for the second PDE of (3.12) leads to

(3.17) \begin{equation} \left(D_t^2\psi_2-v_xD_x^2\psi_2-v_{xx}D_x\psi_2\right)\bigg|_{v_{tt}=v_xv_{xx},w_{tt}=-vv_t+v_xw_{xx}+w_xv_{xx} }=C_2(tv_t+v).\end{equation}

It can be shown that for any dependence $\psi_2[v,w]$ , constraints on $C_2$ exist, leading to $C_2=0$ . Consequently, there is no higher-order FS symmetry corresponding to the unstable point symmetry $\hat{X}^0_2$ admitted by the wave equation (3.11) with $\epsilon=0$ .

3.5. A relation between BGI and FS approximate symmetries

The computational example of Section 3.3 above illustrated the fact that BGI and FS frameworks can yield rather different approximate point symmetry classifications for the same PDE with a small parameter. However, in certain situations, the two approaches can lead to related results. We now show that for a specific class of (1+1)-dimensional PDEs, a stable BGI approximate point symmetry always correspond to a stable FS approximate local symmetry.

Consider the following class of PDEs on u(x,t), written in the Kovalevskaya form with respect to an independent variable t:

(3.18) \begin{equation} \dfrac{\partial ^n u}{\partial t^n} =F_0[u],\qquad F_0[u]\equiv F_0(x,t,u,\partial u, \partial^2 u,...,\partial^k u),\end{equation}

and its perturbed version with a small parameter $\epsilon$ :

(3.19) \begin{equation} \dfrac{\partial ^n u}{\partial t^n}=F_0[u]+\epsilon F_1[u],\qquad F_1[u]\equiv F_1(x,t,u,\partial u, \partial^2 u,...,\partial^\ell u).\end{equation}

A local BGI approximate symmetry of a PDE (3.19) has the form (2.22)

(3.20) \begin{equation}\hat{X} = \hat{X}^0+\epsilon \hat{X}^1 = \left(\zeta_0[u]+\epsilon \zeta_1[u]\right)\dfrac{\partial}{\partial u}\, .\end{equation}

As per Theorem 2.1 and Remark 2.2, the O(1) term in (3.20) corresponds to a local symmetry

(3.21) \begin{equation}\hat{X}^0 = \zeta_0[u]\,\dfrac{\partial}{\partial u}\end{equation}

of the unperturbed equation (3.18).

In order to compute FS approximate symmetries of a PDE (3.19), we substitute $u(x,t)=v(x,t)+\epsilon w(x,t)+o(\epsilon)$ into (3.19) and split the orders of $\epsilon$ to get the Fushchich–Shtelen system

(3.22) \begin{align} \dfrac{\partial ^n v}{\partial t^n} & =F_0[v],\nonumber \\[4pt]\dfrac{\partial ^n w}{\partial t^n} & ={F_0}_vw+{F_0}_{v_i}w_i+{F_0}_{v_{ij}}w_{ij}+...+{F_0}_{v_{i_1i_2...i_k}}w_{i_1i_2...i_k}+F_1[v]. \end{align}

The evolutionary generator of a local FS approximate symmetry has the form (2.29). The determining equations (2.4) for exact local symmetries of (3.22) are

(3.23a) \begin{align} &\hat{Z}^{(n)}\left(\dfrac{\partial ^n v}{\partial t^n}- F_0\right)=0,\\[-10pt]\nonumber \end{align}

(3.23b) \begin{align} &\hat{Z}^{(n)}\left(\dfrac{\partial ^n w}{\partial t^n}-{F_0}_vw-{F_0}_{v_i}w_i-{F_0}_{v_{ij}}w_{ij}-...-{F_0}_{v_{i_1i_2...i_k}}w_{i_1i_2...i_k}-F_1\right)=0,\end{align}

holding on solutions of (3.22).

Theorem 3.1. If (3.20) is a BGI approximate local symmetry generator of a PDE (3.19) having the specific form

(3.24) \begin{equation}\hat{X}=\left(\zeta_0(x,t)+\epsilon \zeta_1(x,t,u,u_x,u_t)\right)\dfrac{\partial}{\partial u}\end{equation}

and additionally, $F_0[v]$ in (3.19) satisfies the following system of equations

(3.25) \begin{align} & \zeta^0{F_0}_{uu}+\zeta_i^{0^{(1)}}{F_0}_{uu_i}+\zeta_{i_1i_2}^{0^{(2)}}{F_0}_{uu_{i_1i_2}}+...+\zeta_{i_1i_2...i_k}^{0^{(k)}}{F_0}_{uu_{i_1i_2...i_k}}=0,\nonumber\\[3pt] & \zeta^0{F_0}_{uu_i}+\zeta_i^{0^{(1)}}{F_0}_{u_iu_i}+\zeta_{i_1i_2}^{0^{(2)}}{F_0}_{u_iu_{i_1i_2}}+...+\zeta_{i_1i_2...i_k}^{0^{(k)}}{F_0}_{u_iu_{i_1i_2...i_k}}=0,\nonumber\\[3pt]& \vdots \\[3pt]& \zeta^0{F_0}_{uu_{i_1i_2...i_k}}+\zeta_i^{0^{(1)}}{F_0}_{u_iu_{i_1i_2...i_k}}+\zeta_{i_1i_2}^{0^{(2)}}{F_0}_{u_{i_1i_2}u_{i_1i_2...i_k}}\nonumber\\[3pt] & \qquad \qquad +...+\zeta_{i_1i_2...i_k}^{0^{(k)}}{F_0}_{u_{i_1i_2...i_k}u_{i_1i_2...i_k}}=0,\nonumber\end{align}

then

(3.26) \begin{equation}\hat{Z}=\zeta^0(x,t)\dfrac{\partial}{\partial v}+\zeta^1(x,t,v,v_x,v_t) \dfrac{\partial}{\partial w}\end{equation}

is a FS approximate local symmetry of the perturbed PDE (3.19) to the point symmetry generator $\hat{X}^{0}=\zeta^0\,\partial/\partial v$ of the unperturbed PDE (3.18).

Proof. We need to show that under the stated conditions, the determining equations (2.20) for BGI approximate symmetries of (3.19) are equivalent to the determining equations for FS approximate symmetries. Since the first PDE of the Fushchich–Shtelen system (3.22) is the same as the unperturbed equation (3.18), the first FS determining equation (3.23a) is satisfied for any $\zeta^0$ and $\zeta^1$ as long as $\zeta^0$ is an exact point symmetry component of (3.18).

The second FS determining equation (3.23b) with $\psi_1=\zeta^0$ , $\psi_2=\zeta^1$ can be rewritten as

(3.27) \begin{equation} \left(\zeta_t^{1^{(n)}}-\zeta^1{F_0}_{v}-\zeta_i^{1^{(1)}}{F_0}_{v_i}-\zeta_{i_1i_2}^{1^{(2)}}{F_0}_{v_{i_1i_2}}-...-\zeta_{i_1i_2...i_k}^{1^{(k)}}{F_0}_{v_{i_1i_2...i_k}}\right)\bigg|_{{\partial ^n v}/{\partial t^n}=F_0}=G,\end{equation}

where

\begin{align*} G&=w\left(\zeta^0{F_0}_{vv}+\zeta_i^{0^{(1)}}{F_0}_{vv_i}+\zeta_{i_1i_2}^{0^{(2)}}{F_0}_{vv_{i_1i_2}}+...+\zeta_{i_1i_2...i_k}^{0^{(k)}}{F_0}_{vv_{i_1i_2...i_k}}\right)\nonumber\\[4pt] &\quad +w_i\left(\zeta^0{F_0}_{vv_i}+\zeta_i^{0^{(1)}}{F_0}_{v_iv_i}+\zeta_{i_1i_2}^{0^{(2)}}{F_0}_{v_iv_{i_1i_2}}+...+\zeta_{i_1i_2...i_k}^{0^{(k)}}{F_0}_{v_iv_{i_1i_2...i_k}}\right) +...\nonumber\\[4pt] &\quad +w_{i_1i_2...i_k}\left(\zeta^0{F_0}_{vv_{i_1i_2...i_k}}+\zeta_i^{0^{(1)}}{F_0}_{v_iv_{i_1i_2...i_k}}+\zeta_{i_1i_2}^{0^{(2)}}{F_0}_{v_{i_1i_2}v_{i_1i_2...i_k}}+... \right.\\[4pt]& \quad \left. +\,\zeta_{i_1i_2...i_k}^{0^{(k)}}{F_0}_{v_{i_1i_2...i_k}v_{i_1i_2...i_k}}\right) + \zeta^0{F_1}_v+\zeta_i^{0^{(1)}}{F_1}_{v_i}+...+\zeta_{i_1i_2...i_\ell}^{0^{(\ell)}}{F_1}_{v_{i_1i_2...i_\ell}}.\nonumber\end{align*}

As $\zeta^0$ and $F_0[v]$ satisfy (3.25), G reduces to

(3.28) \begin{equation} G= \zeta^0{F_1}_v+\zeta_i^{0^{(1)}}{F_1}_{v_i}+...+\zeta_{i_1i_2...i_\ell}^{0^{(\ell)}}{F_1}_{v_{i_1i_2...i_\ell}}.\end{equation}

Now, we proceed to check the determining equation (2.20) of BGI approximate symmetries for (3.19). The left-hand side of (2.20) simplifies to

\begin{equation*} \left(\zeta_t^{1^{(n)}}-\zeta^1{F_0}_{u}-\zeta_i^{1^{(1)}}{F_0}_{u_i}-\zeta_{i_1i_2}^{1^{(2)}}{F_0}_{u_{i_1i_2}}-...-\zeta_{i_1i_2...i_k}^{1^{(k)}}{F_0}_{u_{i_1i_2...i_k}}\right)\bigg|_{{\partial ^n u}/{\partial t^n}=F_0}\end{equation*}

which is equivalent to the left-hand side of (3.27). Now, the right-hand side of (2.20), the function H, is the coefficient of $\epsilon$ in

(3.29) \begin{equation} -\hat{X}^{0^{(n)}}\left(\dfrac{\partial ^n u}{\partial t^n}-F_0-\epsilon F_1\right)\bigg|_{{\partial ^n u}/{\partial t^n}=F_0+\epsilon F_1}.\end{equation}

Since $\zeta^0=\zeta^0(x,t)$ , none of the terms in (3.29) contains $\partial^n u/\partial t^n$ . Hence, the coefficient of $\epsilon$ in (3.29) is

(3.30) \begin{equation} H=\hat{X}^{0^{(n)}} F_1= \zeta^0{F_1}_u+\zeta_i^{0^{(1)}}{F_1}_{u_i}+...+\zeta_{i_1i_2...i_\ell}^{0^{(\ell)}}{F_0}_{u_{i_1i_2...i_\ell}}.\end{equation}

The latter is equivalent to G (3.28). It follows that the determining equation (3.5) of FS symmetries for the system (3.22) and the determining equation (2.20) of BGI approximate symmetries for the PDE (3.19) are equivalent. Hence, $\hat{Z}$ (3.26) is a FS approximate local symmetry of the system (3.22).

The above theorem states that when a point symmetry of an unperturbed PDE yields a BGI approximate point symmetry but not a FS approximate point symmetry of the perturbed PDE, under the conditions of the theorem, there exists a corresponding higher-order FS approximate symmetry of the perturbed PDE instead. This is illustrated by the following example.

Example 3.1 Consider again the PDE (3.11) $u_{tt}+\epsilon uu_t=u_xu_{xx}$ . Using Table 1, we observe that $X^0_1=\partial/\partial u$ is unstable as a FS approximate point symmetry but it is a stable point symmetry in the sense of BGI; the corresponding BGI generator is given by

\begin{equation*}\hat{X}_1=\left(1-\dfrac{\epsilon}{10}t(tu_t+4u)\right)\dfrac{\partial}{\partial u},\end{equation*}

with H (3.30) given by (3.15) with $C_1=1$ . Consider now the Fushchich–Shtelen system (3.12) for the PDE (3.11). Using determining equation (3.5) for exact local symmetries of (3.12), one can find that

\begin{equation*}\hat{Z}_1=\dfrac{\partial}{\partial v}-\dfrac{1}{10}t(tv_t+4v)\dfrac{\partial}{\partial w}\end{equation*}

is a higher-order FS approximate symmetry generator of the PDE (3.11).

Example 3.2. The perturbed wave equation

(3.31) \begin{equation} u_{tt}+\epsilon uu_t=e^{u}u_{xx},\end{equation}

admits an approximate point symmetry with the evolutionary form

(3.32) \begin{equation} \hat{X}=\left(2-xu_x-\epsilon t\left(tu_t+4\right)\right)\dfrac{\partial}{\partial u}\end{equation}

corresponding to the stable point symmetry $\hat{X}^0=\zeta^0\,{\partial}/{\partial u}=(2-xu_x)\,{\partial}/{\partial u}$ . Here, $\zeta^0=2-xu_x$ does not satisfy the conditions of Theorem 3.1 since it involves $u_x$ . It turns out that $\hat{X}^0$ is unstable as a FS approximate point symmetry of (3.31). Indeed, it is easy to check that

(3.33) \begin{equation}\hat{Z}=\left(2-xv_x\right)\dfrac{\partial}{\partial v}-t\left(tv_t+4\right)\dfrac{\partial}{\partial w}\end{equation}

is not a local symmetry of the Fushchich–Shtelen system of the PDE (3.31) given by

(3.34) \begin{equation} v_{tt}-v_xv_{xx} =0, \qquad w_{tt}+vv_t-e^vw_{xx}-e^vwv_{xx} =0.\end{equation}

However, a direct computation shows that a slightly modified version of (3.33) given by

\begin{equation*}\hat{Z}=\left(2-xv_x\right)\dfrac{\partial}{\partial v}-\left(t(tv_t+4)+xw_x\right)\dfrac{\partial}{\partial w}\end{equation*}

yields the corresponding first-order FS approximate symmetry.

4.1. Exact point symmetries of the wave equation family (4.6)

The exact symmetry generator for the PDE (4.6) has the form

(4.7) \begin{equation} Y=\xi^1(x,t,u;\,\epsilon)\dfrac{\partial}{\partial t}+\xi^2(x,t,u;\,\epsilon)\dfrac{\partial}{\partial x}+\eta(x,t,u;\,\epsilon)\dfrac{\partial}{\partial u}.\end{equation}

The following cases arise, holding for an arbitrary $\epsilon$ and non-constant Q.

1. 1. In the general case of arbitrary $Q(u_x)$ , one has the five-dimensional Lie group of point symmetries generated by

   (4.8) \begin{equation}Y_1=\dfrac{\partial}{\partial t},\quad Y_2=\dfrac{\partial}{\partial x},\quad Y_3=\dfrac{\partial}{\partial u},\quad Y_4 = t\dfrac{\partial}{\partial u},\quad Y_5=t\dfrac{\partial}{\partial t}+x\dfrac{\partial}{\partial x}+u\dfrac{\partial}{\partial u},\end{equation}
   describing respectively translations in t, x, u, the Galilei transformation in the direction of the displacement u, and a homogeneous space-time scaling.

2. 2. In the case when $Q(u_x)=u_x$ , the Lie algebra (4.8) is extended by a point symmetry generator

   (4.9) \begin{equation} Y_6=x\dfrac{\partial}{\partial u}-\epsilon\left(x\dfrac{\partial}{\partial x}+\dfrac{3t}{2}\dfrac{\partial}{\partial t}\right). \end{equation}

4.2. BGI approximate point symmetry classification of the wave equation family (4.6)

The BGI approximate point symmetry generator for the PDE (4.6) has the form

(4.10) \begin{equation}X= X^0+\epsilon X^1=X^0 + \epsilon \left(\xi_1^1(x,t,u)\dfrac{\partial}{\partial t}+\xi_1^2(x,t,u)\dfrac{\partial}{\partial x}+\eta_1(x,t,u)\dfrac{\partial}{\partial u}\right),\end{equation}

where, according to Theorem 2.1, the freedom in $X^0$ does not exceed that in $X^0_\infty$ (4.3). From the determining equations (2.19) for BGI point symmetries, the following cases arise.

1. 1. $Q(u_x)$ arbitrary: the O(1) and $O(\epsilon)$ components of the generator (4.10) are given by

   (4.11a) \begin{align} X^0=&(C_1t+C_2)\dfrac{\partial}{\partial t}+(C_1x+C_3)\dfrac{\partial}{\partial x}+(C_1u+C_4t+C_5)\dfrac{\partial}{\partial u},\\[-10pt]\nonumber \end{align}
   (4.11b) \begin{align} X^1=&(\lambda_1+\lambda_2)\dfrac{\partial}{\partial t}+ \left(\lambda_1-\lambda_2\right)\dfrac{\partial}{\partial x}+(C_6u+\lambda_3+\lambda_4)\dfrac{\partial}{\partial u}, \end{align}

   where $C_i={\rm const}$ , $\lambda_1$ and $\lambda_3$ are arbitrary functions of $x+t$ , and $\lambda_2$ and $\lambda_4$ are arbitrary functions of $x-t$ . Consequently, the nonlinear wave equation (4.6) for an arbitrary $Q(u_x)$ admits the approximate point symmetries

   (4.12) \begin{align} X_1 & =t\dfrac{\partial}{\partial t}+x\dfrac{\partial}{\partial x}+u\dfrac{\partial}{\partial u},\quad X_2=\dfrac{\partial}{\partial t},\quad X_3=\dfrac{\partial}{\partial x},\quad X_4=t\dfrac{\partial}{\partial u},\quad X_5 = \dfrac{\partial}{\partial u},\nonumber\\[5pt]X_\infty & =\epsilon\left[\left(\lambda_1+\lambda_2\right)\dfrac{\partial}{\partial t}+\left(\lambda_1-\lambda_2\right)\dfrac{\partial}{\partial x}+\left(C_6u+\lambda_3+\lambda_4\right)\dfrac{\partial}{\partial u}\right],\end{align}
   which are, respectively, the re-numbered exact point symmetries (4.8), and a trivial approximate symmetry $X_\infty$ corresponding to the infinite symmetry set (4.3) of the linear wave equation (4.2). The difference between the freedom in (4.3) and the space spanned by $X_1,\ldots,X_5$ in (4.12) corresponds to unstable point symmetries of the linear wave equation.

2. 2. $Q(u_x)=A\ln(u_x+B)+C$ , where A, B and C are arbitrary constants: here, the nonlinear wave equation (4.6) admits an additional genuine approximate symmetry given by

   (4.13) \begin{equation} X_g=(u+Bx)\dfrac{\partial}{\partial u}-\epsilon \, \dfrac{At}{2}\dfrac{\partial}{\partial t}.\end{equation}

3. 3. $Q(u_x)=u_x$ : the exact symmetry generator (4.3) of the linear wave equation (4.2) reduces to

   (4.14) \begin{equation}X^0=(C_1t+C_2)\dfrac{\partial}{\partial t}+(C_1x+C_3)\dfrac{\partial}{\partial x}+(C_1u+\beta_1+\beta_2)\dfrac{\partial}{\partial u},\end{equation}
   and the $O(\epsilon)$ approximate symmetry components have the form
   (4.15) \begin{align}\xi_1^1 & = \lambda_1(x+t)+\lambda_2(x-t)-\dfrac{1}{4}\displaystyle\int_{}^{t}\left(\beta_2'(x-t)+2z\beta_1''\left(\left(t-2z\right)+x\right)\right)dz, \nonumber\\[5pt]\xi_1^2 & =H(x,t), \qquad \eta_1 =\left( C_4+\dfrac{\beta_2'(x-t)-\beta_1'(x+t)}{4}\right)u+\lambda_3(x+t)+\lambda_4(x-t),\end{align}
   where H(x,t) is an arbitrary solution of the PDEs: $H_t={\xi^1_1}_x$ , $H_x={\xi^1_1}_t+\dfrac{1}{2}\left({\beta_1}_x+{\beta_2}_x\right)$ . In this third case, the point symmetries of the unperturbed linear wave equation (4.2) with arbitrary $\beta_1(x+t)$ and $\beta_2(x-t)$ remain stable and yield genuine approximate symmetries with $O(\epsilon)$ components given by the terms in (4.15) that contain $\beta_1$ and $\beta_2$ . In particular, a genuine approximate symmetry similar to (4.13) is admitted:
   (4.16) \begin{equation} Y_g=x\dfrac{\partial}{\partial u}-\epsilon \, \dfrac{t}{2}\dfrac{\partial}{\partial t}.\end{equation}

4.3. FS approximate point symmetry classification of the wave equation family (4.6)

For the perturbed PDE (4.6) with the solution form

(4.17) \begin{equation}u(x,t)=v(x,t)+\epsilon w(x,t)+o(\epsilon),\end{equation}

the Fuschich-Shtelen system (2.25) reads

(4.18) \begin{equation} v_{tt}=v_{xx},\qquad w_{tt} =w_{xx}+ Q(v_x)v_{xx}.\end{equation}

The basic group of point equivalence transformations of the PDE system (4.18) is given by

(4.19) \begin{align}t & =C_1\tilde{t}+C_2,\qquad\qquad\qquad\ \ x=C_3\tilde{x}+C_4,\nonumber\\v & =C_5\tilde{v}+C_6\tilde{x}+C_7\tilde{t}+C_8,\quad w=C_5\tilde{w}+C_9\tilde{x}+C_{10}\tilde{t}+C_{11},\\c^2 & = \dfrac{C_3^2}{C_1^2},\qquad\qquad\qquad\qquad\ \ Q(u_x) = \dfrac{C_3^2}{C_1^2}\,\tilde{Q}(\tilde{u}_{\tilde{x}}),\nonumber\end{align}

in terms of arbitrary constants $C_1, \ldots, C_{11}$ . We now find exact point symmetries of the system (4.18) that correspond to FS approximate point symmetries of the PDE (4.6). The infinitesimal generator of such symmetries has the form

(4.20) \begin{equation} Z= \lambda^1(x,t,v,w)\dfrac{\partial}{\partial x}+ \lambda^2(x,t,v,w)\dfrac{\partial}{\partial t}+\phi_1(x,t,v,w) \dfrac{\partial}{\partial v}+ \phi_2(x,t,v,w) \dfrac{\partial}{\partial w}. \end{equation}

The solution of the determining equations (2.4) leads to the following classification.

1. 1. $Q(v_x)$ arbitrary:

   (4.21) \begin{align} Z_1 & =\dfrac{\partial}{\partial t},\qquad Z_2=\dfrac{\partial}{\partial x},\qquad Z_3=\dfrac{\partial}{\partial v},\qquad Z_4=t\dfrac{\partial}{\partial v},\nonumber\\Z_5 & = t\dfrac{\partial}{\partial t}+x\dfrac{\partial}{\partial x}+v\dfrac{\partial}{\partial v}+w\dfrac{\partial}{\partial w},\qquad Z_6 = v\dfrac{\partial}{\partial w},\\Z_\infty & =\left(\beta_1(x+t)+\beta_2(x-t)\right)\dfrac{\partial}{\partial w}, \nonumber \end{align}
   where $\beta_1=\beta_1(x+t)$ and $\beta_2=\beta_2(x-t)$ are arbitrary functions. In this general case, no genuine FS approximate symmetries arise. Indeed, the generators $Z_1,\ldots,Z_5$ mimic the exact point symmetry generators (4.8), and $Z_6, Z_\infty$ are trivial FS symmetries. Including the above general symmetries, and modulo the equivalence transformations (4.19), the system (4.18) admits additional point symmetries in the following cases.

2. 2. $Q(v_x)=v_x^s$ , $s\neq 0$ :

   (4.22) \begin{equation} Z_7=v\dfrac{\partial}{\partial v}+(s+1)w\dfrac{\partial}{\partial w}. \end{equation}

3. 3. $Q(v_x)=e^{v_x}:$

   (4.23) \begin{equation} Z_7'=x\dfrac{\partial}{\partial v}+w\dfrac{\partial}{\partial w}. \end{equation}

The symmetries given by $Z_7$ and $Z_7'$ are genuine FS approximate point symmetries of the perturbed PDE (4.6). In particular, $Z_7$ corresponds to a scaling which is different for O(1) and $O(\epsilon)$ components of the approximate solution (4.17).

4.4. Summary

For an arbitrary Q, the perturbed one-dimensional wave equation (4.6) admits five exact symmetries given by (4.8) and it has these five symmetries and a trivial approximate symmetry as BGI approximate symmetries (4.12). For $Q=u_x$ , the equation (4.6) admits (4.8) and an additional exact symmetry given by (4.9). For BGI classification with $Q(u_x)=u_x$ , the PDE (4.6) has an infinite set of BGI approximate symmetries with approximate symmetry components given by (4.15). Note that the exact symmetry generator $Y_6$ in (4.9)

\begin{equation*}Y_6=x\dfrac{\partial}{\partial u}-\epsilon\left(\dfrac{3t}{2}\dfrac{\partial}{\partial t}+x\dfrac{\partial}{\partial x}\right)\end{equation*}

can be obtained from the BGI approximate components (4.15) by taking

\begin{equation*}\beta_1=\dfrac{x+t}{2},\quad \beta_2=\dfrac{x-t}{2},\quad \lambda_1=-\dfrac{11}{16}(x+t),\quad \lambda_2=\dfrac{11}{16}(x-t).\end{equation*}

It follows that the BGI approximate symmetry classification of the wave equation (4.6) includes the exact symmetry classification of (4.6) but corresponds to a subset of exact point symmetries (4.3) of the unperturbed (linear) wave equation (4.2). An additional case appears in the BGI approximate symmetry classification when $Q=A\ln(u_x+B)+C$ , with a corresponding additional approximate symmetry given by (4.13). This case does not arise in the FS symmetry classification.

For an arbitrary Q, the PDE (4.6) admits exact point symmetry generators (4.8) and trivial approximate symmetries given by (4.21). In comparison with the exact and BGI symmetry classifications of (4.6), two different cases appear in FS approximate symmetry classification: $Q=v_x^s$ , $s\neq 0$ and $Q=e^{v_x}$ . For $Q=v_x^s$ , the PDE (4.6) admits an additional FS approximate symmetry given by (4.22). (By contrast, in the exact and BGI symmetry classifications of (4.6), this case appears only when $s=1$ .) For $Q=e^{v_x}$ , a stable exact symmetry $x \,\partial/\partial u$ of the linear wave equation (4.2) yields a genuine FS approximate symmetry of (4.6) given by (4.23).

5.1. Approximate solutions of the PDE (5.2) in the FS framework

Here we use Fuschich-Shtelen approximate symmetries to construct an approximate solution for the PDE (5.2) in the usual FS form

(5.3) \begin{equation}u(x,t)=v(x,t)+\epsilon w(x,t)+o(\epsilon).\end{equation}

In the first-order of precision in $\epsilon$ , the equation (5.2) is equivalent to the Fuschich-Shtelen system (4.18) with $Q(v_x)=v_x^s$ :

(5.4) \begin{equation} v_{tt}=v_{xx},\qquad w_{tt} =w_{xx}+ v_x^sv_{xx},\end{equation}

which admits the symmetry generator (4.22). The corresponding characteristic equations are given by

(5.5) \begin{equation}\dfrac{dt}{0} =\dfrac{dx}{0} =\dfrac{dv}{v}=\dfrac{dw}{(s+1)w}.\end{equation}

Consequently, if v(x,t) is any solution for the first equation of the system (5.4), then the invariant solution following from the characteristic equations (5.5) is given by $w(x,t)=v^{s+1}\phi(x,t)$ . Consider traveling wave solutions of the first equation in (5.4):

(5.6) \begin{equation} v=g\left( x \pm t\right).\end{equation}

Substituting (5.6) and $w= g^{s+1}\phi$ into the second PDE of (5.4), one gets to the PDE in $\phi$

(5.7) \begin{equation} g^{s+1}\left(\phi_{tt}-\phi_{xx}\right)+2(s+1)g^sg'\left(\pm\phi_{t}-\phi_{x}\right)-(g')^{s}g'' =0.\end{equation}

When $s\ne -1$ , the PDE (5.7) has a general solution

\begin{equation*} \phi=g^{-s-1}\left[h \pm \dfrac{t(g')^{s+1}- \int_{}^{t}(g')^{s+1}\left(\pm (2r-t)+x\right)dr}{2(s+1)}\right],\end{equation*}

where $h=h(x,t)$ satisfies $h_{tt}=h_{xx}$ . Similarly, when $s=-1$ , the solution form changes to

\begin{equation*} \phi=h\pm \dfrac{1}{2}t\ln\left(|g'|\right)\,.\end{equation*}

In light of the above results, the higher-order solution part w has the form

(5.8) \begin{equation} w= \begin{cases} h\pm \dfrac{t(g')^{s+1}- \int_{}^{t}(g')^{s+1}\left(\pm (2r-t)+x\right)dr}{2(s+1)}, & s\neq -1, \\[2ex] h\pm \dfrac{t\ln\left(|g'|\right)}{2}, & s=-1. \end{cases}\end{equation}

Finally, when $s\neq -1$ , the perturbed equation (5.2) has the approximate solution (5.3) given by

(5.9a) \begin{align} u(x,t) = g(x\pm t)+\epsilon\left[ h(x,t)\pm \dfrac{ t(g')^{s+1}- \int_{}^{t}(g')^{s+1}\left(\pm (2r-t)+x\right)dr}{2(s+1)}\right]+o(\epsilon)\,.\end{align}

When $s=-1$ , the approximate solution takes the form

(5.9b) \begin{equation} u(x,t)=g(x\pm t)+\epsilon \left(h(x,t)\pm \dfrac{t\ln\left(|g'\left( x\pm t\right)\right|)}{2}\right)+o(\epsilon).\end{equation}

Example 5.1. As a specific example that will be used below, we consider the PDE (5.1) describing shear waves in a fiber-reinforced solid, re-scaled to the form (5.2) with $s=2$ :

(5.10) \begin{equation} u_{tt} = \left(1 + \epsilon \left(u_x\right)^2\right) u_{xx},\end{equation}

and assume in (5.1) that $\beta/\alpha\sim\epsilon\ll 1$ , which corresponds to weak fiber effects. We also choose $v(x,t)=\exp({-(x-t)^2})$ . Then, the solution (5.9a) of the PDE (5.2) with $h=0$ reduces to

(5.11) \begin{align} u(x,t)&=e^{-(x-t)^2}+\dfrac{\epsilon}{6} \biggl[8t\left( x-t\right) ^{3}e^{-3\left( x-t\right) ^{2}}+\dfrac{1}{9}\bigg(\left( 12tx-6t^{2}-6x^{2}-2\right) e^{-3\left( x-t\right) ^{2}}\nonumber\\[5pt] &+\left( 12tx+6t^{2}+6x^{2}+2\right) e^{-3\left( x+t\right)^{2}}\bigg)\biggl] + o(\epsilon). \end{align}

Note that for any fixed t, the approximate solution (5.11) approaches zero as $x \rightarrow \infty$ . Also, for any $x \in ({-}\infty,\infty)$ , the solution (5.11) is bounded as $t \rightarrow \infty$ . The solution (5.11) is not a purely right-traveling wave solution but describes an evolving wave form (see Section 5.2 below). In particular, the PDE (5.10) is known to have breaking wave-type solutions [Reference Cheviakov and Ganghoffer13].

5.2. Numerical simulations of (5.10) and breaking waves

We now compute numerical solutions of the wave equation (5.10) in Example 5.1 in order to model its breaking wave behavior (see Section 5.3 below for details) and provide a reference for comparison with the approximate solutions developed in Section 5.1. Gaussian-type initial conditions corresponding to a right-traveling wave and periodic boundary conditions

(5.12) \begin{equation} u(x,0)=e^{-x^2},\quad u_t(x,0)=2xe^{-x^2},\quad u({-}L,t)=u(L,t),\quad L>0,\quad t \in [0,T],\end{equation}

are posed in the space-time domain $x\in[{-}L,L]$ , $t\geq 0$ , and the equation (5.10) is solved using an explicit finite-difference cross-stencil scheme with constant spatial and temporal steps $\tilde{h}$ , $\tilde{\tau}$ . Following [Reference Cheviakov, Dorodnitsyn and Kaptsov12], we use a conservative finite-difference scheme developed for the PDE (5.10):

(5.13) \begin{align} U_{t\check{t}} & - U_{x\bar{x}} - \epsilon\,\dfrac{U_x^3 - U_{\bar{x}}^3}{\tilde{h}} =0, \quad\tilde{h} = \frac{2L}{M},\quad \tilde{\tau} = \frac{T}{N}, \nonumber\\[3pt]x_m & = -L + m\tilde{h},\qquad\qquad\qquad\! m = 0,\ldots,M, \\[3pt]t_n & = 0 + n\tilde{\tau}, \qquad\qquad\quad\qquad\! n = 0,\ldots,N,\nonumber\end{align}

with $U=U^n_m$ approximating the value of u(x, t) at the mesh node $(x_m, t_n)$ . Here, $U_{t\check{t}}$ and $U_{x\bar{x}}$ represent the second-order central differences, $U_x$ the first-order forward difference and $U_{\bar{x}}$ the first-order backward difference:

(5.14) \begin{align} U_{t\check{t}} & = \dfrac{U^{n+1}_m - 2U^n_m + U^{n-1}_m}{\tilde{\tau}^2}, U_{x\bar{x}} = \dfrac{U^n_{m+1} - 2U^n_m + U^n_{m-1}}{\tilde{h}^2},\nonumber\\U_x & = \dfrac{U^n_{m+1} - U^n_m}{\tilde{h}}, U_{\bar{x}} = \dfrac{U^n_m - U^n_{m-1}}{\tilde{h}}.\end{align}

The numerical solutions provide a good agreement with approximate solutions (5.11), for a broad range of values of $\epsilon$ , from the initial dimensionless time to the time close to the wave breaking. The breaking time increases approximately as $\epsilon^{-1}$ as $\epsilon$ decreases (see Section 5.3 below). Here we present sample computation results for a relatively large value of the small parameter, $\epsilon=0.5$ . The computation is performed from $t=0$ to $t=4$ close to the wave breaking time. A comparison of the numerical solution and the approximate solution (5.11) of PDE (5.10) with initial and boundary conditions (5.12) at several time snapshots is presented in Figure 1(a). The relative difference at the time step $t_n$ between the approximate and numerical solutions is calculated using 2-norms according to the formula

(5.15) \begin{equation} E_n = E(t_n) = \frac{||u_{approx} - u_{num}||_2}{||u_{approx}||_2}\end{equation}

and is shown in Figure 1(b).

Figure 1. (a) Numerical and approximate profiles of u according to the PDE (5.10) ( $\epsilon = 0.5$ ) with initial conditions (5.12) for $L=10$ , $h=0.01$ , $\tau =h/8$ , and $t = 0, 2, 4, 6, 8$ . (b) Relative difference (5.15) between numerical and approximate solutions.

5.3. Estimates of wave breaking times using approximate and numerical solutions

As discussed in [Reference Cheviakov and Ganghoffer13], the nonlinearity $(1 + \epsilon \left(u_x\right)^2)$ leads to greater characteristic speed values at points where $|u_x|$ is larger. This can lead to the intersection of characteristic curves, which corresponds to the wave breaking and the formation of a shock. This behavior can be studied using the method of characteristics. While for linear hyperbolic PDEs, such as the constant-coefficient wave equation $u_{tt} = c^2u_{xx}$ in the simplest case, characteristic curves can be found in terms of explicit formulas such as $x=x_0\pm ct$ and lead to explicit exact solutions, the situation is significantly more complex for nonlinear hyperbolic PDEs. Using the method described in [Reference Haberman25, Reference Zauderer49], one can show that the equation (5.10) can be reduced to the first-order characteristic form

(5.16) \begin{equation} u_t = \pm \dfrac{1}{2\sqrt{\epsilon}}\left(\sqrt{\epsilon}\, u_x\sqrt{1+\epsilon\left(u_x\right)^2} + \ln{\left(\sqrt{\epsilon}\, u_x + \sqrt{1+\epsilon\left(u_x\right)^2}\right)}\right) \end{equation}

on the characteristic curves

(5.17) \begin{equation} \dfrac{dx}{dt} = \pm \sqrt{1+\epsilon\left(u_x\right)^2}\,. \end{equation}

In physical terms, the part of the wave with a time derivative given by (5.16) moves at a finite velocity given by (5.17). The integration of (5.17) yields a constant of integration $x_0$ that corresponds to the point on the characteristic curve where $t = t_0$ is some initial time. Thus, $x = x(x_0, t)$ , from which $u_t(x(x_0, t), t)$ in (5.16) is a function of $x_0$ and t. Similarly, $u_x(x(x_0, t), t)$ in (5.17) is a function of $x_0$ and t.

The formation of a shock where the solution becomes multi-valued takes place when characteristic curves intersect. Without explicit knowledge of $u_x(x, t)$ , no explicit solution $x(x_0,t)$ of (5.17) is available. To estimate the breaking time $T_{b}$ , we use time-progressing linear approximations to characteristic curves, in conjunction with the finite-difference numerical solution of the PDE (5.10) described in Section 5.2. At each time layer $t_n$ in (5.13), linearized characteristics are launched forward in time from each grid point $(x_m, t_n)$ . The smallest time of the intersection of such characteristics estimates the wave breaking time.

For example, when the numerical computation has reached the time layer $t=t_n$ , the linearized characteristics are launched from each spatial grid point $x_m$ , $m = 1,\ldots,M-1$ . In particular, the right-traveling characteristics are approximated by the lines

(5.18) \begin{equation} x = x_m + t\sqrt{1+\epsilon\left(u_x(x_m, t_n)\right)^2}, \end{equation}

where $u_x=U_x$ is the first-order forward finite difference (5.14). To approximate the breaking time $t = \tau$ numerically, we solve (5.18) for the time $\tau$ when two different characteristics intersect. Given two starting points, $x_{m_1}$ and $x_{m_2}$ , we get the system

(5.19) \begin{equation} x = x_{m_1} + \tau\sqrt{1+\epsilon\left(u_x(x_{m_1}, t_n)\right)^2},\quad x = x_{m_2} + \tau\sqrt{1+\epsilon\left(u_x(x_{m_2}, t_n)\right)^2}.\end{equation}

Solving for $\tau$ yields

(5.20) \begin{equation} \tau = \dfrac{x_{m_1} - x_{m_2}}{m_2-m_1}, \end{equation}

where

(5.21) \begin{equation}m_i = \sqrt{1+\epsilon\left(u_x(x_{m_i}, t_n)\right)^2},\qquad i=1, 2\end{equation}

are the slopes of the characteristic lines. We choose $x_{m_1}$ and $x_{m_2}$ to be adjacent grid points, $x_{m_1}+\tilde{h} = x_{m_2}$ . The numerator of (5.20) is constant, so the approximate breaking time corresponds to the largest denominator of (5.20). We determine $x_{m_1}$ corresponding to the largest difference between the slopes $m_1$ and $m_2$ , then solve for $\tau$ .

The meaning of $\tau$ is thus the estimated time from $t_n$ to the wave breaking time $T_{b}$ ; one consequently has an estimate

(5.22) \begin{equation}T_{b}\sim t_n + \tau.\end{equation}

As the wave evolves, the slopes of the linearized characteristics will change and therefore so will $\tau$ . To account for this, at each time step, we repeat the calculation for $\tau$ . We use the first-order forward finite-difference approximation to compute $u_x$ in (5.21) at each time step. Figure 2(a) shows a plot of the value of the time to the break $\tau$ versus the time at which it was calculated, for several values of $\epsilon$ .

Figure 2. (a) Estimates (solid color) and numerical (black dashed) values for the time to break $\tau$ (5.20) for $\epsilon=0.1$ , $0.2$ , … , 1 (right to left). (b) Numerical wave profiles of $u_{xx}$ ( $\epsilon = 0.5$ ) for $t = 0, 1, 2$ .

Alternatively, one can numerically estimate the wave breaking time by defining it as the time when $\text{min}(u_{xx}) \leq \delta$ for some negative number $\delta$ . Choosing for example $\delta = -30$ , and using the second-order central difference approximation to the derivative for $u_{xx}$ in the numerical solution, we calculate the numerical break time for each $\epsilon$ . Using Richardson extrapolation of the break time values found with spatial step sizes $\tilde{h} = 0.01$ and $\tilde{h} = 0.005$ and temporal step size $\tilde{\tau} = 0.00125$ , we found the break times for each $\epsilon$ in the limit as $\tilde{h} \rightarrow 0$ as shown in Figure 2(a). This approach uses the fact that as the wave approaches the breaking time, the second spatial derivative $u_{xx}$ tends to negative infinity (see Figure 2(b)).

We determined the actual, non-linearized characteristic curves by numerically integrating (5.17). The curves are shown in Figure 3.

Figure 3. (a) Characteristic curves found by numerical integration of (5.17) with $\epsilon = 0.5$ . The blue horizontal line is the breaking time as determined by Richardson extrapolation. (b) The same plot as (a) with fewer characteristic lines shown. The thick red characteristic lines correspond to the earliest intersection.

It is interesting that one can also approximately determine the wave breaking time from the approximate solution (5.11) by finding the time when the second spatial derivative $u_{xx}$ of the approximate solution develops an additional root, as shown in Figure 4. This corresponds to an additional inflection point in the wave itself, which can be observed in Figure 1(a) closer to the wave breaking.

The numerically determined break time values ( $\tau_{num}$ ) and the approximate-determined break time values ( $\tau_{approx}$ ) are given in Table 2. Both sets of data are also plotted in Figure 5. The results suggest that the wave breaking time behaves like $\tau \sim \epsilon^{-1}$ .

Table 2. Numerical and approximate wave breaking time estimates for the PDE (5.10) vs. the small parameter values ( $\epsilon$ )

Figure 4. Wave profiles of the approximate solution $u_{xx}$ ( $\epsilon = 0.5$ ) for $t = 0, 1, 3, 3.5$ . Note the development of extra roots as time increases.

Figure 5. Numerical and approximate-derived break times $\tau$ as a function of $\epsilon$ .