The Integrability and Modification to an Auxiliary Function Method for Solving the Strain Wave Equation of a Flexible Rod with a Finite Deformation

Abstract

Our study focuses on the governing equation of a finitely deformed flexible rod with strain waves. By utilizing the well-known Ablowita–Ramani–Segur (ARS) algorithm, we prove that the equation is non-integrable in the Painlevé sense. Based on the bifurcation theory for planar dynamical systems, we modify an auxiliary equation method to obtain a new systematic and effective method that can be used for a wide class of non-linear evolution equations. This method is summed up in an algorithm that explains and clarifies the ease of its applicability. The proposed method is successfully applied to construct wave solutions. The developed solutions are grouped as periodic, solitary, super periodic, kink, and unbounded solutions. A graphic representation of these solutions is presented using a $3D$ representation and a $2D$ representation, as well as a $2D$ contour plot.

1. Introduction

To develop fundamental phenomena and implementations, many non-linear physical structures have been connected to non-linear equations in multiple fields, including plasma physics, fluid dynamics, wave propagation, fluid flow, thermodynamics, non-linear networks, optical including fiber, mechanics, and soil stability. In non-linear sciences, as well as many other disciplines of research, non-linearity plays an important role. In addition, recent years have seen increasing attention directed at finding an exact solution for partial non-linear equations. The reason is that the existence of the exact solutions enables the researchers to explain, understand, and have much knowledge required to analyze the properties of physical phenomena described by non-linear partial differential equations. This motivates scholars to introduce new approaches or develop existing ones. In the literature, there are several efficient and trustworthy approaches for exploring analytic solutions. Several effective and trustworthy approaches have been introduced in the literature for finding the analytic solutions such as the homogeneous balance method [1,2], Hirota bilinear method [3,4], the Backlund transformation method [5,6], the inverse scattering method [7,8], Semi-inverse variational principle [9,10], algebraic method [11,12], the first integral method [13], an extended mapping technique [14], the Riccati–Bernoulli sub-ODE and $exp(G^{\prime }/G)-$expansion method [15], the extended exponential function method [16], Lie symmetry analysis [17,18], a modified F-expansion method [19], unified auxiliary equation method [20], three algebraic methods; $1/G^{\prime }$, modified $G^{\prime }/G^2$ and new extended direct algebraic methods [21], the blackuctive perturbation method [22], bifurcation theory [23,24,25,26,27,28], and for other several methods, see, e.g., [29,30,31,32,33,34,35].

It is well known that non-linear elastic waves have also been used to address a variety of technological matters. There are various sources of non-linearities in solid structures, such as geometrical and physical non-linearities, kinetic non-linearities, and boundary constraints. Solitary and shock wave solutions occur when non-linearity interacts with dispersion or dissipation effects to produce steady traveling wave solutions. Increasingly, problems involving non-linearity are being analyzed and solved qualitatively using the non-linear evolution equation. As a consequence, we analyze some dynamic analysis of

$$\frac{{\partial }^{2}u}{\partial t^2}-c_0^2\frac{{\partial }^{2}u}{\partial x^2}=\frac{{\partial }^2}{\partial x^2}\left[\frac{c_0^2}{2}(3u^2+u^3)+\frac{{\nu }^2{R}^2}{2}\left(\frac{{\partial }^{2}u}{\partial t^2}-c_1^2\frac{{\partial }^{2}u}{\partial x^2}\right)\right],$$

(1)

which is called the strain wave equation [36]. The longitudinal wave velocity and shear wave velocity are denoted by $c_0=\sqrt{\frac{E}{\rho }}$, $c_1=\sqrt{\frac{\mu }{\rho }}$, respectively. $\mu$, $\rho$, E, and $u(x,t)$ refer to the material shear modulus, the rod density per unit volume, elastic modulus, and distribution of the displacement at $x,t$, respectively. The model in Equation (1) explains a double non-linear wave equation for an elastic circular-rod wave-guide. It involves the axial displacement gradient and accounts for the transverse Poisson effect, which causes the longitudinal and shear waves to propagate together. Hamilton’s principle of least action is applied to derive the model [36]. In [37], the authors described the non-linear dispersive dissipative evolution of strain waves in a solid with certain limiting conditions. They constructed a bounded periodic solution and a kink solution. In [36], the authors constructed a solitary wave solution for the governed equation by applying the hyperbolic secant function finite expansion method. In [38], the authors applied the Jacobi elliptic sine function and the third kind of Jacobi elliptic function, finite expansion, to obtain exact periodic solutions of double non-linearity and double dispersion wave equations (cubic non-linear item) and corresponding truncation equations (square non-linear items) of the circular-rod-waveguide. In [39], the authors applied qualitative analysis and utilized the principle of homogeneous balance to solve the governed equation with the Jacobi elliptic function expansion method to construct solitary wave solutions, and they also illustrate that the existence of shock wave solutions is possible under certain conditions. Solely, both solitary and shock wave solutions are found which are commonly not capable of producing an exact periodic solution. Shock, periodic, and solitary solutions of the derived equation have been found by utilizing the Jacobi elliptic function expansion method [40]. Based on the bifurcation analysis, some new solutions, dependence on the initial conditions, influence on the physical parameters, and the existence of quasi-periodic solution are investigated by the current authors [41].

