Approximating Solutions of Non-Linear Troesch’s Problem via an Efﬁcient Quasi-Linearization Bessel Approach

: Two collocation-based methods utilizing the novel Bessel polynomials (with positive coefﬁcients) are developed for solving the non-linear Troesch’s problem. In the ﬁrst approach, by expressing the unknown solution and its second derivative in terms of the Bessel matrix form along with some collocation points, the governing equation transforms into a non-linear algebraic matrix equation. In the second approach, the technique of quasi-linearization is ﬁrst employed to linearize the model problem and, then, the ﬁrst collocation method is applied to the sequence of linearized equations iteratively. In the latter approach, we require to solve a linear algebraic matrix equation in each iteration. Moreover, the error analysis of the Bessel series solution is established. In the end, numerical simulations and computational results are provided to illustrate the utility and applicability of the presented collocation approaches. Numerical comparisons with some existing available methods are performed to validate our results.


Introduction
Mathematical modeling through differential equations is of great importance in different branches of science and engineering with their ability to provide more realistic simulations to real life phenomenons. The non-linear Troesch's model problem has been used in the investigation of the theory of gas porous electrodes [1,2]. This model, as a two-point boundary value problem (BVP) has also found application in the confinement of a plasma column by radiation pressure [3]. To be more precise, the following two-point non-linear BVP is known as Troesch's problem [4,5] ψ (τ) = α sinh(αψ(τ)), (1) with the following boundary conditions Here, the constant parameter α is positive. In the closed-form, the exact solution of (1) has been obtained with the help of the elliptic function Sc(·|·) of Jacobian type [6] in the form defined as Sc(α|x) = sin β cos β , where β and α are related through the integral representation α = β 0 dt 1−x sin 2 t . For α > 1, Troesch's problem is known to be an unstable two-point BVP due to existence of a pole, which is located approximately at τ s = ln(8/ψ (0))/α [4]. These singularities hinder most of the numerical methods to capture the solutions of the Troesch's problem, especially for large values of the parameter α.
Over the past decades, a considerable attention has been drawn to the numerical solutions of (1). Among other existing approaches, let us mention the following methods: the shooting schemes [5][6][7][8][9], the embedding technique [10], the combined approach based on modified decomposition and Laplace transform [11], the modified homotopy perturbation algorithm [12], the finite element method based on B-splines [13], the sinc-Galerkin methods [14,15], the discontinuous Galerkin finite element scheme [16], the stochastic algorithm based on artificial neural network [17], the shifted Jacobi-Gauss collocation scheme [18], the finite difference scheme based on a Shishkin mesh [19], the shifted Bessel Tau method [20], the iterative method based on Green's function and optimal homotopy analysis technique [21], the one-step hybrid block method [22], and the wavelet-homotopy analysis method [23].
In this manuscript, an efficient and reliable method based on the quasi-linearization technique combined with Bessel functions is proposed to obtain an accurate approximate series solution of Troesch's model problem (1). Additionally, the direct Bessel collocation approach is also developed for this model problem. In the former method, we solve the original model via converting it into a sequence of linearized subproblems in conjunction with an iteration parameter r, which will be set at most 5. In the latter approach, we solve it directly via the Bessel matrix algorithm, which converts the governing differential equation into a system of non-linear algebraic equations. The Bessel polynomials first appeared in the study of the solution of classical wave equation and were systematically introduced in [24]. However, they are recently adapted to solve differential equations through a collocation procedure [25][26][27][28]. The main feature of the novel Bessel polynomials is that all coefficients are positive integers. It is also worth mentioning that the Bessel functions used here are different from the Bessel functions of the first kind, which were previously considered in the literature, see cf. [29][30][31][32]. Collocation-based approaches along with (orthogonal) basis functions are very common in solving differential and integral equations due to their simplicity and efficiency. Among many different bases used in the literature, let us consider the methods based on the Legendre, Chebyshev, Jacobi, Genocchi, and Hermit functions reported in [33][34][35][36][37][38].
The rest of this paper is organized as follows: in Section 2, a brief discussion on the Bessel polynomials is given. Hence, the technique of quasi-linearization for the model problem is discussed. Section 3 is devoted to the methodology of the direct Bessel-collocation procedure for the Troesch model problem. Converting the boundary conditions into a matrix form is also done afterwards. The second and efficient Bessel collocation approach based on the quasi-linearlization technique is given in Section 4. An error analysis of the Bessel solution is given in Section 5. Numerical experiments and discussions of the presented methods are reported in Section 6. Finally, a conclusion is made in the last section.

