Three new approaches for solving a class of strongly nonlinear two-point boundary value problems

Three new and applicable approaches based on quasi-linearization technique, wavelet-homotopy analysis method, spectral methods, and converting two-point boundary value problem to Fredholm–Urysohn integral equation are proposed for solving a special case of strongly nonlinear two-point boundary value problems, namely Troesch problem. A quasi-linearization technique is utilized to reduce the nonlinear boundary value problem to a sequence of linear equations in the first method. Second method is devoted to applying generalized Coiflet scaling functions based on the homotopy analysis method for approximating the numerical solution of Troesch equation. In the third method we use an interesting technique to convert the boundary value problem to Urysohn–Fredholm integral equation of the second kind; afterwards generalized Coiflet scaling functions and Simpson quadrature are employed for solving the obtained integral equation. Introduced methods are new and computationally attractive, and applications are demonstrated through illustrative examples. Comparing the results of the presented methods with the results of some other existing methods for solving this kind of equations implies the high accuracy and efficiency of the suggested schemes.


Introduction
Troesch problem arises from a system of nonlinear ordinary differential equations which occur in the theory of gas porous electrodes [1] and investigation of the confinement of a plasma column by radiation pressure [2]. This problem is a special type of nonlinear two-point boundary value problems, defined as follows: subject to the following conditions: where the Troesch parameter λ is a positive constant. The existence of the solution of (1)- (2) was proved in [3] for λ < 5. The closed form solution to this problem in terms of the Jacobian elliptic function was given [4] as where ν (0) = 2 √ 1μ, and μ satisfies the transcendental equation where the Jacobian elliptic function sc(λ | μ) = tan α, and α and λ are related through the following integral: It is clear that ν(x) has a singularity located at a pole of sc (λx|μ) or approximately at [5,6] x s = 1 λ ln 8 ν (0) , which implies that, if ν (0) > 8e -λ , then the singularity lies within the integration range. Because of difficulty in direct solving of the nonlinear boundary value problems, many researchers have paid considerable attention to numerical solving of these problems. For example, the numerical solution of one-dimensional Bratu's problem is solved by an operational matrix of the derivative of Chebyshev wavelets in [7]. Authors of [8] used a manipulation in the cubic splines to develop a collocation method and the generalized Newton method for solving the nonlinear Troesch problem. Sinc-Galerkin approach is applied to numerical solving of second-order nonlinear Dirichlet-type boundary value problems in [9]. In [10] and [11] some approaches based on B-splines are presented for solving the Troesch problem. The homotopy perturbation technique is presented in [4], and the reproducing kernel Hilbert space method is introduced in [12]. Sinc-Galerkin method based upon double exponential transformation is proposed to solve Troesch's problem in [13]. Wavelet theory is a relatively new and emerging area in mathematical research. In recent decades, wavelets have been extensively used in applied mathematics and a wide range of engineering disciplines, because they permit accurate representation of a variety of functions and operators. Wavelets were applied in numerically solving the integro-differential and differential equations [14,15]. Also the hybrid methods, combination of some analytical and numerical methods, were successfully employed to achieve the solutions of initial and boundary value problems; for instance, see [16,17,26,27]. In this paper we apply three different methods based on a quasi-linearization technique, the homotopy analysis method based Coiflet orthogonal scaling functions, and an efficient approach for transforming the aforementioned equation to an integral equation.
The outline of the paper is as follows: Sect. 2 is devoted to the brief definition of Coiflet scaling functions and function approximation by them. Suggested computational methods for solving Troesch equation are discussed in Sect. 3. Indeed, in this section three different methods are introduced based on quasi-linearization, operational matrices of derivative and collocation and Galerkin methods and integral equation transform technique. Convergence analysis is discussed in Sect. 4. To show the accuracy and performance of the presented methods, some test problems are given in Sect. 5. The conclusion of the paper is described in Sect. 6.

