Lie-Group Shooting/Boundary Shape Function Methods for Solving Nonlinear Boundary Value Problems

: In the numerical integration of the second-order nonlinear boundary value problem (BVP), the right boundary condition plays the role as a target equation, which is solved either by the half-interval method (HIM) or a new derivative-free Newton method (DFNM) to be presented in the paper. With the help of a boundary shape function, we can transform the BVP to an initial value problem (IVP) for a new variable. The terminal value of the new variable is expressed as a function of the missing initial value of the original variable, which is determined through a few integrations of the IVP to match the target equation. In the new boundary shape function method (NBSFM), we solve the target equation to obtain a highly accurate missing initial value, and then compute a precise solution. The DFNM can ﬁnd more accurate left boundary values, whose performance is superior than HIM. Apparently, DFNM converges faster than HIM. Then, we modify the Lie-group shooting method and combine it to the BSFM for solving the nonlinear BVP with Robin boundary conditions. Numerical examples are examined, which assure that the proposed methods together with DFNM can successfully solve the nonlinear BVPs with high accuracy.


Introduction
There are numerous numerical methods to tackle the boundary value problems (BVPs) in [1][2][3][4][5][6][7][8][9][10]. The present paper will derive a very powerful numerical solver with a derivativefree Newton method to solve the target equation in terms of the Robin type boundary conditions. Recently, Liu and Li [11] designed the bases of solution by transforming the nonhomogeneous Robin boundary conditions to the homogeneous ones by using the homogenization function technique, and then seeking the polynomial bases to satisfy the homogeneous Robin boundary conditions. They introduced the energetic Robin boundary functions as the bases, whose processes are quite complicated.
For the first-order system of ordinary differential equations (ODEs), the group-preserving scheme (GPS) has been developed by Liu [12] for the solutions of initial value problems (IVPs). The main difference between the GPS and the traditional numerical integration methods is that those schemes are all formulated directly in the usual Euclidean space, while the GPS is formulated in the Minkowski space. A major advantage of the formulation of ODEs in the Minkowski space is that we can develop group preserving schemes to retain the inherent Lie-symmetry of the augmented ODEs in the Minkowski space. Recently, Xu and Wu [13] combined the GPS with a midpoint technique to improve the accuracy and convergence property. Liu [14] extended and modified the GPS for solving secondorder BVPs, based on the Lie-symmetry of the proper orthochronous Lorentz group, who introduced one-step GPS by utilizing the closure property of the Lie-group and labelled it as the Lie-group shooting method (LGSM). Moreover, Liu [15] developed the LGSM for solving nonlinear singularly perturbed boundary value problems. Recently, Hajiketabi and Abbasbandy [16] developed a simple, efficient and accurate Lie-group shooting method for solving nonlinear boundary value problems.
The paper is outlined as follows. In Section 2, the boundary shape function for automatically satisfying the Robin boundary conditions is introduced, where two cases are considered. In Section 3, we transform the BVP to a new IVP for the new variable y(x), with its terminal value y(1) being expressed as a function of the missing initial value of the original variable u(x). This new version is a modification of the original boundary shape function method (BSFM). In Section 4, we develop an iterative algorithm relying on the new boundary shape function method (NBSFM), and the half-interval method (HIM) is used to solve the target equation. In Section 5, the Lie-group shooting method is combined with the BSFM for developing a new iterative algorithm of LGSBSFM. Here, a main contribution is developing a derivative-free Newton method to solve the target equation. In Section 6, we test numerical examples by using these algorithms. We introduce the general Robin boundary value problems in Section 7, the corresponding NBSFM is derived, and two numerical examples are given. Finally, Section 8 makes the conclusions.

Boundary Shape Function
Consider the boundary value problem (BVP): A more general BVP with two-side Robin boundary conditions will be considered in Section 7. We first discuss the right boundary condition with u(1) = c 2 . For u (1) = c 2 , the processes are similar. For Equation (2), we consider two cases: (a) For cases (a) and (b), we have, respectively, A translation function G(x) in terms of s 1 (x) and s 2 (x) reads as where y(x) is a new variable. Theorem 1. The boundary shape function u(x) is given by which satisfies a 1 u(0) Proof. Inserting x = 0 into Equation (9) and using Equation (8), Taking the differential of Equation (9), inserting x = 0 and using Equation (8), we have A combination of Equations (11) and (12) leads to by Equations (4) and (5), which becomes Inserting x = 1 into Equation (9) and using Equation (8), which, due to Equations (4) and (5), becomes The proof of Equation (10) is ended.