As a result, Section 2 includes the study of the integrability of the governing equation via Painlevé analysis. Section 3 presents the proposed method for developing traveling wave solutions for non-linear partial differential equations, based on bifurcation analysis, and an algorithm to demonstrate its applicability. By applying the newly proposed technique to the governing equation, we obtain some new wave solutions in Section 4. By presenting 3D and 2D representations as well as 2D contours, these solutions are clarified graphically in Section 5. The results are summarized in Section 6.

2. Painlevé Analysis

The integrability of Equation (1) is examined in the current section by utilizing the Painlevé singularity analysis. The ARS algorithm is a powerful technique to study the integrability of partial differential equations [42]. To have a self-contained article, we introduce briefly the ARS algorithm in Appendix A. It is applied successfully in several works, see e.g., [43,44,45,46]. We prove the next theorem

Theorem 1. 

The strain wave Equation (1) of a flexible rod with a finite deformation is not integrable in Painlevé sense.

Proof. 

The generalized Laurent series for the independent variable $u(x,t)$ in the neighborhood of a non-characteristic singular manifold $\psi (x,t)$ with ${\psi }_t\ne 0$ and ${\psi }_x\ne 0$ admits the form

$$u(x,t)={\psi }^p(x,t)\sum _{j=0}^{\infty }{u}_j(x,t){\psi }^j(x,t),$$

(2)

where $u_0(x,t)\ne 0$ and $p\in {\mathbb{Z}}^{-}$ are to be calculated. The leading order term is postulated to be

$$u(x,t)=u_0(x,t){\psi }^p(x,t).$$

(3)

Inserting Equation (2) into Equation (1), p is calculated by balancing the most dominate terms. Hence, we find $p=-1$. The coefficient of the leading term $u_0(x,t)$, which is determined by comparing the coefficient of ${\psi }^{-5}$, is given by

$$u_0(x,t)=\pm \frac{R\nu }{c_0}\sqrt{2(c_1^2{\psi }_x^2-{\psi }_t^2)}.$$

(4)

The resonances, which are the powers of $\psi$ at which the free functions appear in the Laurent series, are determined by plugging the expression

$$u(x,t)=u_0(x,t){\psi }^{-1}(x,t)+u_r(x,t){\psi }^{r-1}(x,t),$$

(5)

into Equation (1) and equating the coefficient of $\psi {(x,t)}^{r-5}$ to zero, we get

$$R^2{\nu }^2{u}_r(r+1)(r-3){(r-4)}^2{\psi }_x^2(c_1^2{\psi }_x^2-{\psi }_t^2)=0,$$

(6)

where $u_0(x,t)$ is eliminated by Equation (5). It is clear that $u_r$ are arbitrary functions if $r=-1,3,4,4$. Notice that there is always a resonance $r=-1$ which corresponds to the arbitrariness of the singular manifold $\psi (x,t)$. As the Equation (1) admits more resonances at $j=4,4$ than the required arbitrary functions, we conclude that Equation (1) fails to satisfy the Painlevé test. Hence, the theorem is proved. □

3. Proposed Method

Due to the governed equation for the elastic rod with a finite deformation, Equation (1) being non-integrable in the Painlevé sense, quasi-analytical methods or numerical analysis are needed to discover special solutions such as wave solutions. This motivate us to develop and modify the current methods to be more systematic and effective for constructing wave solutions for non-linear partial differential equations. In this section, we perform a combination of the bifurcation theory [47] and the auxiliary equation method [48]. In a comparison with previous methods [36,37,38,39,40], this combination will lead to some significant advantages as we clarify below.