Generalized Coiflet orthogonal wavelets
The basic concepts and preliminaries of the wavelet theory and multi-resolution analysis are given in literature [18,19,[28][29][30]. Scaling function φ and wavelets of generalized Coiflet-type orthogonal wavelet possess the following properties [20,21]: p k are the low-pass filter coefficients [22], N is the number of vanishing moment, and M n is the nth-order moment of the scaling function. In this paper, we apply the generalized Coiflet type with N = 6 and M 1 = 7. The general difficulty for working with this family of wavelets is having no analytical expression for them. So the values of these functions are calculated [23]. Figure 1 shows the plot of Coiflet scaling function, and Table 1 represents some values of Coiflet scaling function in k = 1, . . . , 16. By the definition of Coiflet scaling function, we can construct a multi-resolution analysis for L 2 (R) with nested subspaces V j as follows: where V j = span{ϕ j,k := ϕ(2 j xk); j, k ∈ Z}. In this paper we introduce a classification for generalized Coiflet scaling functions on [0, 1], which have not been addressed previously. For this purpose we first define the following matrices: and

Boundary scaling adaptation • Left boundary Coiflet scaling functions
We define the boundary near functions at the left boundary for k = 0, . . . , 3 by • Right boundary Coiflet scaling functions At the right end of the interval, for k = 2 j -3, . . . , 2 j , we have Interior scalings For interior Coiflet scaling functions, we have Figure 2 shows the family of Coiflet scaling functions for j = 5.

Function approximation
A function f (x) ∈ L 2 [0, 1] can be approximated in the subspace V j as [24] f where F and j are row and column vectors defined by It should be noted that by using the orthogonality property of Coiflet, the approximation coefficients could be obtained as follows: where ξ j,k = ϕ j,k (x), ϕ j,k (x) . For example, the function sin(x) is approximated by Coiflet in V 4 and V 5 , and the plots of absolute error of approximation are shown in Fig. 3. Also the n-th derivative of function f (x) can be approximated by the n-th derivatives of Coiflet scaling functions as follows: where C is a positive constant that depends only on the function f (x) and low-pass filter coefficients ρ k .
Proof By using the Taylor expansion of f at x, we get Putting y = k 2 j in equation (11) and substituting it in equation (8), we have we know that ϕ has N vanishing moments, that is, We know Supportϕ = [0, 3N -1], thus Support k = [k, k + 3N -1], that is, k has compact support. So there exists a positive constant C such that , the accuracy of derivative approximation can be estimated by Proof It is sufficient to put f (n) (x) = g(x) and repeat the above proof.

Methods
In this section we introduce three computational methods for numerically solving a strongly nonlinear two-point boundary value problem, namely Troesch problem. In the first method, some quasi-linearization technique is applied to reduce nonlinear boundary value problem (1)-(2) to a sequence of linear equations. Generalized Coiflet scaling functions via HAM and Galerkin method are used to reduce the nonlinear equation to some algebraic system in the second approach. Utilizing an interesting technique, the third approach is divided to transform the given nonlinear boundary value problem to some Urysohn-Fredholm integral equation of the second kind, and the obtained integral equation is solved applying Coiflet scaling functions and the Simpson quadrature rule.

Approach I: quasi-linearization technique
Consider boundary value problem (1)- (2), in this method we use a quasi-linearization technique based on Taylor expansion for reducing the main nonlinear problem to a sequence of linear algebraic equations. For this purpose, we choose a reasonable initial approximation for the function u(x), call it u 0 , and we develop ω(x, ν) = sinh(λν(x)) by the Taylor expansion around ν 0 , so we get or in a general form we can rewrite the current equation for n = 0, 1, 2, . . . (n = iteration index) By truncating the infinite series in equation (13) for arbitrary M, we get where (14) in equation (1), we have By truncating the obtained equation for arbitrary M, we get where ξ 1 , ξ 2 , . . . and ξ are determined similarly to η. According to current equations, nonlinear boundary value problem (1) converts to equations (15), a sequence of semilinear differential equations. It is interesting to point out that for suitable choices of ν 0 equations (15) may be reduced to linear differential equations. The obtained system of semilinear differential equations with boundary conditions (2) can be solved directly or numerically.

Approach II: wavelet homotopy analysis method
In this approach we use HAM based on generalized Coiflet scaling functions for solving equation (1). In order to simplify the calculation, we put u(x) = e λυ(x) , then strongly nonlinear problem (1)-(2) can be written as The basic idea of HAM is constructing a continuous variation between the known initial guess u 0 and the unknown solution u. Let q ∈ [0, 1] denote the embedding parameter of homotopy, c 0 the convergence-control parameter, L the auxiliary linear and N the nonlinear operators of HAM, respectively. Then a continuous variation (with respect to the embedding parameter q from the initial guess u 0 (x) to the solution u(x)) can be built by means of the so-called zeroth-order deformation equation as with α(0; q) = 1, α(1; q) = e λ , and Now we write Maclaurin series of α(x; q) with respect to q, so ∂q j | q=0 , and D j is the jth homotopy derivative operator. By this assumption that the convergence control parameter c 0 and the auxiliary linear operator L are chosen that guarantee the convergence of series (19) at q = 1, we get If we truncate the infinite series (20) in suitable m, we obtain the mth-order approximation as Taking the operator D m on both sides of the zeroth-order deformation equation (18) and its boundary conditions, we have and and χ 0 = χ 1 = 0 and χ j = 1, j > 1. Therefore nonlinear problem (1)-(2) reduces to a system of recursive equations. For solving this system, we utilize the Galerkin method with generalized Coiflet scaling functions as test and weight functions. For this purpose, first U m and R m are written in generalized Coiflet scaling functions terms in scale j as equation (8) U m (x) =  Substituting these approximations in equation (22), we get Now we utilize the Galerkin method to solve the obtained system. By multiplying the both sides of equation (23) in ϕ j,i , i = 0, 1, . . . , 2 j , and integrating in the interval [0, 1], we get The matrix form of system (24) is as where From the boundary conditions it is clear that U k (0) = U k (1) = 0. It should be noticed that because of orthogonality and having compact support of generalized Coiflet scaling functions, the matrices and are so sparse. Figure 4 shows the gray scale plots of matrix for j = 4, 5. By solving the system of equations (25) for m = 1, 2, 3, . . . step by step and substituting obtained U m in equation (21), we attain the approximate solution for u(x). Finally, in order to gain the solution of problem (1)-(2), we put ν(x) = 1 λ Ln(u(x)).

Approach III: converting to integral equation
In this method, using an interesting technique, nonlinear boundary value problem (1)- (2) is converted to a second kind nonlinear Fredholm integral equation of Urysohn type. For this purpose, we integrate equation (1) two times, so we get Reducing the double integral to a single integral, we have For evaluation ν (0) in equation (26), we put x = 1 and get By partitioning the integration interval to [0, x] and [x, 1] and substituting the result in equation (26), we have where

Equation (27) is a nonlinear second-order Fredholm integral equation of the Urysohn type.
For solving this equation, we apply generalized Coiflet scaling functions. Using equation (8), we approximate the unknown function as ν(x) = ν j (x) and substitute it in equation (27). Now we utilize the collocation method with mesh points τ i = i 2 j , i = 0, 1, . . . , 2 j , so we have The integral term in equation (28) can be found by applying some quadrature rule. We use the Simpson integration rule for evaluating the mentioned term. For this purpose, we consider the following partition for integration interval [0, 1]: τ , γ 2k , ν j (γ 2k ) .
Therefore nonlinear boundary value problem (1)-(2) converts to some linear or nonlinear algebraic system. In this paper the obtained system is solved using Wolframe Mathematica Programming 10.

Convergence analysis
In this section, by using Theorem 2.1, we present a theorem for error bound of approximate solution obtained by approach III.
where α is constant and depends on the function ν, and j is the resolution scale of generalized Coiflet scaling functions.
Proof Substituting ν * in equation (27) and subtracting two relations, we get By considering the definition of , we prove the theorem for two following cases.
By the mean value theorem we can write where A = sup{|λν (t) cosh(λν (t))|; 0 ≤ t ≤ 1}, so we have Considering Theorem 2.1, we can write It is clear that • Suppose t ∈ [x, 1], then Similar to the proof of the pervious case, we can write Again considering Theorem 2.1, we can write It is clear that Therefore, by putting α = AC 3 , the proof is complete.

Case 2. Approach I for ν 0 = 1
Numerical solutions of equation (1) for different values of λ, initial solution ν 0 = 1, and iteration index n is evaluated by approach I, and results are given in Figs. [11][12][13]. As the pervious case results for λ ≤ 1 are so close together that it is impossible detect them by plot, thus in Fig. 11 we show the error plot of ε 0.1,0.25 (x), ε 0.25,0.5 (x), ε 0.5,0.75 (x) and ε 0.75,1 (x) for iteration index n = 1. In Fig. 12 the numerical solutions of equation (1) obtained by approach I are shown for the Troesch parameter λ = 1, 2, 3, 4, 5. Figure 13 is useful for geometric understanding of the numerical solutions for λ = 6, . . . , 20. As we can see in this plot for 5 ≤ λ, almost approximated solutions have a linear form and are close together. Also it is considerable that the solutions function ν(x) is increasing with respect to the Troesch parameter λ.

Numerical results of approach II
The performance of this approach is analyzed for the auxiliary linear operator L(α(x; q)) = ∂ 2 α(x;q) ∂ 2 x , the resolution level j = 3, 4, 5, the iteration index m = 4, 5, and the convergencecontrol parameter c 0 = -0.5. The attained results of approach II for λ = 0.5 are given in Table 2 and are compared with the solutions of decomposition method (DM) of [4] and the spline collocation method (SCM) of [8]. Table 3 contains the numerical solutions for scales m = 4, 5 and the Troesch parameter λ = 1 and also the results of variational iteration method (VIM) of [4] and the coupled Laplace transform and modified decomposition method (LT-MDM) of [11]. In Table 4 numerical solutions of problem (1)-(2) for λ = 5 are given in some arbitrary points and compared with the results of Fortran code (FD) [25] and B-spline collocation method (B-SCM) [25]. Also the effect of resolution scale j and iteration index m for λ = 1 are presented in Tables 5-6 by means averaged square error E(m, j), defined as follows: where ν is the exact and ν m,j is the approximated solution obtained by approach II at scale j for iteration index m. As we can see by choosing higher scale j, the relevant E(m, j) is de- creasing exponentially, while the required CPU time is increasing linearly. Also the results of Table 6 show that, for different values of m, the CPU time is nearly fixed. It is interesting that the required CPU time for j = 5 and m = 10 is only 1.2 seconds. Of course this result was predictable because of sparsity of operational connection coefficients matrices.

Numerical results of approach III
Using approach III, we solve nonlinear boundary value problem (1)-(2) in scale j = 5 for different values of the Troesch parameter λ and compare our findings with the results of some other existing methods. The absolute errors of numerical results evaluated by approach III are shown in Tables 7-8 for different values of the Troesch parameter. In Table 7 absolute errors relevant to approach III for m = 4, 5, the modified homotopy perturbation technique (MHPT) [4], and the Laplace transform and a modified decomposi- tion technique (LT-MDM) [11] are reported for Troesch parameters λ = 0.5 and 1. Table 8 contains absolute errors of numerical solutions evaluated by approach III, sinc-Galerkin method based upon double exponential transformation (DESG) [13], and B-SCM [25] for λ = 5. For having a comparison between the results of presented methods, L 2 -errors of numerical solutions are given in Table 9, where L 2 -error is defined as follows. Suppose that ν and ν * are the exact and numerical solutions of problem (1), respectively, we put

Conclusion
In this paper three new and applicable methods are implemented on a class of strongly nonlinear differential equations. The proposed methods are based on some simple quasilinearization technique, the generalized Coiflet scaling functions based homotopy analysis method, and Galerkin and also some interesting approach for converting a boundary value problem to an integral equation. A common property of the presented method was reduc-    ing the nonlinear equation to a system of algebraic equations. Comparison between our results and the findings of some other existing methods for solving this kind of problems shows high accuracy and efficiency of the methods. It is considerable that the results of 2.0 × 10 -9 0.5 3.5 × 10 -8 6.5 × 10 -11 --0. 6 6.6 × 10 -8 7.4 × 10 -11 --0. 7 1.1 × 10 -7 2.8 × 10 -11 --0. 8 5.9 × 10 -7 1.7 × 10 -11 1.4 × 10 -2 1.6 × 10 -7 approach I for iteration index n = 1 are as accurate as some methods for large iteration index. Also because of some significant properties of generalized Coiflet scaling functions, such as orthogonality, having compact support, and vanishing moments, the operational matrices in approaches II and III are so sparse, and consequently relevant required computational time and memory is so low.

Future remarks
The presented methods in this study are attractive and can be extended for similar highorder nonlinear differential equations. Also, they could applied on the fractional-order differential and integro-differential equations with little additional work. Also the mentioned approaches could be implemented for solving nonlinear dynamic systems, including differential equations. Further research along these lines is under progress and will be reported in due time.