A Reliable Combination of Double Laplace Transform and Homotopy Analysis Method for Solving a Singular Nonlocal Problem with Bessel Operator

: In this article, we present a numerical iterative scheme for solving a non-local singular initial-boundary value problem by combining two well known efﬁcient methods. Namely, the homotopy analysis method and the double Laplace transform method. The resulting scheme is tested on a set of test examples to illustrate its efﬁciency, it generates the exact analytical solution for each one of these examples. The convergence of the resulting numerical solutions of these examples is tested both graphically and numerically

Unfortunately, theoretical analytical solutions of this type of problems mostly can not be determined, especially for nonlinear problems.Thus, in the literature various numerical techniques have been developed by several researchers to determine approximate solutions for these problems.In [10], finite difference methods have been used to address a one-dimensional parabolic partial differential equation with initial and nonlocal boundary conditions that involve nonlocal integral terms.Meanwhile, a combination of finite difference and orthogonal function approximation techniques has been utilized to solve a one-dimensional parabolic equation featuring two integral conditions, as demonstrated in [11].Furthermore, a hybrid approach integrating finite difference and spectral methods has been proposed for solving a one-dimensional wave equation accompanied with an integral condition, as presented in [12].Moreover, the Chebyshev spectral method is used to solve a class of local and nonlocal elliptic boundary value problems in [13].Also, the finite element method is employed to solve a nonlocal problem of Kirchhoff type in [14].
On the other hand, several computational techniques that do not require discretization have been developed in the last decades.For example the differential transform method is applied in [15,16].The Adomian decomposition method given in [17][18][19], which is developed to treat nonlinearity in partial differential equations.In addition to that Laplace decomposition method is given in [20,21].One more analytical method is the variational iteration method given in [22,23].Another method is the homotopy perturbation presented in [24][25][26].A powerful analytical technique which is widely used by many researchers is the homotopy analysis method introduced by Liao [27][28][29][30][31].
In this work, we combined two well known techniques, the homotopy analysis method and the double Laplace transform method, to develop a numerical scheme, named as homotopy analysis double Laplace transform method (HADLTM), to solve a singular onedimensional parabolic equation subject to Dirichlet conditions and a non-local condition of integral type.To our knowledge it is the first time in the literature where such combination between these two methods is applied.
Here we will consider the following problem: where L is the operator ∂ ∂t − 1 x ∂ ∂x x ∂ ∂x , and and w(x) are known functions.Problem (1) can be written in an operator form as: Kη = (L, )η =: ℱ , in which K is an unbounded operator that satisfies K : A → H, with domain all functions of η in the set: D((Ω) = {η : η, η t , η x , and η xx ∈ L 2 x (Ω)} which satisfy the given boundary conditions, and A is a Banach space of functions η equipped with a finite norm as [32]: with values in the Hilbert space H consisting of the functions ℱ = ( , w) ∈ L 2 x (Ω) × W 1 x (Ω) equipped with the norm [32]: The double Laplace transform of a function f (x, t) is given as [33][34][35]: where x and t are positive real variables, and r and s are complex variables.On the other hand, the double Laplace transform of its time-partial derivative; , is given by: where F is the double Laplace transform of f .The rest of the article is organized as follows: In Section 2, we recall some existence results of the solution of problem (1).In Section 3, we present the development of the numerical scheme for solving problem (1) based on the HADLTM.In Section 4, we provide several examples to test the applicability and efficiency of the developed scheme.Finally, we present some comments and conclusions in Section 5.

Existence and Uniqueness of the Solution
Here, we present some existence and uniqueness results of the solution of problem (1).First, we recall the two sided a prior estimates: Theorem 1 ([32,36]).For every function η ∈ D(K), we have: Theorem 2 ([32,36]).For every function η ∈ D(K), the following energy inequality holds true: for some positive constant C independent of η.
As pointed in it [32], in view of Theorem 1, it follows that the operator K : A → H is continuous.Thus, Theorem 2 implies that the set R(K) ⊂ H; the range of K, is closed.Hence, the inverse operator K −1 exists and is continuous.To verify the existence of the solution of (1), we need to verify that the set Im(K); the image of K, coincides with the whole Hilbert space H. Theorem 3 ([32,36]).For any two functions ∈ L 2 x (Ω) and w ∈ W 1 x (Ω), there exists a unique solution η = K −1 Ψ of problem (1) that satisfies the inequality where Ψ = ( , w), and C is a positive constant which does not depend on η.

Method Development
Consider a general partial differential equation as: where θ is an unknown function in x and t, R is a linear differential operator, Ñ represents a nonlinear differential operator, and f is a given function.Applying double Laplace transform on both sides of (3), implies: or equivalently, where Θ(r, s) is the double Laplace transform of θ(x, t).According to the homotopy analysis method [27], we define the operator: where p ∈ [0, 1], and φ is a real valued function of x, t and p.Thus, the zeroth-order deformation equation will be on the form: where h = 0 is an auxiliary parameter, p ∈ [0, 1] is an embedding parameter, θ 0 (x, t) is an initial guess to start with to get the solution θ(x, t), and φ(x, t; p) is an unknown function.
In view of Equation (4), it is clear that for p = 0 and p = 1 we obtain: Thus, as p increases from 0 to 1, the function φ(x, t; p) deforms from θ 0 (x, t) to the exact solution θ(x, t).
Then, using the Taylor series expansion of φ(x, t; p) with respect to p gives: where As Liao mentioned in [30] if the auxiliary parameter h, the inverted operator L, and the initial guess θ 0 (x, t) are chosen properly, then the power series (5) will converge at p = 1 to one of the solutions of the original equation, and this solution is given in a series form as: as pointed out by Liao in [28], the parameter h helps one in controlling and adjusting the convergence region of the series solution.The values of this parameter can be determined through the h−curve.Now, differentiating Equation (4) k-times with respect to p, dividing by k!, then setting p = 0, produces the kth order deformation equations as: where and Then, applying the inverse double Laplace transform to both sides of ( 6), the components θ k (x, t) of the HADLTM can be determined recursively through the formula: where

Application of the Method
The presence of the integral condition in problem (1) complicates the computations.Therefore, to avoid this difficulty in exploring the applicability and efficiency of the HADLTM for solving this problem, we consider the forthcoming equivalent problem (9), in which the integral condition is replaced by classical conditions.
In particular 1 0 x θ(x, 0) dx = c = 0, as the initial condition θ(x, 0) = η(x) satisfies the compatibility condition, which implies provided that 1 0 x f (x, t)dx = 0, where f , η and ϕ are given functions.Now, in view of (2), taking double Laplace transform for both sides of (9), we obtain: Thus, we define the operator N [φ(x, t; q)] as: Hence, the kth order deformation equation is given by: Hence, the series solution is given as: where To test the efficiency of the scheme (10) in handling the numerical solutions of problems of the type (9), it is applied to the forthcoming set of test examples.
Example 1.Consider the homogeneous PDE: satisfying the following initial and boundary conditions: In view of (10) we obtain: In general we obtain: Now, if we choose the parameter h within the range −2 < h < 0, then the series (11) converges to the exact solution given by: Figure 1 displays the h-curve corresponding to the 15th order truncated series solution at x = 0.75, from which it appears that the values of the parameter h required for the convergence of the series solution are located in the interval −1.8 < h < −0.2.
Figure 2 shows the plots of the truncated series solution using different number of terms together with the corresponding exact solution of Example 1 at x = 0.7 and h = −0.6.It shows that the truncated series solution of order j = 5 is very close to the exact solution, which shows the rapid convergence of the proposed method.Table 1 shows the absolute error Er [j] = |θ(x, t) − θ [j] (x, t)| due to the difference between the exact solution and the approximate solution of distinct orders of Example 1 at different values of x and t.It illustrates the rapid convergence of the approximate solution obtained by using this method.Table 1.Absolute error at h = −0.9 and different values of j, x, and t. t 0.1 0.5 1 5 x j Er [j]  Er [j]  Er [j]  Er [j]   0. x j Er [j]  Er [j]  Er [j]  Er [j]   0.
subject to the conditions:

Solution.
Let θ 0 (x, t) = θ(x, 0) = 1 + x 4 4 − ln(x).Then, in view of (10) we get: Proceeding on this manner we obtain: Hence, if the parameter h satisfies −2 < h < 0, then the series ( 12) converges to the exact solution given by: Figure 3 shows the h-curve corresponding to the 16th order truncated series solution at x = 0.75, it appears that the valid values of the parameter h that lead to a convergent series solution are located in the interval −1.8 < h < −0.2. Figure 4 shows plots of the truncated series solution using different number of terms together with the corresponding exact solution of Example 2 at x = 0.8 and h = −0.7.It shows that the truncated series solution of order j = 5 almost coincide with the exact solution, which shows the rapid convergence of the proposed method.Table 2 shows the absolute error Er [j] = |θ(x, t) − θ [j] (x, t)| due to the difference between the exact solution and the approximate solution of distinct orders of Example 2 at different values of x and t.It illustrates the rapid convergence of the approximate solutions generated by this method.
subject to the constraints: (10), we obtain: In general we get: hence, the series solution takes the form: Thus, if the parameter h satisfies −2 < h < 0, then the series ( 12) converges to the exact solution: Figure 5 shows the h-curve corresponding to the 15th order truncated series solution at x = 0.7, from which it appears that the valid values of the parameter h that lead to a convergent series solution are located in the interval −1.8 < h < −0.2. Figure 6 shows the plots of the truncated series solution using different number of terms together with the corresponding exact solution of Example 3 at x = 0.3 and h = −0.7.It shows that the truncated series solution of order j = 5 is very close to the exact solution, which demonstrates the rapid convergence of the proposed method.
Table 3 presents the absolute error Er [j] = |θ(x, t) − θ [j] (x, t)| due to the difference between the exact solution and the approximate solution of distinct orders of Example 3 at different values of x and t.It illustrates the rapid convergence of the approximate solutions generated by this method.x j Er [j]  Er [j]  Er [j]  Er [j]   0.
subject to the constraints: , in view of (10), we have: Proceeding on this manner we get: hence, the series solution is given as: It follows that if the parameter h lies in the range −2 < h < 0, then the series (14) converges to the analytical solution: θ(x, t) = θ 0 (x, t) + sin(t) = sin(t) + x 2 2 − ln(x), Figure 7 shows the h-curve corresponding to the 15th order truncated series solution at x = 0.7, it appears that the valid values of the parameter h that lead to a convergent series solution are located in the interval −1.8 < h < −0.2. Figure 8 shows the plots of the truncated series solution using different number of terms together with the corresponding analytical exact solution of Example 4 at x = 0.7 and h = −0.7.It shows that the truncated series solution of order j = 5 is very close to the exact solution, which indicates the rapid convergence of the proposed method.Table 4 presents the absolute error Er [j] = |θ(x, t) − θ [j] (x, t)| due to the difference between the exact solution and the approximate solution of various orders of Example 4 at different values of x and t.It illustrates the rapid convergence of the approximate solution obtained by using this method.x j Er [j]  Er [j]  Er [j]  Er [j]   0.

Conclusions
In this article, a numerical scheme is developed to solve a mathematical model subject to a non-local condition of integral type by combining two well known methods; the double Laplace transform and the homotopy analysis method.The derived scheme is applied on a set of test examples to demonstrate it's reliability and efficiency.This scheme generates the exact solution for each one of these example.The convergence of the obtained approximate solutions is tested graphically as shown in Figures 2, 4, 6, and 8.These figures show rapid convergence of the numerical solutions towards the exact solution just after few iterations.

Figure 1 .
Figure 1.The h-curve corresponding to the 15th order series solution at t = 0.

1 ΘFigure 2 .
Figure 2. Approximate solution θ [j] (x, t) with various values of j together with the exact solution θ(x, t) of Example 1.

Figure 3 .
Figure 3.The h-curve corresponding to the 16th order approximate solution at t = 0.

1 ΘFigure 4 .
Figure 4. Approximate solution θ [j] (x, t) with various values of j together with the exact solution θ(x, t) of Example 2.

Figure 5 .
Figure 5.The h-curve corresponding to the 15th order series solution at t = 0.

1 ΘFigure 6 .
Figure 6.Approximate solution θ [j] (x, t) with various values of j together with the exact solution θ(x, t) of Example 3.

Figure 7 .
Figure 7.The h-curve for the 15th order series solution at t = 0.

1 ΘFigure 8 .
Figure 8. Approximate solution θ [j] (x, t) with various values of j together with the exact solution θ(x, t) of Example 3.

Table 3 .
Absolute error at h = −1.1 and different values of j, x, and t.

Table 4 .
Absolute error at h = −1.1 and different values of j, x, and t.