*Article* **Numerical Investigation of Nonlinear Shock Wave Equations with Fractional Order in Propagating Disturbance**

**Jiahua Fang 1, Muhammad Nadeem 2,\*, Mustafa Habib <sup>3</sup> and Ali Akgül <sup>4</sup>**


**Abstract:** The symmetry design of the system contains integer partial differential equations and fractional-order partial differential equations with fractional derivative. In this paper, we develop a scheme to examine fractional-order shock wave equations and wave equations occurring in the motion of gases in the Caputo sense. This scheme is formulated using the Mohand transform (MT) and the homotopy perturbation method (HPM), altogether called Mohand homotopy perturbation transform (MHPT). Our main finding in this paper is the handling of the recurrence relation that produces the series solutions after only a few iterations. This approach presents the approximate and precise solutions in the form of convergent results with certain countable elements, without any discretization or slight perturbation theory. The numerical findings and solution graphs attained using the MHPT confirm that this approach is significant and reliable.

**Keywords:** Mohand transform; homotopy perturbation method; shock wave equation

#### **1. Introduction**

In recent decades, various fractional models in science and technology have been designed in terms of nonlinear partial differential equations (PDEs), such as plasma physics, fluid dynamics, nonlinear optics, quantum mechanics, solid-state physics, mathematical biology and chemical kinetics [1–3]. Fractional differential equations have been widely used to model complex phenomena in various branches of science and engineering, such as wave propagation, lattice vibration, optical fiber, nanotechnology and biology [4,5]. The scientific theory of shock waves played a role in the problems of motion of gases and compressible liquids in the second half of the 19th century. They are described by nonlinear hyperbolic PDEs and can be written in their simplest form as [6]

$$D^a\_\wp \theta(\mathfrak{F}, \wp) + f\Big(\mathfrak{F}(\mathfrak{F}, \wp)\Big)\_{\mathfrak{F}} = 0, \qquad \mathfrak{F} \in \mathbb{R}, \ \wp > 0 \tag{1}$$

with the initial condition

$$
\theta(\odot, 0) = \theta\_0(\odot), \qquad \odot \in \mathbb{R}. \tag{2}
$$

The shock wave equation is a nonlinear PDE and has given an important contribution to various studies, such as those of explosions, traffic flow, glacier waves and airplanes breaking the sound barrier. Goswami et al. [7] used an effective scheme based on the Sumudu transform and the homotopy perturbation method to find the numerical solutions of time fractional Schrodinger equations with harmonic oscillator. Singh and Gupta [8] presented the homotopy perturbation method (HPM) to examine the numerical solution of the time fractional shock wave equation and wave equation. Allan and Khaled [9] employed

**Citation:** Fang, J.; Nadeem, M.; Habib, M.; Akgül, A. Numerical Investigation of Nonlinear Shock Wave Equations with Fractional Order in Propagating Disturbance. *Symmetry* **2022**, *14*, 1179. https://doi.org/ 10.3390/sym14061179

Academic Editor: Francisco Martínez González

Received: 3 May 2022 Accepted: 4 June 2022 Published: 8 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the Adomian decomposition method to provide the analytical solution of the shock wave equation. Das and Kumar [10] proposed a method for calculating the approximate solution of the shock wave equation and shallow water equation with time derivatives. Later, many researchers [11–14] have developed different strategies to achieve the approximate solution of nonlinear shock wave equations of fractional order.

A differential problem of symmetry is a modification that generates the differential equation continuously in such a way that these symmetries can help to achieve the solution of the differential equation. Solving these equations is sometimes easier than solving the Volterra integro-differential equations [15]. Symmetries can be identified by solving a set of connected ordinary differential equations. PDEs of fractional order are PDEs whose symmetry condition is separated into two segments of integer order and fractional order, and the linear scheme of fractional PDEs reveals a wide dimensional trivial solution continuously. Various numerical and analytical approaches have been demonstrated to attain the semianalytical solution of nonlinear PDEs, such as the (G /G)-expansion method [16], the neural network approach [17], the variational iteration method [18], the Exp-function method [19], the homotopy perturbation method [20], the homotopy analysis method [21], residual power series [22], the residual power series method [23], the quasi-wavelet method [24], the Haar wavelet method [25] and the two-scale approach [26]. New developments of the HPM can be found in [27,28].