3.1. Method Description

Considering that the general form of the partial differential equation takes the form:

$$\mathcal{K}(u,u_t,u_x,u_{tt},u_{xt},u_{xx},\dots )=0.$$

(7)

We employ the transformation

$$u(x,t)=\phi \left(\xi \right),\xi =k(x-\omega t),$$

(8)

where $\xi$ is the wave number while w and k are constants characterizing the wave speed and wave number respectively. Inserting the transformation (8) into Equation (7), we obtain an ordinary differential equation

$$\mathcal{P}(\phi ,-k\omega {\phi }^{\prime },k{\phi }^{\prime },k^2{\omega }^2{\phi }^{″},-k^2\omega {\phi }^{″},k^2{\phi }^{″},\dots )=0,$$

(9)

where primes refer to derivatives with respect to the wave variable $\xi$. This equation is called an reduced equation. We assume Equation (9) has a solution in the form

$$\phi \left(\xi \right)=\sum _{j=0}^N{a}_j\tau {\left(\xi \right)}^j,$$

(10)

where $a_j$ are constants required to be determined and $N\in {\mathbb{Z}}^{+}$ is calculated by balancing the highest power of the non-linear term with the highest order derivative terms in Equation (9), besides our assumption that $\tau \left(\xi \right)$ satisfies the auxiliary differential equation

$${\tau }^{\prime }=\sqrt{\beta ({\tau }^4-\delta {\tau }^2)+2e},$$

(11)

where $\beta ,\delta ,e$ are arbitrary constants. Thus, the key to constructing a solution for the main partial differential Equation (7) is to find a solution for the reduced Equation (9). Thus, separating the variables of Equation (11) gives the one-differential form

$$\frac{d\tau }{\sqrt{P_4\left(\tau \right)}}=d\xi ,$$

(12)

where

$$P_4\left(\tau \right)=\beta ({\tau }^4-\delta {\tau }^2)+2e.$$

(13)

The parameters $\beta ,\delta ,e$ must fall within a suitable range for the integration of both sides of Equation (12). There are various ways to accomplish this, including the discriminating approach and bifurcation. The bifurcation theory is more significant because it enables us to specify the kind of solutions before constructing them through type of the phase orbits. For instance, the existence of periodic orbits, homoclinic orbits, and heteroclinic orbits indicates the presence of periodic, solitary and kink (or anti-kink) solutions. It has been applied successfully in several works, see e.g., [49,50,51] as well as being applied to find the solutions for the stochastic partial differential equation [52,53] and investigating the chaotic behavior for perturbed partial differential equations containing perturbed periodic term [54]. In this work, we apply the bifurcation theory [47] by using the Hamiltonian concepts to find the solution of Equation (11) that can be expressed as

$$\frac{1}{2}{\tau }^{{\prime }^2}-\frac{\beta }{2}({\tau }^4-\delta {\tau }^2)=e.$$

(14)

Equation (14) is equivalent to the conserved quantity for a one-dimension conservative Hamiltonian system which is derived from the Hamiltonian function

$$H=\frac{1}{2}{\tau }^{{\prime }^2}+U\left(\tau \right),$$

(15)

where U is the potential function

$$U\left(\tau \right)=-\frac{\beta }{2}({\tau }^4-\delta {\tau }^2).$$

(16)

Notice that the constant e does not appear in the potential function (16) and moreover, it enters the Hamiltonian (15) linearly. Hence, it can be taken as the energy constant for the Hamiltonian (15) [55,56].

Now, we are going to apply the bifurcation analysis for the Hamiltonian (15). Firstly, we find the equilibrium points for the Hamiltonian (15) which appear as the critical points for the potential function (16), i.e., they are $({\tau }_0,0)$, where ${\tau }_0$ is a solution for $U^{\prime }\left({\tau }_0\right)=-\beta \tau (2{\tau }_0^3-\delta \tau )=0$. Thus, if $\delta >0$, there exist three equilibrium points ${\mathcal{E}}_0=(0,0),{\mathcal{E}}_{1,2}=(\pm \sqrt{\frac{\delta }{2}},0)$ while if $\delta <0$, there is a single equilibrium point ${\mathcal{E}}_0=(0,0)$. The equilibrium point can be assorted according to whether it is either the maximum or minimum point for the potential function (16). Therefore, we calculate