Bessel Functions
Let us assume that k is a positive integer. The Bessel polynomials are known to be the solutions of the following second-order differential equation which has naturally appeared in the study of classical wave equation [24], see also [25]. The first two Bessel polynomials are B 0 (τ) = 1 and B 1 (τ) = 1 + τ. The remaining Bessel polynomials can be determined via the following recurrence formula In an explicit form, we can write these polynomials as Obviously, all coefficients of each B k (τ) are positive and the constant term for each one is 1. It can be proved that the set of these polynomials forms an orthogonal system on the unit circle. The corresponding weighting function is ω(τ) = e − 2 τ . The key of the presented method is that the unknown solution ψ(τ) and its secondorder derivative in (1) are expressed in terms of a truncated series of Bessel polynomials with expansion coefficients as Hence, our aim will be reduced to find the expansion coefficients d k for k = 0, 1, . . . , M. Let us introduce the vector of unknown as Thus, we may write the relation (7) in the following representation form being Ω Ω Ω M (τ) the vector of Bessel basis defined as It can be easily shown that Ω Ω Ω M (τ) in (8) is expressed as a product of the monomial basis functions with a lower-triangular matrix as where Π Π Π M (τ) = 1 τ τ 2 . . . τ M and C C C can be represented as . By utilizing (8) and (9), we may represent the approximated solution ψ M (τ) in (7) as follows

Quasi-Linearization Approach
The technique of quasi-lineraization aims to transform the original model problem (1) into a sequence of linear equations. Afterwards, we will apply the Bessel collocation procedure to the resultant linearized models. In this respect, the basic ideas underlying the solutions of our model problems via quasi-linearization method (QLM) are briefly described, see [39][40][41] for more detailed discussions and applications.

Direct Approach
To proceed, we partition the domain [0, 1] uniformly into M subdomains. Below, the following collocation points are used We next substitute the collocation points (14) into (10), we arrive at Our next aim is to represent the second-order derivatives of ψ(τ) at the collocation points in a matrix representation form. For this purpose, we need to calculate d 2 dτ 2 Π Π Π M (τ). It can be easily checked that where the matrix E E E denotes the operational matrix of differentiation and is defined as .
By applying differentiation on (16), we obtain Lemma 1. The matrix representation form of d 2 dτ 2 ψ M (τ) at the collocation points (14) is given by where the matrix E E E and the vector ∆ ∆ ∆ are defined in (16) and (15), respectively, and Proof. We differentiate (10) two times with respect to τ to obtain We combine relations (17) and (19) to obtain Now, it is sufficient to insert the collocation points (14) into (20). By considering the definitions of ∆ ∆ ∆ and Ψ Ψ Ψ (2) the proof is concluded.
We now proceed by placing the collocation points (14) into the model problem to obtain In the matrix notation form, we are able to write the preceding equations as where the matrix A A A and the right-hand-side vector Z Z Z are Suppose that the approximate solution of (1) can be written as (7) or (10). By placing the relations (15) and (18) into (22), we obtain the following fundamental matrix equation It can be clearly seen that the relation (23) is a non-linear matrix equation, which can be solved for the vector of unknowns D D D M as the Bessel coefficients by using, e.g., Newton-type solvers. However, we are left with the task of entering the boundary conditions (2) into the fundamental matrix equation.

Boundary Conditions in the Matrix Forms
To obtain the approximated solution of our model problem (1) by solving the fundamental matrix Equation (23), we need to take into account the boundary conditions (2). First, we consider the initial condition ψ(0) = 0. For this purpose, we let τ → 0 in (10) to obtain For the end condition ψ(1) = 1, we tend τ → 1 in (10) to arrive at

Remark 1.
We note that the computational complexity of evaluating d 2 dτ 2 Π Π Π M (τ) at the collocation points via (17) is on the order O((M + 1) 2 ). In the following, we propose an alternative to reduce the involved complexity within O(M + 1). The following algorithm will be invoked to calculate the second-order derivative of the basis functions Π Π Π M (τ) directly. Indeed, Algorithm 1 takes Π Π Π M (τ) as input and the output is Π Π Π

Algorithm 1 The computation of
We are now able to enter the boundary conditions (2) into the fundamental matrix Equation (23). In this respect, we replace the first row as well as the last row of the augmented matrix [Σ Σ Σ; Z Z Z] in (23) by the row matrices [ Σ Σ Σ 0 ; 0] and [ Σ Σ Σ 1 ; 1]. Therefore, the resultant modified fundamental matrix becomes [ Σ Σ Σ; Z Z Z]. After solving this modified matrix equation, the unknown coefficients d k , k = 0, 1, . . . , M will be computed and, therefore, the approximate solution ψ M (τ) for the model problem in (1) will be determined.