The aim of this paper is to present the idea of the MT coupled with the HPM for the numerical investigation of nonlinear shock wave equations of fractional order. The obtained results are expressed in terms of series with easily computable components. This series solution converges to the exact solution rapidly. This study is summarized as follows: In Section 2, we demonstrate some basic preliminary concepts. In Section 3, a new strategy is sorted out to handle nonlinear expressions. In Section 4, some numerical examples are demonstrated to determine the competence of the proposed strategy, and at last, some results are discussed with our conclusions in Sections 5 and 6.

#### **2. Preliminary Concepts**

**Definition 1.** *Let ϑ*(℘) *be a function precise for* ℘ ≥ 0 *[29]; then, we have*

$$\mathcal{X}[\theta(\wp)] = V(r) = \int\_0^\infty \theta(\wp) e^{-r\wp} d\wp.$$

*which is said to be a Laplace transform, where* ℘ *is a function (i.e., a function of the time domain), defined on* [0, ∞)*, to a function of r (i.e., of the frequency domain).*

**Definition 2.** *If V*(*r*) *symbolizes the Laplace transform of ϑ*(℘)*, then*

$$\vartheta(\wp) = \mathcal{X}^{-1} V(r)\_{\prime}$$

*is termed as the inverse Laplace transform of V*(*r*)*.*

**Definition 3.** *Mohand and Mahgoub [30,31] developed the MT to facilitate ordinary and PDEs. Let the MT be expressed with the help of operator* M(.)*. Then* =⇒

$$\mathcal{A}\left[\theta(\wp)\right] = \mathcal{S}(r) = r^2 \int\_0^\infty \theta(\wp) e^{-r\wp} d\wp, k\_1 \le r \le k\_2, \qquad k\_1, k\_2 \in \mathbb{N}$$

*where k*<sup>1</sup> *and k*<sup>2</sup> *are constants. On the other hand, if S*(*r*) *is the MT of ϑ*(℘)*, then ϑ*(℘) *is said to be the inverse of S*(*r*)*, so*

$$\mathcal{A} \mathcal{A}^{-1} \{ S(r) \} = \theta(\wp) \quad \implies \quad \mathcal{A}^{-1} \text{ is the inverse MT.}$$

*One may see that the Laplace transform and the Mohand transform differ in the function of r (i.e., the frequency domain).*

**Lemma 1.** *The MT of a function of fractional order is [32]*

$$\mathcal{M}\left\{\mathcal{S}^{\mathfrak{A}}(\wp)\right\} = r^{\mathfrak{A}}S(r) - \sum\_{k=0}^{n-1} \frac{\mathfrak{u}^{k}(0)}{r^{k} - (\mathfrak{a}+1)}, \quad 0 < \mathfrak{a} \le n$$

**Proposition 1.** *Let* M{*ϑ*(℘)} = *S*(*r*)*; then, the MT of ϑ* (℘) *has the following properties:*