$$U^{″}\left({\mathcal{E}}_0\right)=\beta \delta ,U^{″}\left({\mathcal{E}}_{1,2}\right)=-2\beta \delta .$$

(17)

Assuming the condition of the equilibrium point is met, we find:

(a)

If $\delta <0$, the Hamiltonian (15) has a single equilibrium point ${\mathcal{E}}_0=(0,0)$ which is either center point if $\beta <0$ or saddle point if $\beta >0$. The phase portrait for these cases are depicted in Figure 1a,b.

(b)

If $\delta >0$, the Hamiltonian (15) owns three equilibrium points ${\mathcal{E}}_{0,1,2}$. For $\beta >0$, ${\mathcal{E}}_0$ is a center point while ${\mathcal{E}}_{1,2}$ are saddle points while for $\beta <0$, the two equilibrium points ${\mathcal{E}}_{1,2}$ are center and ${\mathcal{E}}_0$ is saddle point. The phase portrait for both cases of $\beta$ is depicted by Figure 1c,d.

The value of the conserved quantity (14) at the equilibrium points are

$$e_0=e\left({\mathcal{E}}_0\right)=0,e_1=e\left({\mathcal{E}}_{1,2}\right)=\frac{\beta {\delta }^2}{8}.$$

(18)

Depending on the bifurcation conditions on the parameters $\beta ,\delta ,e$, we can integrate both sides of Equation (12) to construct the possible solutions for the auxiliary Equation (9). These solutions are accumulated in

Table 1. Possible bounded real solutions for Equation (9).

Case	$\mathit{\delta }$	$\mathit{\beta }$	Interval of	Range of e	Solution
			Real Propagation
1.	−	−	$(-|r_2|,0)\bigcup (0,|r_2\left|\right)$	$(e_0,\infty )$	$\tau \left(\xi \right)=\pm A\mathrm{sd}({\Omega }_1\xi ,k_1).$
2.	+	+	$(-|r_2|,0)\bigcup (0,|r_2\left|\right)$	$(e_0,e_1)$	$\tau \left(\xi \right)=\pm r_2\mathrm{sn}(r_1\sqrt{\beta }\xi ,\frac{r_2}{r_1}).$
3.	+	−	$(-|r_2|,-|r_1\left|\right)\bigcup \left(\right|r_1|,|r_2\left|\right)$	$(e_1,e_0)$	$\tau \left(\xi \right)=\pm r_2\mathrm{dn}(r_2\sqrt{-\beta }\xi ,\frac{1}{k_1}).$
4.	+	−	$(-|r_2|,|r_2\left|\right)$	$(e_0,\infty )$	$\tau \left(\xi \right)=\pm r_2\mathrm{cn}({\Omega }_1\xi ,k_1).$

where $\mathrm{sn}(z,k), \mathrm{dn}(z,k), \mathrm{cn}(z,k), \mathrm{sd}(z,k)=\frac{\mathrm{sn}(z,k)}{\mathrm{dn}(z,k)}$ are Jacobian elliptic functions [57], $A=\frac{r_2\sqrt{-r_1^2}}{\sqrt{r_2^2-r_1^2}}$, ${\Omega }_1=\sqrt{-\beta (r_2^2-r_1^2)}$, $k_1=\sqrt{\frac{r_2^2}{r_2^2-r_1^2}}$, and $r_2^2>0$, $r_1^2<0$ in Cases 1 and 4.

,

Table 2. Soliton and kink (or anti-kink) solutions for Equation (9).

Case	$\mathit{\delta }$	$\mathit{\beta }$	Interval of	Range of e	Solution
			Real Propagation
1.	+	+	$(-|r_1|,0)\bigcup (0,|r_1\left|\right)$	$e_1$	$\tau \left(\xi \right)=\pm \sqrt{\frac{\delta }{2}}\tan \mathrm{h}\left(\sqrt{\frac{\delta \beta }{2}}\xi \right).$
2.	+	−	$(-|r_2|,0)\bigcup (0,-|r_2\left|\right)$	$e_0$	$\tau \left(\xi \right)=\pm \sqrt{\delta }\sec \mathrm{h}\sqrt{-\delta \beta }\xi .$

and

Table 3. Possible unbounded real solutions for Equation (9).

Case	$\mathit{\delta }$	$\mathit{\beta }$	Interval of	Range of e	Solution
			Real Propagation