The Initial Value Problem
For Equation (1), if very accurate initial conditions u(0) and u (0) are determined, one can obtain a very precise solution of u(x) by using a higher order numerical integrator to integrate the ODE. Therefore, it is an important topic to precisely determine B 0 := u(0) and A 0 := u (0) by an efficient numerical method. Based on Theorem 1, we can develop an iterative algorithm to solve Equation (1) subjected to Equations (2) and (3). Inserting Equation (9) for u(x) into Equation (1), we can achieve where G(x) involves y(1) as unknown constant in Equation (8), but y(0) and y (0) are given initial values.
In order to calculate the unknown value y(1) in Equation (16) and achieve a deeper analysis, we carry out the following works. If a 1 = 0, we assume that is a constant to be determined and u(0) is given by is a constant to be determined and u (0) is given by Under the condition a 1 = 0 for case (a), it follows from Equations (6), (9) and (17) that For case (b), it follows from Equations (7), (9) and (17), that On the other hand, under the condition b 1 = 0 and for case (a), it follows from Equations (6), (9) and (18) that For case (b), it follows from Equations (7), (9) and (18) that The above equations reveal that y(1) can be computed from A 0 or B 0 . Inserting these y(1) into Equation (16), we can obtain a definite ODE in terms of A 0 or B 0 : which, upon giving y(0) and y (0), are initial value problems (IVPs) endowed with one unknown parameter A 0 or B 0 to be determined. Notice that in the original BSFM developed in [17], y(1) in Equation (8) is determined by an iterative method, which repeatedly integrates Equation (16) to obtain the new value of y(1) until convergence. To improve the slower convergence and slightly lower accuracy of the original BSFM, we recast Equation (16) to Equation (23) or (24), which may be called a new boundary shape function method (NBSFM). Powerful methods will be developed below to determine A 0 or B 0 by solving the target equation. Deng et al. [18] have extended the BSFM to solve the nonlinear BVP under nonlinear boundary conditions.