$$\begin{aligned} \text{(a)} \quad \mathcal{A}\{\;\theta'(\wp)\} &= r\mathbb{S}(r) - r^2\theta(0);\\ \text{(b)} \quad \mathcal{A}\{\;\theta''(\wp)\} &= r^2\mathcal{S}(r) - r^3\theta(0) - \theta^2\theta'(0);\\ \text{(c)} \quad \mathcal{A}\{\;\theta^n(\wp)\} &= r^n\mathcal{S}(r) - r^{n+1}\theta(0) - r^n\theta'(0) - \dots - r^2\theta^{n-1}(0). \end{aligned}$$

**Definition 4.** *The fractional derivative [15] in the Caputo sense is*

$$D\_{\tau}^{\alpha}\theta(\mathcal{G},\wp) = \begin{cases} \frac{\partial^{n}\theta(\mathcal{G},\wp)}{\partial\wp^{n}}, & \alpha \in \mathbb{N} \\\\ \frac{1}{\Gamma(n-\alpha)}\int\_{0}^{\wp}(t-\phi)^{n-\alpha-1}\theta^{n}(\phi)\partial\phi, & n-1 < n < n \end{cases}$$

#### **3. Idea of MHPT**

In this section, we construct the idea of the MHPT to find the approximate solution of fractional problems. Therefore, consider a differential equation of fractional order

$$D^{\mathfrak{A}}\_{\wp}\theta(\mathfrak{F},\wp) + R\theta(\mathfrak{F},\wp) + N\theta(\mathfrak{F},\wp) = \mathfrak{g}(\mathfrak{F},\wp),\tag{3}$$

$$
\theta(\odot, 0) = h(\odot),
\tag{4}
$$

where *D<sup>α</sup>* <sup>℘</sup> <sup>=</sup> *<sup>∂</sup><sup>α</sup> <sup>∂</sup>*℘*<sup>α</sup>* is an operator with fractional order *<sup>α</sup>*; *<sup>ϑ</sup>* is the function in the direction of spital and time ℘; *R* is the linear; *N* represents the nonlinear differential operator; and *g*(, ℘) is the source term. Employing the MT in Equation (3), we obtain

$$\mathcal{A}\mathcal{A}\left[D^{a}\_{\phi}\theta(\odot,\wp) + R\theta(\odot,\wp) + N\theta(\odot,\wp)\right] = \mathcal{A}\left[\mathcal{g}(\oslash,\wp)\right],\tag{5}$$

using the differentiation property of the MT, we obtain

$$r^{\mu}\left[R(r) - r\theta(0)\right] = -\mathcal{A}\left[R\theta(\odot,\wp) + N\theta(\odot,\wp)\right] + \mathcal{A}\left[\mathcal{S}(\heartsuit,\wp)\right].$$

which leads to

$$R(r) = r\theta(0) - \frac{1}{r^a} \mathcal{A} \left[ R\theta(\odot, \wp) + N\theta(\odot, \wp) + \mathfrak{g}(\odot, \wp) \right].$$

Using the initial condition (4), we obtain

$$R(r) = rh(\odot) - \frac{1}{r^{\mathfrak{a}}} \mathcal{A} \left[ R\theta(\odot, \wp) + N\theta(\odot, \wp) + \mathfrak{g}(\odot, \wp) \right],$$

thus, operating the inverse MT, we obtain

$$\mathcal{A}(\heartsuit,\wp) = G(\heartsuit,\wp) - \mathcal{A}\ell^{-1}\left[\frac{1}{r^a}\mathcal{A}\ell\left[R\theta(\heartsuit,\wp) + N\theta(\heartsuit,\wp)\right]\right],\tag{6}$$

which is called the recurrence relation of *ϑ*(, ℘), where

$$G(\mathbb{S}, \wp) = \mathcal{A}^{-1} \left[ rh(\mathbb{S}) + \mathcal{A} \left\{ \mathfrak{g}(\mathbb{S}, \wp) \right\} \right].$$

The approximate solution of Equation (3) can be expressed in terms of the power series

$$\mathfrak{d}(\mathbb{S}, \wp) = \sum\_{n=0}^{\infty} p^n \mathfrak{d}\_n(\mathbb{S}, \wp)\_{\prime} \tag{7}$$

and

$$N\theta(\heartsuit,\wp) = \sum\_{n=0}^{\infty} p^n H\_n \theta(\heartsuit,\wp),\tag{8}$$

.

where *p* ∈ [0, 1] is an embedding parameter and considered as a small parameter, whereas *ϑ*0(, ℘) is an initial guess of Equation (3). The following strategy can be operated to acquire He's polynomials as

$$H\_{\boldsymbol{n}}(\boldsymbol{\theta}\_{0} + \boldsymbol{\theta}\_{1} + \cdots + \boldsymbol{\theta}\_{n}) = \frac{1}{n!} \frac{\partial^{n}}{\partial p^{n}} \left( \boldsymbol{N} \left( \sum\_{i=0}^{\infty} p^{i} \boldsymbol{\theta}\_{i} \right) \right)\_{p=0} . \quad \boldsymbol{n} = 0, 1, 2, \cdots$$

With the help of Equations (7) and (8), we can obtain Equation (6) as

$$\sum\_{n=0}^{\infty} p^n \mathfrak{d}\_n(\mathbb{S}, \wp) = G(\mathbb{S}, \wp) - p.\mathcal{A}^{-1} \left[ \frac{1}{r^a} \mathcal{A} \left\{ R \left( \sum\_{n=0}^{\infty} p^n \mathfrak{d}\_n(\mathbb{S}, \wp) \right) + \sum\_{n=0}^{\infty} p^n H\_n \mathfrak{d}\_n(\mathbb{S}, \wp) \right\} \right]$$

Equating the similar components of *p*, we obtain

$$\begin{aligned} p^0: \mathfrak{d}\_0(\mathbb{S}, \wp) &= G(\mathbb{S}, \wp), \\ p^1: \mathfrak{d}\_1(\mathbb{S}, \wp) &= -\mathcal{A}^{-1} \left[ \frac{1}{r^a} \mathcal{A} \left\{ R\mathfrak{d}\_0(\mathbb{S}, \wp) + H\_0 \right\} \right], \\ p^2: \mathfrak{d}\_2(\mathbb{S}, \wp) &= -\mathcal{A}^{-1} \left[ \frac{1}{r^a} \mathcal{A} \left\{ R\mathfrak{d}\_1(\mathbb{S}, \wp) + H\_1 \right\} \right], \\ p^3: \mathfrak{d}\_3(\mathbb{S}, \wp) &= -\mathcal{A}^{-1} \left[ \frac{1}{r^a} \mathcal{A} \left\{ R\mathfrak{d}\_2(\mathbb{S}, \wp) + H\_2 \right\} \right], \\ \vdots & \vdots & \vdots \end{aligned} \tag{9}$$

Thus, we can generate Equation (7) in the collection of orders as

$$\theta(\heartsuit,\wp) = \theta\_0(\heartsuit,\wp) + p^1\theta\_1(\heartsuit,\wp) + p^2\theta\_2(\heartsuit,\wp) + + p^3\theta\_3(\heartsuit,\wp) + \cdots \tag{10}$$

Let *p* = 1; the analytical solution of Equation (3) is

$$\theta(\bigotimes\_{N} \wp) = \lim\_{N \to \infty} \sum\_{n=0}^{N} \theta\_{\mathfrak{n}}(\bigotimes \wp). \tag{11}$$

We put forward this strategy in the strength of upcoming mathematical applications.

**Theorem 1.** *Consider that and ζ are two Banach spaces with I* : → *ζ as nonlinear operator, such that ϑ*; *ϑ*<sup>∗</sup> ∈; , *I*(*ϑ*) − *I*(*ϑ*∗) ≤ *Kϑ* − *ϑ*∗, 0 < *K* < 1*. According to the Banach contraction theorem, I has a unique fixed point ϑ, i.e., Iϑ* = *ϑ. Let us recall Equation (11); we have*

$$\mathfrak{d}(\mathbb{S}, \wp) = \lim\_{N \to \infty} \sum\_{n=0}^{N} \mathfrak{d}\_{\mathfrak{n}}(\mathbb{S}, \wp)\_{\mathfrak{n}} \tag{12}$$

*and let us assume that* <sup>0</sup> = *ϑ*<sup>0</sup> ∈ S*p*(*ϑ*)*, where* S*p*(*ϑ*) = {*ϑ*<sup>∗</sup> ∈ : *ϑ* − *ϑ*∗ < *p*}*; then, we have*

$$\begin{aligned} (B\_1) \lhd\_n &\in \mathcal{S}\_p(\theta)\_\prime \\ (B\_2) \lim\_{n \to \infty} \lhd\_n &= \theta . \end{aligned}$$

**Proof.** (*B*1) In view of the mathematical induction for *n* = 1, we have

$$\|\|\mathbb{G}\_1 - \theta\_1\|\| = \|\|T(\mathbb{G}\_0 - T(\theta))\|\| \le K \|\|\theta\_0 - \theta\|\|.$$

Consider that the result is true for *n* = 1, so

$$||\mathbb{G}\_{n-1} - \theta|| \le K^{n-1} ||\theta\_0 - \theta||.$$

Thus, we have

$$\|\left(\mathbb{S}\_n - \theta\right)\| = \left\|\left|T(\mathbb{S}\_{n-1} - T(\theta))\right|\right\| \le K||\mathbb{S}\_{n-1} - \theta|| \le K^\eta ||\theta\_0 - \theta||.$$

Hence, using (*B*1), we have

$$||\heartsuit\_n - \theta|| \le K^n ||\theta\_0 - \theta|| \le K^n p < p\_\prime \,\_2\text{I}$$

where *p* is a contact point of a super norm *S*, which shows *<sup>n</sup>* ∈ S*p*(*ϑ*). *<sup>B</sup>*2: Since *<sup>n</sup>* <sup>−</sup> *<sup>ϑ</sup>* <sup>≤</sup> *<sup>K</sup>n*||*ϑ*<sup>0</sup> <sup>−</sup> *<sup>ϑ</sup>* and lim*n*→<sup>∞</sup> *<sup>K</sup><sup>n</sup>* <sup>=</sup> 0. Therefore, we have lim*n*→∞*ϑ<sup>n</sup>* − *ϑ* = 0 ⇒ lim*n*→<sup>∞</sup> *ϑ<sup>n</sup>* = *ϑ*.

#### **4. Numerical Examples**

In this segment, we deal with the MHPT to present the analytical and numerical solutions of time fractional shock wave equations and time fractional wave equations. The obtained results of these two problems show the performance and high accuracy of the suggested approach. The graphical results declare that this approach has good agreement.

#### *4.1. Example 1*

Consider the time fractional shock wave equation

$$D^{\mathfrak{a}}\_{\wp}\theta + \left(\frac{1}{c\_0} - \frac{\gamma + 1}{2}\frac{\theta}{c\_0^2}\right)D\_{\mathfrak{I}}\theta = 0, \qquad (\mathfrak{I}, \wp) \; \mathfrak{e} \; R \times [0, T], \qquad 0 < \mathfrak{a} \le 1,\tag{13}$$

where *c*<sup>0</sup> and *γ* are constants, and *γ* is the specific heat. If *c*<sup>0</sup> = 2, and *γ* = 1.5, the study case under consideration relates to the flow of air, as

$$
\frac{\partial^a \theta}{\partial \wp^a} + \left(\frac{1}{2} - \frac{5}{16} \theta\right) \frac{\partial \theta}{\partial \heartsuit} = 0,\tag{14}
$$

with the initial condition

$$
\theta(\text{\S}, \text{0}) = e^{-\frac{\text{\S}^2}{2}}.\tag{15}
$$

Taking the MT of Equation (14), we obtain

$$\mathcal{M}\left[\frac{\partial^{\alpha}\theta}{\partial\phi^{\alpha}} + \left(\frac{1}{2} - \frac{5}{16}\theta\right)\frac{\partial\theta}{\partial\zeta^{\zeta}}\right] = 0.$$

Using the definition of the MT, we can write it as

$$R(r) = r\theta(0) - \frac{1}{r^{\mathfrak{a}}} \mathcal{A} \left[ \left( \frac{1}{2} - \frac{5}{16} \mathfrak{d} \right) \frac{\partial \theta}{\partial \mathbb{S}} \right].$$

The inverse MT is

$$\mathcal{A}(\bigotimes , \wp) = \theta(\bigotimes , 0) - \mathcal{A}^{-1} \left[ \frac{1}{r^a} \mathcal{A} \left\{ \left( \frac{1}{2} - \frac{5}{16} \vartheta \right) \frac{\partial \vartheta}{\partial \heartsuit} \right\} \right] \rho$$

which is the recurrence relation of Equation (14); now, using Equation (7) together with the HPM, we obtain

$$\sum\_{n=0}^{\infty} p^n \theta\_n(\mathfrak{Z}, \wp) = \theta(\mathfrak{Z}, 0) - p. \mathcal{A}^{-1} \left[ \frac{1}{r^a} \mathcal{A} \ell \left\{ \frac{1}{2} \sum\_{n=0}^{\infty} p^n \frac{\partial \theta\_n}{\partial \mathbb{S}} - \frac{5}{16} \sum\_{n=0}^{\infty} p^n \theta\_n \frac{\partial \theta\_n}{\partial \mathbb{S}} \right\} \right], \tag{16}$$

by comparing, we can obtain the iterations

$$\begin{split} p^0: \mathfrak{G}\_0(\mathbb{G}, \wp) &= \mathfrak{G}(\mathbb{G}, 0), \\ p^1: \mathfrak{G}\_1(\mathbb{G}, \wp) &= -\mathcal{A}^{-1} \left[ \frac{1}{r^\alpha} \mathcal{A} \left\{ \frac{1}{2} \frac{\partial \theta\_0}{\partial \mathbb{S}} - \frac{5}{16} \theta\_0 \frac{\partial \theta\_0}{\partial \mathbb{S}} \right\} \right], \\ p^2: \mathfrak{G}\_2(\mathbb{G}, \wp) &= -\mathcal{A}^{-1} \left[ \frac{1}{r^\alpha} \mathcal{A} \left\{ \frac{1}{2} \frac{\partial \theta\_1}{\partial \mathbb{S}} - \frac{5}{16} \left( \theta\_0 \frac{\partial \theta\_1}{\partial \mathbb{S}} + \theta\_1 \frac{\partial \theta\_0}{\partial \mathbb{S}} \right) \right\} \right], \\ \vdots & \vdots \end{split}$$

which give the solutions

$$\begin{aligned} \theta\_0(\mathfrak{D}, \wp) &= e^{-\frac{\mathfrak{D}^2}{2}}, \\ \theta\_1(\mathfrak{D}, \wp) &= \left[ \frac{1}{2} \text{xe}^{-\frac{\mathfrak{D}^2}{2}} - \frac{5}{16} \text{xe}^{-\mathfrak{D}^2} \right] \frac{t^\mu}{\Gamma(a+1)}, \\ \theta\_2(\mathfrak{D}, \wp) &= \frac{1}{256} \left[ -25e^{-\frac{3\mathfrak{D}^2}{2}} + 80e^{-\mathfrak{D}^2} - 64e^{-\frac{\mathfrak{D}^2}{2}} + 753^2e^{-\frac{3\mathfrak{D}^2}{2}} - 1603^2e^{-\mathfrak{D}^2} - 643^2e^{-\frac{\mathfrak{D}^2}{2}} \right] \frac{t^{2\alpha}}{\Gamma(2a+1)}, \\ \vdots \end{aligned}$$

Proceeding with a similar process, the other elements of *ϑ<sup>n</sup>* can be calculated, and the series solutions are thus completely obtained. This series converges to the exact solution for high iterations. Finally, the analytical solution of *ϑ*(, *t*) can be obtained by using Equation (10), which is in full agreement with [6,13].

*4.2. Example 2*

Again, assume the time fractional wave equation

$$D^{\alpha}\_{\wp}\theta + \theta D\_{\heartsuit}\theta - D\_{\heartsuit\heartsuit\wp}\theta = 0,\tag{17}$$

with the initial condition

$$\theta(\text{\textquotedblleft},0\text{\textquotedblright}) = \text{3\textquotedblleft}\text{\textquotedblright}\left(\frac{\text{\textquotedblleft}-15}{2}\right)\_{\text{\textquotedblleft}}\tag{18}$$

According to the HPTM, the recurrence relation of Equation (17) can be written as

$$\theta(\mathcal{G}, \wp) = \theta(\mathcal{G}, 0) - \mathcal{M}^{-1} \left[ \frac{1}{r^{\mathfrak{a}}} \mathcal{M} \left\{ \vartheta \frac{\partial \vartheta}{\partial \mathbb{G}} - \frac{\partial}{\partial \wp} \left( \frac{\partial^2 \vartheta}{\partial \mathbb{G}^2} \right) \right\} \right],$$

Now, using Equation (7) together with the HPM, we obtain

$$\sum\_{n=0}^{\infty} p^n \theta\_{\mathbb{II}}(\mathfrak{Z}, \wp) = \theta(\mathfrak{Z}, 0) - p \mathcal{A} \ell^{-1} \left[ \frac{1}{r^a} \mathcal{A} \ell \left\{ \sum\_{n=0}^{\infty} p^n \theta\_n \frac{\partial \theta\_n}{\partial \mathbb{S}^1} - \frac{\partial}{\partial \wp} \left( \frac{\partial^2}{\partial \mathbb{S}^2} \sum\_{n=0}^{\infty} p^n \theta\_n \right) \right\} \right], \tag{19}$$

by comparing, we can obtain the iterations

$$\begin{split} p^{0} &= \theta\_{0}(\mathbb{S}, \wp) = \theta(\mathbb{S}, 0), \\ p^{1} &= \theta\_{1}(\mathbb{S}, \wp) = -\mathcal{M}^{-1} \left[ \frac{1}{r^{a}} \mathcal{M} \left\{ \theta\_{0} \frac{\partial \theta\_{0}}{\partial \mathbb{S}} - \frac{\partial}{\partial \wp} \left( \frac{\partial^{2} \theta\_{0}}{\partial \mathbb{S}^{2}} \right) \right\} \right], \\ p^{2} &= \theta\_{2}(\mathbb{S}, \wp) = -\mathcal{M}^{-1} \left[ \frac{1}{r^{a}} \mathcal{M} \left\{ \theta\_{0} \frac{\partial \theta\_{1}}{\partial \mathbb{S}} + \theta\_{1} \frac{\partial \theta\_{0}}{\partial \mathbb{S}} - \frac{\partial}{\partial \wp} \left( \frac{\partial^{2} \theta\_{1}}{\partial \mathbb{S}^{2}} \right) \right\} \right], \\ &\vdots \end{split}$$

which give the solutions

$$\begin{split} \theta\_{0}(3,\wp) &= 3\operatorname{sech}^{2}\left(\frac{3-15}{2}\right), \\ \theta\_{1}(3,\wp) &= 9\operatorname{sech}^{2}\left(\frac{3-15}{2}\right)\tanh\left(\frac{3-15}{2}\right)\frac{\wp^{a}}{\Gamma(1+a)}, \\ \theta\_{2}(3,\wp) &= \left[-\frac{27}{2}\operatorname{sech}^{8}\left(\frac{3-15}{2}\right) + 81\operatorname{sech}^{6}\left(\frac{3-15}{2}\right)\tanh^{2}\left(\frac{3-15}{2}\right)\right] \frac{\wp^{2a}}{\Gamma(1+2a)} \\ &- \left[\frac{63}{2}\operatorname{sech}^{6}\left(\frac{3-15}{2}\right)\tanh\left(\frac{3-15}{2}\right) - 36\operatorname{sech}^{4}\left(\frac{3-15}{2}\right)\tanh^{3}\left(\frac{3-15}{2}\right)\right] \frac{\wp^{2a-1}}{\Gamma(2a)}, \\ \vdots \end{split}$$

Proceeding with a similar process, the other elements of *ϑ<sup>n</sup>* can be calculated, and the series solutions are thus completely obtained. This series converges to the exact solution for high iterations. Finally, the analytical solution of *ϑ*(, ℘) can be obtained by using Equation (10) as

$$\theta(\heartsuit,0) = \mathfrak{Z}\operatorname{sech}^2\left(\frac{\heartsuit - 15 - \wp}{2}\right),\tag{20}$$

which is in full agreement with [6,13].

#### **5. Results and Discussion**

In this segment, we demonstrate the physical interpretations of the illustrated problems. We observe that the HPTM is fully capable of handling time fractional shock wave equations. Figure 1a–d show the surface solutions of *ϑ*(, ℘) for various time fractional equations in Brownian motion, and it is observed that *ϑ*(, ℘) reduces with the growth of and ℘ for *α* = 0.25, 0.50, 0.75 and 1. Figure 2a–d show the surface solutions of *ϑ*(, ℘) for the analytical solution obtained by the MHPT and the exact solution for various values of and ℘, respectively. It is observed that *ϑ*(, ℘) increases with the increase in and decreases with the increase in ℘ for *α* = 0.25, 0.50, 0.75 and 1.

**Figure 1.** The surface solutions of *u*(, ℘) with respect to and ℘ for distinct values of *α*. (**a**) Surface solution of *ϑ*(, ℘) when *α* = 0.25. (**b**) Surface solution of *ϑ*(, ℘) when *α* = 0.50. (**c**) Surface solution of *ϑ*(, ℘) when *α* = 0.75. (**d**) Surface solution of *ϑ*(, ℘) when *α* = 1.

**Figure 2.** The surface solutions of *ϑ*(, ℘) with respect to and ℘ for different values of *α*. (**a**) Surface solution of *ϑ*(, ℘) when *α* = 0.25. (**b**) Surface solution of *ϑ*(, ℘) when *α* = 0.50. (**c**) Surface solution of *ϑ*(, ℘) when *α* = 0.75. (**d**) Surface solution of *ϑ*(, ℘) when *α* = 1.

#### **6. Conclusions**

In this paper, we successfully apply the HPTM to achieve the approximate and analytical solutions of nonlinear time fractional shock wave and wave equations. This study demonstrates the importance of fractional derivatives and the technique of dealing with the recurrence relation. Since the MT is limited to linear problems only, whereas the HPM is applicable to nonlinear problems, we conclude that the MHPT is the best tool to provide significant results for both linear and nonlinear problems. The MHPT is here directly applied to obtain the series solutions. The present scheme shows higher efficiency and fewer computations than other approaches studied in the literature. All the iterations were calculated with the help of MAPLE Software. The solution graphs show that this approach is suitable for a broad variety of nonlinear fractional differential equations in science and engineering. In future work, this approach could further be extended to solve various nonlinear obstacle problems.

**Author Contributions:** Conceptualization, formal analysis and funding acquisition, J.F.; investigation, methodology, software and writing—original draft, M.N.; validation and visualization, M.H.; supervision and writing—review and editing, A.A. All authors have read and agreed to the submitted the manuscript.

**Funding:** This work was supported by the Foundation of Yibin University, China (grant No. 2019QD07).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