1.	−	+	$(-\infty ,-|r_1\left|\right)\bigcup \left(\right|r_1|,\infty )$	$(-\infty ,e_0)$	$\tau \left(\xi \right)=\pm A_2\mathrm{ds}({\Omega }_2\xi ,k_2).$
2.	−	+	${\mathbb{R}}^{*}$	$e_0$	$\tau \left(\xi \right)=\pm \sqrt{-\delta }\csc \mathrm{h}\sqrt{-\delta \beta }\xi .$
3.	−	+	$\mathbb{R}$	$(e_0,\infty )$	$\tau \left(\xi \right)=\pm \sqrt{-r_1^2}\mathrm{tn}(\sqrt{-r_2^2\beta }\xi ,\frac{1}{k_1}).$
4.	+	+	$(-\infty ,-|r_1\left|\right)\bigcup \left(\right|r_1|,\infty )$	$(-\infty ,e_0)$	$\tau \left(\xi \right)=\pm r_1\mathrm{nc}({\Omega }_2\xi ,k_2).$
5.	+	+	$(-\infty ,-|r_1\left|\right)\bigcup \left(\right|r_1|,\infty )$	$e_0$	$\tau \left(\xi \right)=\pm \sqrt{\delta }\sec \sqrt{\delta \beta }\xi .$
6.	+	+	$(-\infty ,-|r_1\left|\right)\bigcup \left(\right|r_1|,\infty )$	$(e_0,e_1)$	$\tau \left(\xi \right)=\pm r_1\mathrm{dc}(r_1\sqrt{\beta }\xi ,\frac{r_2}{r_1}).$
7.	+	+	$(-\infty ,-|r_1\left|\right)\bigcup \left(\right|r_1|,\infty )$	$e_1$	$\tau \left(\xi \right)=\pm \sqrt{\frac{\delta }{2}}\cot \mathrm{h}\left(\sqrt{\frac{\delta \beta }{2}}\xi \right).$
8.	+	+	$\mathbb{R}$	$(e_1,\infty )$	$\tau \left(\xi \right)=\pm \frac{|r_1\left|\mathrm{nd}\right(|r_1|\sqrt{\beta }\xi ,k_3)}{\mathrm{tn}\left(\right|r_1|\sqrt{\beta }\xi ,k_3)}.$

where $\mathrm{ds}(z,k)=\frac{\mathrm{dn}(z,k)}{\mathrm{sn}(z,k)}, \mathrm{tn}(z,k)=\frac{\mathrm{sn}(z,k)}{\mathrm{cn}(z,k)}, \mathrm{nc}(z,k)=\frac{1}{\mathrm{cn}(z,k)}, \mathrm{dc}(z,k)=\frac{\mathrm{dn}(z,k)}{\mathrm{cn}(z,k)}, \mathrm{nd}(z,k)=\frac{1}{\mathrm{dn}(z,k)}$ are Jacobian elliptic functions [57], $A_2=\sqrt{r_1^2-r_2^2}$, ${\Omega }_2=\sqrt{\beta (r_1^2-r_2^2)}$, $k_2=\sqrt{\frac{-r_2^2}{r_1^2-r_2^2}}$, $k_3=\frac{Im\left(r_1\right)}{|r_1|}$, and $r_1^2>0$, $r_2^2<0$ in Cases 1 and 4, while $r_1^2<0$, $r_2^2<0$ in Case 3.

to avoid ambiguity.

• Equation (13) can be expressed as

  $$P_4\left(\tau \right)=\beta ({\tau }^2-r_1^2)({\tau }^2-r_2^2),$$

  (19)

  where

  $$r_1^2=\frac{\beta \delta +\sqrt{{\beta }^2{\delta }^2-8\beta e}}{2\beta },r_2^2=\frac{\beta \delta -\sqrt{{\beta }^2{\delta }^2-8\beta e}}{2\beta }.$$

  (20)

The zeros of $P_4\left(\tau \right)$ are then $\tau =\pm \sqrt{r_1^2}$ and $\tau =\pm \sqrt{r_2^2}$ which may be real or complex according to the values of the parameters $\beta ,\delta$ and e.

We sum up our proposed method in the following algorithm to demonstrate how easily it can be applied.

Algorithm

An algorithm is presented to demonstrate the applicability of the proposed method. We consider a non-linear partial differential equation in the form (7). To construct possible wave solutions, we follow the steps shown below:

• Step 1: Applying the wave transformation (8) to Equation (7), we obtain an ordinary differential Equation (9) which is the reduced equation.
• Step 2: Assuming the reduced Equation (9) has a solution

  $$\phi \left(\xi \right)=\sum _{j=0}^N{a}_j\tau {\left(\xi \right)}^j,$$

  (21)

  where $\tau \left(\xi \right)$ is a solution of the auxiliary equation

  $${\tau }^{\prime }\left(\xi \right)=\sqrt{\beta ({\tau }^4-\delta {\tau }^2)+2e}.$$

  (22)
• Step 3: Determining the positive integer N by balancing the highest power of the non-linear term with the highest order derivative terms in the auxiliary Equation (9). After some calculations, the highest order derivative term ${\phi }^{″}$ and highest power of the non-linear term ${\phi }^m$ with positive integer $m>1$, in the reduced equation are given by

  $$O\left({\phi }^{″}\right)=N+2,O\left({\phi }^m\right)=mN.$$

  (23)

Equating the last expression, we obtain

$$N=\frac{2}{m-1}.$$

(24)

Because N is a positive integer, we have two possible cases:

• Case I: If $m=2,3$, then we have $N=2,1$, respectively, and go directly to the next step.
• Case II: If $m\ge 4$, we perform the transformation $\phi \left(\xi \right)=G{\left(\xi \right)}^{2/(m-1)}$ to the reduced equation and assume its solution admits the form (21), i.e.,

  $$G\left(\xi \right)=\sum _{j=0}^N{a}_j\tau {\left(\xi \right)}^j.$$
  Then, go to the next step.

• Step 4: Inserting the assumed solution (21) in Equation (7), using the Equation (22), and equating all the coefficients of all powers of $\tau$ in the resulting polynomial to zero, we obtain a non-linear algebraic system in $a_j$ and the physical parameters of the problem.
• Step 5: Solving the last system to find $a_j$ in terms of the physical parameters and utilizing Table 1, Table 2 and Table 3 to construct the solutions for the reduced equation and consequently, to the non-linear partial differential equation under consideration.

The proposed method has many advantages as a result of its dependence on bifurcation theory in its formulation. This is illustrated in the following:

(a)