The Iterative Algorithm of NBSFM
Because the initial values y(0) = y (0) = 0 are given, we can integrate Equation (23) or (24) from x = 0 to x = 1, with a given A 0 or B 0 to obtain y(1) and y (1); hence, inserting x = 1 into Equation (9) and using Equations (8) and (3) yields which can determine A 0 or B 0 by solving a scalar equation, where y(1) and y (1) are obtained by integrating Equation (23) or (24). The curve of u(1; A 0 ) vs. A 0 or u (1; B 0 ) vs. B 0 is a curve of a target function. We can adjust the value of A 0 or B 0 until u(1; where ε is a given tolerance of the error to mismatch the right-end boundary condition (3), which is an implicit equation to determine A 0 or B 0 . We can apply the half-interval method to solve Equation (25) or (26), where we repeatedly integrate Equation (23) or (24) with N steps from x = 0 to x = 1 to obtain y(1) and y (1), and then u(1; A 0 ) by Equation (25) The half-interval method (HIM) is described as follows. We can decide two initial guesses A 0 and A, such that u(1; A 0 )u(1; A) < 0. Then, we repeat the following process until convergence: In the half-interval method twice integrations to obtain the right-end values of u at each iteration are required.
The process to determine B 0 can be done similarly, by using Equations (24) and (26).

A Combination of LGSM and BSFM
The original Lie-group shooting method developed by Liu [14] was applicable to the nonlinear BVP with simple Dirichlet or Neumann boundary conditions. For solving Equations (1)-(3), we can combine the Lie-group shooting method and the BSFM as a new Lie-group shooting boundary shape function method (LGSBSFM), which is used to solve y(x) in Equation (16). Here, with case (a) and a 1 = 0 as a demonstrative case, we let where A and B are two unknown constants, and c is a given constant. Based on the Lie-group shooting method [14], without going into the details, we can derive A for the case A > 0: wherê and for the case A < 0: For a trial r ∈ [0, 1], we can iteratively solve A from Equations (33), (34), and (37). The advantage is that we can determine the slope of the new variable y(x) in terms of a factor r within a unit interval, which, however, has a drawback of the Lie-group shooting method to iteratively solve A for each r. As a modification, we directly determine A by integrating Equations (27)-(30) to obtain the right-end value y(1, A), and then compare it to the desired one y(1) = c. If |y(1, A) − c| is smaller than a given error tolerance , then the process to find A is finished. Similarly, the iteration process can be performed by HIM to find a suitable value of A by starting from two given values A 0 and A 2 with . On the other hand, according to Equations (19), (29), and (30), we can derive wherein c is eliminated and y (0) = A is used. The initial conditions for u(x) are available as follows: In addition, then, we can integrate Equation (1) using the derived initial conditions u(0) and u (0) to obtain u(x).

A Derivative-Free Newton Method
We can counter the number of iterations for the half-interval method (HIM), which needs 50 iterations for reducing an initial error 1 to a final error with 1/2 50 = 8.882 × 10 −16 . If the error tolerance is = 10 −15 to satisfy the right boundary condition, the HIM will spend 50 × 2N steps in the numerical integrations of the ODE to obtain the right-end values. To reduce the computational burden, a derivative-free Newton method (DFNM) for solving a scalar equation f (x) = 0 is motivated by the Newton iterative method: Unfortunately, the derivative term f (x n ) hinders the use of the Newton method in LGSB-SFM, since the scalar target equation y(1, A) − c = 0 is an implicit function of A.
Suppose that x * is a simple root with f (x * ) = 0 and f (x * ) = 0. In order to get rid of the derivative term f (x n ) in Equation (40), we consider Neglecting the higher-order terms and inserting it into Equation (40), we have On the other hand, we have Neglecting the higher-order terms and replacing we can derive a derivative-free Newton method (DFNM): where Next, we turn our attention to the determination of a and b in Equation (45), whose values will influence the convergence speed. Like the half-interval method, the first step is choosing two initial guesses x 0 and x 2 such that f (x 0 ) f (x 2 ) < 0 to render x * ∈ (x 0 , x 2 ). Then, we take x 1 = (x 0 + x 2 )/2. As the approximations of a and b in Equation (45), we can evaluate them by using the technique of finite difference: In summary, the procedures of the DFNM are given as follows: (i) Give the initial guesses of x 0 and x 2 to render f (x 0 ) f (x 2 ) < 0.
(ii) Compute a and b by Equation (46). (iii) For n = 0,1,. . . , doing The LGSBSFM together with the DFNM for solving u(x) in Equations (1)-(3) is summarized as follows: (i) Give y(0) = c, the initial guesses of A 0 and A 2 to render [y(1, and give , and ∆x = 1/N. A 1 ), and a and b by (iii) Let A = A 0 and for n = 0,1,. . . , doing In each iteration, the integration of Equations (27)- (30) to obtain the right-end value y(1, A n ) is required, which is saving more than twice the integrations at each iteration used in the HIM. For saving computational time, the DFNM is suggested to be used.
For solving a scalar equation f (x) = 0, the numerically computed order of convergence (COC) is approximated by [19,20] where r is a solution of f (x) = 0 and the sequence x n is generated from an iterative scheme.
In the computation of COC, we store the values of A n where n ≤ k 0 − 1, and take r = A k 0 , with k 0 the number of iterations for the convergence.

Examples
To solve Equations (1)-(3), we have developed two iterative algorithms of the NBSFM and LGSBSFM, which can be combined with the HIM or DFNM to solve the target equation. Some examples are given below to explore the performance of these algorithms.

Example 1
The following BVP is given in [21]: An exact solution is (50) Instead of Equation (49), we consider a left Robin boundary condition: Due to a 1 = 2, b 1 = 1, and a 1 − b 1 = 1, this is a case (a) problem in Sections 2 and 3. By using the conventional shooting method, we take where A 0 is an unknown constant to match the target equation u(1) = 1. We can employ HIM to solve the scalar equation u(1) − 1 = 0. Under N = 2000 and = 10 −15 , we find that the conventional shooting method does not converge within 1000 iterations, although the solution of u(x) is accurate with the maximum error (ME) being 2.44 × 10 −13 .
In the NBSFM with HIM, we take s 1 (x) = 1 − x, s 2 (x) = 2x − 1, and Under the parameters y(0) = y (0) = 0, N = 2000, and = 10 −15 , the NBSFM with HIM converges with 42 iterations with ME = 2.55 × 10 −13 . The result is more accurate than that in [14]. If we raise N = 4000, the ME is reduced to 2.31 × 10 −14 , and it is with 48 iterations for convergence. If we employ the NBSFM with DFNM, it is convergent with 14 iterations to obtain the same results. Since only 14 iterations are required, the CPU time is short with 0.6 s. Although with the HIM the CPU time is 0.9 s, its number of iterations is 42. COC = 1.6122 obtained by DFNM reveals that the present method is convergent very fast. Next, we apply the LGSBSFM in Section 5 together with HIM to solve this problem with c = 0, A 0 = −5.1, A 2 = −4.95, N = 5000, and = 10 −15 , which is convergent with 48 iterations, and the numerical solution coincides wth the exact one as shown in Figure 1a, whose errors are plotted in Figure 1b with ME = 1.13 × 10 −14 . If we employ the LGSBSFM with DFNM, it is convergent with 10 iterations to obtain the same results with ME = 1.13 × 10 −14 . The CPU time is very short with 0.46 s, and COC = 1.2546 is obtained, which indicates that the LGSBSFM with DFNM is effective. When we take c = 1, A 0 = −0.8, A 2 = 0, we obtain the second solution as shown in Figure 1a by dashed line. The LGSBSFM with HIM is convergent with 47 iterations and with the right boundary error being 1.44 × 10 −15 . If we employ the LGSBSFM with DFNM, it is convergent with 16 iterations and with the right boundary error being 2.22 × 10 −16 .
In Table 1, we tabulate the absolute errors and compare that obtained in [17]. The accuracy upon comparing to that computed in [17] using the boundary shape function method is raised about five orders. In Table 2, we tabulate the numerical value at different x and compare it to the exact one. The accuracy is up to 10 −15 to 10 −14 . Next, we consider other boundary conditions for Equation (48): Due to a 1 = 1, b 1 = 1, and a 1 − b 1 = 0, this is a case (b) problem in Sections 2 and 3. By using the NBSFM with DFNM, s 1 (x) = 1 − x 2 and s 2 (x) = 1 − x + x 2 and y(1) is given by

Example 2
Consider a reaction problem [22]: For this problem, we assume that B 0 in Equation (18) is a constant to be determined and u (0) is given by u (0) = Pe(B 0 − 1). By using the NBSFM with the DFNM, we take s 1 (x) = 1/Pe and s 2 (x) = 1/Pe + x and y (1) is given by Under the parameters Pe = 1, R = 3.5, n = 2.5, y(0) = y (0) = 0, N = 4000, and = 10 −15 , the NBSFM with HIM by two initial guesses B 0 = 0.6 and B 2 = 0.61 to render u (1, B 0 )u (1, B 2 ) < 0 converges with 45 iterations with the right boundary error being 7.8 × 10 −15 . In contrast, by using the NBSFM with DFNM, it converges with seven iterations with the same right boundary error 7.8 × 10 −15 . The numerical solution is plotted in Figure 3. Since only seven iterations are required, the CPU time is very short with 0.35 s, and COC = 1.0044 is obtained by DFNM.
Then, we apply the LGSBSFM together with HIM or DFNM to solve this problem with c = 0.95, A 0 = −0.76, A 2 = −0.66, N = 10000, and = 10 −15 , which is convergent with 39 iterations for the HIM, and with nine iterations for the DFNM. The numerical first solution coincides with the exact one in Equation (62) with c 0 = 0.1836751649400605 with ME = 3.22 × 10 −11 , which is plotted in Figure 4 with the dashed line. When we take c = 1, A 0 = −0.473, A. 2 = −0459, and LGSBSFM together with DFNM converges with seven iterations. The second numerical solution coincides with the exact one in Equation (62) with c 0 = 0.5330478767543923 with ME = 9.02 × 10 −13 , which is plotted in Figure 4 by the solid line. The accuracy is better than that obtained in [16] about five orders.
With five iterations as shown in Figure 5 for the convergence, we obtain the solution as shown in Figure 5b by a solid line. Since only five iterations are required, the CPU time is very short with 0.49 s, and COC=1.0204 reveals that the present method is convergent fastalthough, for HIM, the CPU time is 0.87 s. The numerical solution coincides with the exact one in Equation (65), and the numerical error is shown in Figure 5b by a dotted-dashed line with ME = 3.84 × 10 −14 . In Table 3, we tabulate the absolute errors and compare that obtained in [29][30][31].

Two-Side Robin Boundary Conditions and Examples
In this section, we extend the NBSFM and DFNM to the following nonlinear BVP under the two-side Robin boundary conditions:

A New Methodology
As that done in [17], two shape functions are determined by Then, the following variable transformation is employed: which guarantees that u(x) satisfies Equation (68) when y(x) is solved.
In [17], the authors determine y(1) and y (1) in Equation (71) by the iteration method until the convergence of y(1) and y (1). Here, we let α := a 2 y(1) where α is an unknown constant to be determined. As before, we let u(0) = A be an unknown constant to be determined, and, after inserting x = 0, it follows from Equation (73) that Like that in Equations (19)- (22), where y(1) is expressed in terms of A 0 or B 0 , for the present BVP with a right Robin boundary condition, α = a 2 y(1) + b 2 y (1) as a combination of y(1) and y (1) is expressed in terms of A := u(0) in Equation (74) or B := u (0) as follows: which is obtained by inserting x = 0 into the differential of Equation (73). Inserting Equation (74) into Equation (73) yields where In addition, then inserting u(x) into Equation (67), a new ODE for y(x) is generated, like that in Equation (16), which through G(x) involves A as a parameter. Then, we employ the DFNM in Section 5.2 to iteratively determine A by matching the Robin type target equation: The numerical procedure is similar if we employ Equation (75) to express α. We can also develop a similar LGSBSFM for solving the BVP with two-side Robin boundary conditions by where c and A are given constants and u(0) = A 0 is determined by matching the Dirichlet type target equation:
By using the LGSBSFM and DFNM with c = −1 and A = −1, it is convergent with eight iterations. The numerical solution coincides with the exact one in Equation (50) with ME = 6.12 × 10 −8 . The CPU time is very short with 0.35 s for eight iterations, and COC = 1.17223 is obtained by DFNM.

Example 6
We calculate a nonlinear singular perturbation problem: The exact solution is Using the original boundary shape function method through 221 iterations, Liu and Chang [17] find the numerical solution with ME = 1.993 × 10 −4 for ε = 0.02. Here, we consider a severely singular case with ε = 0.001. With s 1 (x) = (2 + ε − 2x)/(4 + ε) and s 2 (x) = (1 + x)/(4 + ε), y(0) = −1, y (0) = −1, N = 5000, and = 10 −12 , and (A 0 , A 2 ) = (−0.1, 0.1), the DFNM is convergent with six iterations as shown in Figure 6a. We obtain the solution as shown in Figure 6b by a solid line, which coincides with the exact one in Equation (85), and, in Figure 6b, the numerical error is shown by a dotted-dashed line with ME = 1.35 × 10 −4 . The CPU time is very short with 0.32 s for six iterations, and COC = 1.5375 is obtained by DFNM. In Table 4, we tabulate the numerical value at different x and compare it to the exact value. Besides at the point within the boundary layer, the accuracy is up to 10 −14 to 10 −13 .

Conclusions
In the numerical integration of nonlinear BVP, a key issue is how to determine the missing initial value A 0 or B 0 quickly and accurately. In the present paper, we have further modified the boundary shape function method (BSFM) for a new variable, of which the resulting initial value problem with its ODE involves the missing initial value of the original variable as a parameter. Owing to the implicit and nonlinear property of the target equation, we employed a new derivative free Newton method (DFNM) to solve the target equation, which can find very accurate left boundary values within a few iterations, and its performance is better than the half-interval method. In the new boundary shape function method (NBSFM), we solve a target equation to render a highly precise missing initial value and then the solution obtained is very accurate. Through the tests of numerical examples, we found that the CPU time is less than one second, and the COC is located between one and two. We combined the techniques of the Lie-group shooting method to the BSFM, whose advantage is that we can determine the missing slope of the new variable or missing initial value of the original variable in terms of a simple Dirichlet type target equation. Furthermore, we can easily find the multiple solutions of the considered problems. These numerical results clearly showed that the proposed methods can provide excellent approximation to the true solution of the nonlinear BVP with high accuracy.