QLM-Bessel
Our goal is to solve the two-point BVP (1)- (2) approximately such that the desired solution is represented in terms of the truncated Bessel series (10). Unlike the direct Bessel-collocation approach described in the last section, this task is accomplished for the corresponding approximated quasilinear model problem (12). For this purpose, let us assume that an approximated solution ψ (r) M (τ) is known for the model problems (12) in the iteration r = 0, 1, . . .. In the next iteration, we consider Here, d k , k = 0, 1, . . . , M are the unknown coefficients that have to be sought. Let us define Utilizing (10), we can rewrite the finite series (24) in a matrix representation form where the vector of basis function Π Π Π M (τ) and the matrix C C C are defined in (9). Calling Algorithm 1, we calculate the second-order derivative of the vector of basis functions We then proceed by placing the collocation points (14) into (25) and (26) to write We now are aiming to find an approximate solution represented as (24) for the quasilinear model problem (12) through the Bessel-collocation approach. By collocating (12) at the collocation points (14), we obtain for r = 0, 1, . . .. To proceed, we define p r (t) = −α 2 cosh(αψ r (t)), f r (t) = p r (t) ψ r (t) + α sinh(αψ r (t)). Now, by using the matrix notations defined in (27) and (28), we may write the preceding equations as compactly asΨ Ψ Ψ r+1 + P P P r Ψ Ψ Ψ r+1 = F F F r .
In (29), the matrix P P P r as well as the vector F F F r can be written in the following forms .

Lemma 2.
Suppose that the approximate solution of (12) can be written as (25). Then, we obtain the following fundamental linear matrix equation Proof. By inserting the relations (27) and (28) into (29), we obtain the fundamental matrix equation as desired.
It can be obviously seen that, unlike Equation (22), the fundamental matrix equation (30) is a set of (M + 1) linear equations in terms of (M + 1) unknown coefficients d M to be determined. In each iteration, one requires to take into consideration the boundary conditions (13). The matrix representations of the boundary conditions (13) can be similarly done as for the non-linear counterpart. Besides, to begin the computation, we need to prescribe the initial guess ψ 0 (τ) as an approximation for the solution.
In a similar manner as the direct approach, the first and last rows of the augmented matrix [Σ Σ Σ r ; F F F r ] are replaced by the row matrices [ Σ Σ Σ r,0 ; 0] and [ Σ Σ Σ r,1 ; 1]. The resulting algebraic system of linear equations becomes which can be solved by any classical linear solver. Consequently, the unknown Bessel coefficients in (24) will be determined once we solve this system of equations.