It is utilized only to construct real (not complex) wave solutions for the given NPDE, because these solutions are constructed by integrating the conserved quantities along all possible intervals of real wave propagation. Moreover, the entering of the concept of the interval of real wave propagation intervals enables us to construct all possible wave solutions that are completely different from mathematical and physical points of view. Let us clarify this point. For the choice $\delta >0,\beta >0$ and $e\in ]e_0,e_1[$, there are two solutions. One of them is periodic as illustrated by row 2 in Table 1 while the other is unbounded as outlined by row 6 in Table 3 row. Thus, we have two completely different solutions from mathematical and physics points of view for the same conditions of the parameters.

(b)

It enables us to know in advance the types of solutions. For example in Table 1, the first three cases are periodic solutions because they are related to periodic orbits in phase portrait while the fourth case is a super periodic solution since it corresponds to the super periodic orbit, see Figure 1. Therefore, the assumed solution can be viewed as a combination of periodic solutions or super periodic solutions, respectively. For this reason, we know the obtained solution will be periodic or supper periodic.

(c)

Thanks to employing bifurcation analysis, we were able to isolate all possible bounded solutions that are required and significant in real-world problems.

4. Application

The current section aims to clarify how to apply the proposed method by following the introduced algorithm. Based on the first step, we insert the wave transformation (8) into the Equation (1), and we get

$$({\omega }^2-{c_0}^2)\frac{d^{2}u}{d{\xi }^2}-\frac{c_0^2}{2}\frac{d^2}{d{\xi }^2}(3u^2+u^3)-\frac{{\nu }^2{R}^2{k}^2}{2}\frac{d^2}{d{\xi }^2}[({\omega }^2-{c_1}^2)\frac{d^{2}u}{d{\xi }^2}]=0.$$

(25)

Integrated Equation (25) twice with respect to $\xi$ and neglecting the integration constants gives

$$\frac{d^{2}u}{d{\xi }^2}=cu+3a{u}^2+a{u}^3,$$

(26)

where

$$c=\frac{2({\omega }^2-c_0^2)}{{\nu }^2{R}^2{k}^2({\omega }^2-c_1^2)}\mathrm{and}a=-\frac{c_0^2}{{\nu }^2{R}^2{k}^2({\omega }^2-c_1^2)}.$$

(27)

According to the second step, we postulate Equation (26) has a solution in the form

$$u\left(\xi \right)=\sum _{j=0}^N{a}_j\tau {\left(\xi \right)}^j,$$

(28)

where $a_j$ is a constant required to be determined besides $\tau \left(\xi \right)$ satisfies a first order differential Equation (22). Following Equation (24), we get $N=1$. Hence, the assumed solution becomes

$$u\left(\xi \right)=a_0+a_1\tau \left(\xi \right).$$

(29)

Now, we are going to determine the constants $a_0,a_1,\beta ,\delta ,e$ in terms of the physical parameters. Inserting the expression (29), taking into account the Equation (22), and equating the coefficient of all powers of $\tau \left(\xi \right)$ to zero, we get the following algebraic system

$$\left\{\begin{array}{c}{a}_1(2\beta -a{a}_1^2)=0,\hfill \\ a{a}_1^2(a_0+1)=0,\hfill \\ a_1[c+\beta \delta +6a{a}_0+3a{a}_0^2]=0,\hfill \\ a_0[c+3a{a}_0+a{a}_0^2]=0.\hfill \end{array}\right.$$

(30)

The only working solution of the algebraic system (30) is

$$a=\frac{c}{2},a_0=-1,\delta =\frac{a}{\beta },a_1=\sqrt{\frac{2\beta }{a}},$$

where $a\beta >0$ is a condition for the solution (29) to be real. Thus, we consider the following cases:

• Case A: For $a>0$, we have $\beta >0$ and consequently $\delta >0$. For bounded solutions, the only working case is Case 2 in Table 1. Thus, Equation (1) has the following solution if $e\in ]0,\frac{a}{8\beta }[$

  $$u(x,t)=-1\pm \sqrt{\frac{2\beta }{a}}{r}_2\mathrm{sn}(r_1\sqrt{\beta }k(x-\omega t),\frac{r_1}{r_2}).$$

  (31)

The solution (31) is a new periodic solution for Equation (1). For Table 2, the first case in which $\delta >0,\beta >0$ is the only working case in this table. Thus, Equation (1) has the solution

$$u(x,t)=-1\pm \sqrt{\frac{2\beta }{a}}\sqrt{\frac{\delta }{2}}\tan \mathrm{h}\sqrt{\frac{\delta \beta }{2}}k(x-\omega t).$$

(32)

The solution (32) is a new kink solution for Equation (1). For unbounded wave solutions, the working cases in Table 3 are Cases 4–7 and they are classified according to values of e. Thus, we have

1.

If $e<0$, Equation (1) owns the following new solution

$$u(x,t)=-1\pm r_1\sqrt{\frac{2\beta }{a}}\mathrm{nc}(\sqrt{\beta (r_1^2-r_2^2)}k(x-\omega t),\sqrt{\frac{r_2^2}{r_2^2-r_1^2}}).$$

(33)

2.

If $e=0$, we introduce a novel solution for Equation (1) in the form

$$u(x,t)=-1\pm \sqrt{\frac{2\beta \delta }{a}}\sec \sqrt{\delta \beta }k(x-\omega t).$$

(34)

3.

If $e\in ]0,\frac{\beta {\delta }^2}{8}[$, the governing Equation (1) has a new solution

$$u(x,t)=-1\pm r_1\sqrt{\frac{2\beta }{a}}\mathrm{dc}(r_{1}k\sqrt{\beta }(x-\omega t),\frac{r_2}{r_1}).$$

(35)

4.

If $e>\frac{\beta {\delta }^2}{8}$, then

$$u(x,t)=-1\pm r_1\sqrt{\frac{2\beta }{a}}\frac{\mathrm{nd}(|r_1|\sqrt{\beta }\xi ,k_3)}{\mathrm{sc}(|r_1|\sqrt{\beta }\xi ,k_3)}.$$

(36)

is a new solution for Equation (1).

• Case B: For the choice $a<0$, $\beta$ must be negative, i.e., $\beta <0$ and hence $\delta >0$. Therefore, the only working cases in Table 1 are Case 3 and Case 4. Let us consider them individually.

1.

If $e\in (\frac{a}{8\beta },0)$, Equation (1) has a new solution in the form

$$u(x,t)=-1+\sqrt{\frac{2\beta }{a}}{r}_2\mathrm{dn}(r_2\sqrt{-\beta }k(x-\omega t),\sqrt{1-\frac{r_1^2}{r_2^2}}).$$

(37)

2.

If $e\in ]0,\infty [$, then Equation (1) has the solution

$$u(x,t)=-1+\sqrt{\frac{2\beta }{a}}{r}_2\mathrm{cn}(k\sqrt{-\beta (r_2^2-r_1^2)}(x-\omega t),\frac{r_2}{\sqrt{r_2^2-r_1^2}}).$$

(38)

As we see also in Table 2, the only working case is the second one. Hence, if $e=e_0=0$ Equation (1) has the solution

$$u(x,t)=-1\pm \sqrt{\frac{2\beta \delta }{a}}\sec \mathrm{h}\sqrt{-\delta \beta }k(x-\omega t),$$

(39)

which is a new solitary wave solution for Equation (1).

5. Physical Interpretations

This section aims to illustrate the obtained solutions graphically by introducing $2D$ and $3D$ graphics, as well as the contour plots for different kinds of the obtained solutions using Matlab and Mathematica programs. We presume that $k=0.1,\omega =0.3$ in all the subsequent computations.

For real solutions, the selected values of a and $\beta$ must satisfy the condition $a\beta >0$ and consequently $\delta >0$. Hence, there are two possible choices of the two parameters, a and $\beta$. They are either $a>0,\beta >0$ or $a<0,\beta <0$. Let us consider each of them individually:

Firstly, We consider the solutions introduced in Case A. We assume $c=0.8$ and $\beta =1$, and so $a=0.4$ and $\delta =0.4$. For these values, we have $e_1=0.02$ and the roots of the quartic polynomial (13) satisfy $r_1^2=0.2-\sqrt{0.08}$ and $r_2^2=0.2+\sqrt{0.08}$. By taking $e=0.01$, the periodic solution (31) for the Equation (1) is depicted in Figure 2. It also clarifies one of the advantages of the proposed method, which is the ability to determine the type of solution before it is constructed, simply by knowing the parameters.

Preserving the same values for c, and $\beta$ and choosing different values of e, we will construct several solutions which are completely different from mathematical and physical points of views.

1.

For $e=e_1=0.02$, Figure 3 shows the graphical representation of the solution (32) by presenting the $3D$ and $2D$ representation besides the $2D-$ contour plot.

2.

For $e=-0.5\in ]-\infty ,0[$, we have $r_1^2=0.2-\sqrt{1.1}$ and $r_2^2=0.2+\sqrt{1.1}$. Hence, the quartic polynomial (13) has two real roots and two imaginary roots. Therefore, the solution (33) can be visualized by Figure 4 which includes the $3D$ and $2D$ graphical representation of the solution in addition to the $2D-$ contour plot.

Secondly, we consider the solutions which are presented in Case B. We postulate $c=-0.8$ and $\beta =-1$, so, $a=-0.4$ and $\delta =0.4$. Depending on the value of e, the possible solutions can be constructed.

1.

If $e=0.001\in (0,\infty )$, Equation (1) admits the super periodic solution (38) which is illustrated by Figure 5.

2.

If $e=0$, Equation (1) has the one soliton solution (39) which is clarified by Figure 6.

6. Conclusions

The current work focused on studying the governing equation describing the strain wave equation of a flexible rod with a finite deformation. Based on the singularity analysis, we have proved the non-integrability of the governing equation. We present a new modification to the auxiliary equation method for constructing wave solutions for a large class of non-linear partial differential equations. The formulation of this method relies on the bifurcation theory for planar integrable systems. The formulation procedures furnish more effective advantages for the new method. These advantages are summarized as:

(a)

It is only applied to construct real (not complex) wave solutions for a large class of partial differential equations by integrating the conserved quantities along the intervals of real wave propagation. Moreover, with the same conditions on the parameters and different intervals of real wave propagation, distinct solutions from mathematical and physical points of view are constructed. For instance, for the choice $\delta >0,\beta >0$ and $e\in ]e_0,e_1[$, there are two solutions. One of them is periodic as illustrated by row 2 in Table 1, while the other is unbounded as outlined by row 6 in Table 3 row.

(b)

It enables us to know the kinds of solutions before establishing them.

(c)

We can isolate all bounded wave solutions, which is useful in real applications.

This method has been applied to the governing equation aiming to illustrate its applicability and effectiveness to present new solutions for the governing equation. Some novel wave solutions that are assorted into periodic, super periodic, kink, and solitary wave solutions have been introduced. By introducing $2D$, $3D$, and contour plots, some of these solutions are graphically clarified.