1. Introduction
Nonlinear partial differential equations serve as fundamental tools for describing a diverse array of phenomena across numerous scientific disciplines, including physics, biology, chemistry, economics, epidemiology, engineering, and beyond. Despite their ubiquity, exact solutions to the majority of these equations remain elusive. Consequently, various computational methodologies have been devised by researchers to address this challenge. Notable among these techniques are the Adomian decomposition method pioneered by G. Adomian [
1,
2,
3,
4], the variational iteration method introduced by J. He [
5,
6], the finite difference method developed by Crank–Nicolson in [
7], the homotopy analysis method (HAM) proposed by Liao [
8,
9,
10,
11,
12], the homotopy perturbation method advocated by J. He [
13,
14,
15], and the differential transform method utilized by several researchers [
16,
17].
In recent years, a number of researchers have focused on investigating solutions to these partial differential equations through the integration of diverse techniques with the Laplace transform. Notable examples include the Laplace decomposition method detailed in references [
18,
19], the homotopy perturbation transform method explored in [
20], and the introduction of the homotopy analysis transform method (HATM) as documented in [
21]. Since then, the HATM has been widely implemented by many authors to solve various types of mathematical equations. For example, it was employed in [
22] to solve a fractional equation modeling gas dynamics, and it was utilized in [
23] to solve a fractional Fokker–Planck equation. In [
24], it is applied to solve an integro-differential equation; in [
25], the HATM was used to obtain exact solutions for fractional equations of type Korteweg–de Vries and Korteweg–de Vries–Burger’s; in [
26], it was implemented to solve fractional Navier–Stokes equations; in [
27], a modified version of the HATM was utilized to numerically solve a fractional-order smoking epidemic model; in [
28], the HATM was employed to solve a nonlinear integro-differential equation; etc.
The goal of this study was to derive a numerical solution for a singular nonlinear second-order hyperbolic pseudo-differential equation based on the HATM. Moreover, the effectiveness of the proposed technique was explored through a set of illustrative example as well as a comparative analysis with existing methodologies in the literature. In particular, we employed the HATM to numerically solve the following pseudo-hyperbolic initial value problem:
subject to the constraints
where
,
,
, and
are constants, and
and
are known functions.
In fact, at
, Equation (
1) reduces to a singular one-dimensional wave equation, and at
, it reduces to a singular one-dimensional pseudo-wave equation.
The nonlinear singular one-dimensional pseudo-hyperbolic equation is a type of PDE that involves mixed partial derivatives in space and time, and it exhibits both hyperbolic and singular characteristics. In general, these equations are often used to model complex physical phenomena. In our model, the equation describes how waves propagate, attenuate, and interact in a medium with radial symmetry, such as sound waves in a cylindrical duct or stress waves in a cylindrical rod. The nonlinear terms in the equation represent complex interactions, like nonlinear elastic responses or external forces acting on the medium. These nonlinear terms can lead to complex behaviors such as wave steepening, shock formation, or solitons. They often represent interactions within the medium or external forces that depend on the wave’s state. More precisely, the equation models a radial wave propagation in an elastic medium, where the Bessel operator reflects the radial nature of wave propagation.
To achieve our goal, we need the Laplace transform for a second-order partial derivative of the function
with respect to
t which is given by
Let us mention that linear and nonlinear versions of Equation (
1) were considered in [
29], where a combination of the double Laplace transform and the Adomian decomposition method was applied to solve these versions.
The rest of this paper is organized as follows. In
Section 2, we provide the basic outlines of the HATM. In
Section 3, the HATM is employed to derive the numerical method for solving model (
1).
Section 4 is devoted to the application of the derived iterative scheme to a set of test examples to evaluate its efficiency. Finally, in
Section 5, we comment on our findings and conclusions.
2. Basic Outline of HATM
The homotopy analysis transform method was first introduced in [
21]. It is an elegant combination of the Laplace transform method and the homotopy analysis method. Its implementation does not require discretization or any other restrictive assumptions, which avoid round-off errors that may affect the reach of a closed-form solution of the problem. Moreover, it does not require decomposing the nonlinear terms in the equation under consideration, which usually leads to tedious algebraic work. Additionally, for some strongly nonlinear problems other techniques may fail to produce a convergent series solution, while the existence of the parameter
ℏ in theh HATM can be used to adjust and control the convergence region to ensure the convergence of a series solution for such problems. This reveals that the HATM is a very effective and accurate method in handling the solution of various nonlinear types of partial differential equations. Thus, our motivation for this work was to determine a numerical solution for the singular problem (
1) based on the HATM that is easy to implement and does not require any restrictive assumptions needed for the implementation of the aforementioned methods. To formulate the basic outlines of the HATM, we consider a general partial differential equation as follows:
subject to the initial conditions
where
represents a linear differential operator,
represents a nonlinear differential operator, and
, and
are known functions.
Taking the Laplace transform
L for both sides of Equation (
3) implies that
that is,
Thus, according to the theory of the HAM [
8], we can define an operator
as follows:
where
and
is a real-valued function in
, and
q.
Then, we can define the zeroth-order deformation equation as follows:
where
is an embedding parameter;
L is the Laplace transform operator;
is an initial approximation of the exact solution
;
is an unknown function of
and
q; and
ℏ is a nonzero auxiliary parameter, which plays a fundamental role in the convergence of the series solution.
It follows from Equation (
4) that
at
and
at
. Thus, as
q moves continuously from 0 to 1, the function
deforms from the initial guess
toward the exact solution
.
Next, the Taylor series expansion of
in powers of
q yields
where
As Liao pointed out in [
11], if the initial guess
and the parameter
ℏ are chosen properly, then power series (
5) will converge at
to a solution of the considered partial differential equation, and this solution is expressed in series form as follows:
In fact, the success of the convergence of the series solution in the HATM depends on three elements, namely, the choice of the operator to be inverted, the value of the parameter ℏ, and the initial guess. The choice of the suitable operator takes into account that the operator is easily inverted. The existence of the parameter ℏ distinguishes the HATM from other semi-analytical methods, as its occurrence in the series solution enables one to adjust the convergence interval as required. But in the absence of easy theoretical techniques that can be implemented to determine its value, an alternative graphical technique based on the h-curve resulting from the plots of the approximate truncated series solutions is used to determine a range of values for this parameter, usually the interval on which the curve is almost horizontal. Finally, the choice of the initial guess depends on the operator to be inverted. It can be chosen from the initial or boundary conditions in the problem according to the inverted operator, and it should satisfy the differential equation under consideration.
Now, differentiating Equation (
4)
k times with respect to
q, setting
, and dividing the result by
produce the following kth-order deformation equation:
where
and
Finally, taking the inverse Laplace transform for both sides of Equation (
6) implies that the components
can be determined by applying the following iterative formula:
where
Let us mention that the accuracy of the numerical solution will be improved as the number of terms in the truncated series solution is increased. On the other hand, the computational process can be terminated when either of the computed components lead to the closed-form solution or the terminated series solution provide the required accuracy of the numerical solution.
4. Numerical Results
In this section, we test the efficiency and accuracy of scheme (
7) in handling the numerical solutions for problems of type (
1). It is notable to mention that all computations in this section were conducted using the software Mathematica 8.04.0., and all the numerical results in the tables were computed using the default machine precision and accuracy of 17 significant digits rounded to 6 significant digits. Now, we apply this scheme to the following set of examples.
Example 1. Consider the singular linear partial differential equationsubject to the following initial constraints:where . Solution.
Let
. Then, in view of (
7), we obtain
⋯
Thus, the series solution is given by
It is noticed that at
, the components of the series solution (
8) are reduced to
Continuing in this manner, we obtain
for all
. Hence, the series solution (
8) reduces to
which is the exact solution for this equation given in [
29].
Figure 1 shows the
h-curve corresponding to the 13th-order series solution. This figure shows that the values of
ℏ that lead to the convergence of the series solution (
9) should be selected from the range
.
Figure 2 shows plots of the truncated series solutions for Example 1 with a different number of terms at
and
. This figure shows the rapid convergence of these approximate solutions toward the exact solution. It illustrates that the curves of the truncated series solutions are getting closer to the exact solution as the number of terms in the truncated series solution is increased.
Table 1 and
Table 2 show the absolute error
resulting from taking the difference between the exact solution and the approximate solution of several orders of Example 1 at different values of
t, and
x. These data illustrate the rapid convergence of the approximate solution just after a few terms. Moreover, it appears from an error analysis of these tables that the convergence of the approximate solutions is more rapid at
than at the other values of
ℏ in the interval of convergence.
Example 2. Consider the singular nonlinear partial differential equationsubject to the following initial conditions: Solution.
Let
. Then, in view of (
7), we obtain
Continuing in this manner, we obtain
for all
. Thus, the series solution is given by
which is the exact solution for this equation presented in [
29], without using Adomain polynomials.
Example 3. Consider the singular nonlinear partial differential equationsubject to the following initial conditions: Solution.
Let
. Then, in view of (
7), we obtain
Thus, the series solution is given by
In fact, at
, the components of the series solution (
9) are reduced to
Hence, series solution (
9) reduces to
Figure 3 shows the
h-curve corresponding to the 10th-order series solution. This figure shows that the values of
ℏ that lead to the convergence of series solution (
9) should be selected from the interval
.
Figure 4 shows plots of the truncated series solutions for Example 3 with a different number of terms at
and
. This figure shows the rapid convergence of these approximate solutions toward the exact solution. It illustrates that the curves of the truncated series solutions are getting closer to the exact solution as the number of terms in the truncated series solution is increased.
The exact solution of Example 3 is
.
Table 3 and
Table 4 show the absolute error
resulting from taking the difference between the exact solution and the approximate solution of several orders at different values of
t, and
x. These data demonstrate the rapid convergence of the approximate solution just after a few terms. Moreover, it appears from an error analysis of these tables that the convergence of the approximate solutions is more rapid at
than at the other values of the parameter
ℏ in the interval of convergence.
It is seen from the tables and figures in the above examples that the error is decreased and the accuracy of the approximate solution is refined as the number of iterations is increased, which reflects the reliability of the approximate solutions obtained using the HATM. On the other hand, these examples illustrate the stability of the derived iterative scheme (
7) numerically and graphically. In both cases, this scheme generates iterations used to construct the approximate series solutions. The plots in
Figure 2 and
Figure 4 show the behavior of the first few terms in the sequence of partial sums of these series solutions,
. These plots display that the graph of
is getting closer and closer to the graph of the exact solution as
n increases, which means that the error in each step is decreased and the global error is bounded, i.e., the scheme is stable, and the sequence of partial sums converges to the analytical solution. In addition to this, the numerical results in
Table 1,
Table 2,
Table 3 and
Table 4 show that the absolute error in the computed values of these partial sums at selected values of
x and
t decreases rapidly as
n increases. This occurs for distinct values of the parameter
ℏ in the interval of convergence. Again, this means that the scheme is stable; otherwise, the global error increases without a bound as
n increases.
5. Conclusions
In this study, an iterative scheme was obtained to find a numerical solution for a singular nonlinear second-order hyperbolic pseudo-differential equation based on the HATM. A series of three examples was given to demonstrate the accuracy and efficiency of the obtained iterative scheme. The error analysis conducted for the computed numerical solutions of these examples was discussed and illustrated both graphically and numerically. The plots in
Figure 2 and
Figure 4 show the behavior of the first few truncated series solutions. These plots are getting closer and closer to the graph of the exact solution as the number of terms in these series increases, which means that the error is reduced in each step, i.e., the derived scheme is stable, and the approximate series converges to the analytical solution. Also, the results in
Table 1,
Table 2,
Table 3 and
Table 4 show that the computed absolute error is inversely proportional to the number of terms in the series solution, which shows that the numerical solutions obtained using the HATM are accurate and reliable. In addition, it is noticed from these tables that the best performance of this method is obtained at
. On the other hand, the numerical solutions obtained using this method were compared with those obtained using the double Laplace transform given in [
29]. It was found that the results obtained using the two methods coincide. But the implementation of the HATM avoids the use of Adomian polynomials required for the application of the other method, as these polynomials require tedious computational work.
Thus, we conclude that the HATM is an efficient method and can be implemented easily for solving this type of problems without any restrictive assumptions. Motivated by these results, we look forward to considering a numerical solution for a two-dimensional singular nonlinear hyperbolic equation based on the HATM in our future work. Finally, we would like to mention that the authors contributed equally to this work.