Error Analysis
In this section, we make an investigation on error analysis for the Bessel series solution. In this respect, we present a theorem giving an upper boundary of errors.
Proof. Utilizing the triangle inequality after adding and subtracting the Maclaurin expansion ψ Mac M (τ) by M-th degree, we write the following expression According to ( Due to the fact that Π Π Π M (τ) ∞ ≤ 1 for τ ∈ [0, 1], we can arrange Equation (34) as follows On the other hand, since the remainder term of the Maclaurin series ψ Mac M (τ) by M-th degree is Combining Equations (33), (35), and (36), the proof of the theorem is completed.

Remark 2.
The upper bound obtained in Theorem 1 can be slightly modified. In fact, according to the remainder of the Taylor's theorem, there exist c τ ∈ (0, 1) such that Thus, we obtain The Accuracy of Methods The exact analytical solution of the non-linear model problem (1) is not known explicitly and only available in a closed-form. In this case, it is important to check the accuracy of the presented collocation methods. Let us assume that E M (τ) and E (r+1) M (τ) be the residual error functions, which are obtained by substituting the truncated Bessel series solutions (7) and (24)

Graphical and Computational Results
Here, we are going to illustrate the effectiveness of the presented direct Bessel collocation as well as QLM-Bessel collocation method for the non-linear model problem (1) through numerical experiments. All experimental computations have been carried out via the MATLAB R2017a software. For validation, comparisons with available existing numerical models are also performed. The corresponding approximate solution obtained via the QLM-Bessel approach (30) using the iteration parameter r = 5 is ψ (6) The initial guess in the latter approach is taken as ψ 0 (τ) = 0. Clearly, the outcomes of both approaches are the same while the direct collocation scheme takes more time, especially when M is increased and may be not be convergent at all. The above approximations along with the exact solutions are visualized in Figure 1. Note, the exact solution at some points τ ∈ [0, 1] is taken from [11]. These values are obtained via solving the closed-form solution (3). In Figure 1, on the right panel, we depict the computed residual error functions (38a) and (38b). Besides M = 5, the results corresponding to M = 8 are further visualized for both direct and QLM-Bessel collocation approaches.  Our results with M = 8 using both direct and linearization approaches are similar up to ten digits as follows ψ (6) 8 (τ), ψ 8 (τ) = 8.2519 × 10 −6 τ 8 + 1.8235 × 10 −5 τ 7 + 1.9425 × 10 −5 τ 6 + 9.4595 × 10 −4 τ 5 + 5.3536 × 10 −6 τ 4 + 0.03995880813 τ 3 + 1.9691 × 10 −7 τ 2 + 0.9590437817 τ − 6.276852 × 10 −17 .
The residual errors achieved via the QLM-Bessel method for Test case 1 using M = 5, 10, 15, and M = 20 are graphically visualized in Figure 2. It can be seen that the approximate solutions are convergent exponentially when M is increased. Finally, for this value of α, we justify our results by comparing to other existing available numerical models. In this respect, we employ M = 9 and report the numerical results in Table 1. We compare our computed solutions with the outcomes of the shifted Bessel Tau (SBT) method [20], the B-spline approach (BSA) [13], the double exponential sinc-Galerkin (DESG) method [15], and the modified non-linear shooting method (MNLSM) [9]. In this case, we first utilize M = 7 in the Bessel-collocation method to obtain In fact, using M = 10, the direct algorithm gives no reliable results. In the QLM-Bessel method, the corresponding approximate solution for 0 ≤ τ ≤ 1 is ψ (6) 10 (τ) = 5.9125 × 10 −4 τ 10 − 0.001726998 τ 9 + 0.0032385917 τ 8 − 0.001625926 τ 7 + 0.0020372912 τ 6 + 0.011214379 τ 5 + 2.4030 × 10 −4 τ 4 + 0.140824035 τ 3 + 4.6403 × 10 −6 τ 2 + 0.8452024383 τ − 5.267094243 × 10 −109 . Figure 3 shows the residual error functions achieved via the QLM-Bessel method for Test case 2. In these experiments, we utilize M = 5, 10, 15, and M = 20. The validation includes a comparison with existing available numerical and experimental results, reported for the methods used in Table 1, and is listed in Table 2.   Let us examine the impact of utilizing various values of parameter α on the obtained results. To this end, we fix M = 10 and let α vary as 0.125, 0.25, 0.5, 1, 2, 3, 4. The residual errors calculated via (38b) for τ ∈ [0, 1] are presented in Figure 4. It can be concluded that the error function E (6) 10 (τ) is a decreasing function with respect to α < 1, while it is inherently an increasing function with regard to α ≥ 1. In all cases, the maximum value of the errors occurred at τ = 1. Indeed, as investigated by Troesch [4], there exists a pole for the model problem (1), which is located at The values of τ s approach to τ = 1 if one increases α > 1. The solution profiles at diverse α = 0.125, 0.25, 0.5, 1, 2, 3, 4, and α = 5 are presented in Figure 5. Finally, we validate our results by comparing our results with an optimal iterative scheme, which is based on the homotopy analysis method and the Green's function technique (HAMG) [21].
In this respect, the residual errors E  Table 3. Looking at Table 3, one can infer that by increasing α the magnitude of errors achieved by our approach are getting smaller compared to the HAMG.
10 (τ) may be not so accurate as for the smaller values of α. To this end, we need to take M sufficiently large to obtain the desired approximations. The impacts of diverse M on the numerical solutions as well as the residual errors are shown in Figure 6. We use M = 10, 15, and M = 20 in the QLM-Bessel collocation method. Obviously, all approximated solutions are very close together while the corresponding residual errors are decreasing with respect to M.
Next, we compare the results obtained via the QLM-Bessel approach for this test case. In this respect, we employ M = 15, 25, and M = 35. In Table 4, comparisons have been made with the results of the Fortran code called TWPBVP [13] as well as with the B-spline approach (BSA) [13], the double exponential sinc-Galerkin (DESG) method [15], and the shifted Jacobi-Gauss collocation method (SJGCM) [18]. The behavior of numerical solutions obtained by the QLM-Bessel taking various (large) values of parameter α = 1, 3,5,8,10,15,20,30,40,50, and α = 100 are depicted in Figure 7. Due to difficulties arising in the approximate solutions for the larger values of α, we take M = 50 to keep the accuracy at some acceptable level.

Conclusions
An efficient collocation method based on the Bessel matrix representation combined with the quasi-linearization technique has been presented to find the solution of the two-point non-linear Troesch's problem arising in the modeling of a plasma confinement problem and gas porous electrodes. In addition, a direct Bessel collocation method has been developed for the underlying model problem. The error analysis of Bessel functions was established theoretically and the related convergence verified experimentally through numerical simulations on diverse values of the sensitive parameter α. Moreover, the comparisons of numerical solutions on different well-established schemes were performed to validate our presented approximation algorithms. It is found that the simulation results agree quite well in general with outcomes of existing schemes. The present QLM-Bessel method is broadly applicable and can be applied to solve diverse nonlinear model problems in engineering and science.