Boundary shape function method for nonlinear BVP, automatically satisfying prescribed multipoint boundary conditions

It is difficult to exactly and automatically satisfy nonseparable multipoint boundary conditions by numerical methods. With this in mind, we develop a novel algorithm to find solution for a second-order nonlinear boundary value problem (BVP), which automatically satisfies the multipoint boundary conditions prescribed. A novel concept of boundary shape function (BSF) is introduced, whose existence is proven, and it can satisfy the multipoint boundary conditions a priori. In the BSF, there exists a free function, from which we can develop an iterative algorithm by letting the BSF be the solution of the BVP and the free function be another variable. Hence, the multipoint nonlinear BVP is properly transformed to an initial value problem for the new variable, whose initial conditions are given arbitrarily. The BSF method (BSFM) can find very accurate solution through a few iterations.


Introduction
Many engineering problems can be modeled by ordinary differential equations (ODEs). When they are subjected to prescribed boundary conditions, we encounter the boundary value problems (BVPs), which manifest themselves in many applications, for instance, engineering, control theory, and optimization. For the details of the conditions for the existence and uniqueness of solutions of second-order BVPs, we refer to [1][2][3].
The multipoint BVPs arise when the states of an ODE system are measured at many points, which are important in many areas of engineering applications. The multipoint BVPs have attracted a lot of researchers. For three-point BVP of the second-order ODE, Ahmad et al. [4] adopted the quasilinearization method to obtain a monotone sequence, which converged quadratically to a solution. After that, Henderson [5] developed a double fixed-point theorem, applied to yield the existence of at least two nonnegative solutions for the second-order three-point BVP. Then, Sun and Liu [6] investigated the existence of a nontrivial solution for the three-point BVP. Several sufficient conditions for the existence of nontrivial solution were obtained by using Leray-Schauder nonlinear alternative.
numerical approximation of solutions are relatively rare [24][25][26][27]. It is desired that the numerical solution of the nonlinear multipoint BVP can exactly satisfy the prescribed boundary conditions, but in the case of nonseparable multipoint boundary conditions, it might be a difficult task. A common disadvantage of the above-mentioned literature is that it did not investigate the robustness of proposed schemes. Moreover, they are not guaranteed to satisfy the multipoint boundary conditions, automatically.
In the paper, a novel method based on the new concept of shape function and boundary shape function is derived for solving the second-order nonlinear BVP under nonseparable multipoint boundary conditions. We arrange the paper as follows. In Sect. 2, we introduce two shape functions and a boundary shape function, which is designed for automatically satisfying the boundary conditions prescribed at several points, where some main results are shown. In Sect. 3, an iterative algorithm, namely the boundary shape function method (BSFM), is developed, and some examples are given in Sect. 4. The conclusions are described in the last section.

Boundary shape function
For the solution of the following nonlinear second-order boundary value problem (BVP), endowed with prescribed nonseparable multipoint boundary conditions: we propose a new iterative method. In above, u(x 1 ), u (x 1 ), . . . , u(x m ), u (x m ) are respectively the values of u(x) and u (x) at m different points x 1 < · · · < x m . Here, [x 1 , x m ] is an interval of our problem. Since the boundary conditions are specified at m distinct points, this problem is called an m-point BVP; L 1 and L 2 are linear operators, acting on [u(x 1 ), u (x 1 ), . . . , u(x m ), u (x m )] by where c ij , i = 1, 2, j = 1, . . . , 2m are given constant coefficients, not all being zeros.
Proof Beginning with we prove the existence of s 1 (x). From Eqs. (4)- (6) and (8), a and b are determined by Obviously, when the rank of the coefficient matrix is two, there exists a unique solution for a and b.
In the situation when the rank of the above coefficient matrix is one, we can, instead of Eq. (8), take which results to a consistent system: such that we have a unique solution for a and b. Therefore, there exists a solution (a, b), and hence the solution s 1 (x) in Eq. (8) or Eq. (10) exists. Similarly, we can do it for s 2 (x).

Theorem 2.2 For a given free function f
, if s 1 (x) and s 2 (x) satisfy Eqs. (6) and (7), then can be defined and satisfies Proof In Theorem 2.1, the existence of s 1 (x) and s 2 (x) renders the existence of B(x), wherein f (x) ∈ C 1 [x 1 , x m ] is a given free function.
Applying the linear operator L 1 to Eq. (12) on both sides and using the linearity property, we have which, with the help from Eqs. (6) and (7), becomes Similarly, applying the linear operator L 2 on both sides of Eq. (12) and using the linearity property, we have which, with the help from Eqs. (6) and (7), becomes Thus, the proof of Eqs. (13) and (14) is completed. Theorem 2.2 is crucial, from which the treatment of very complex nonseparable multipoint boundary conditions for the nonlinear BVP becomes easy because the boundary shape function B(x) is guaranteed to satisfy the multipoint boundary conditions exactly and automatically.

The numerical algorithm
Utilizing the boundary shape function (BSF), the developed iterative algorithm to solve Eqs. (1)-(3) is given below. According to Theorem 2.2, B(x) given in Eq. (12) satisfies the multipoint boundary conditions in Eqs. (2) and (3). Thus, we may transform u(x) to y(x) by For any function y(x) ∈ C 2 [x 1 , x m ], u(x) automatically satisfies the multipoint boundary conditions (2) and (3). Inserting Eq. (15) for u(x) into Eq. (1), we achieve which can be viewed as an initial value problem (IVP), whose initial values are given arbitrarily, say, y(x 1 ) = y (x 1 ) = 0. The function H is given by where There are a number of unknown parameters y( Letting from Eq. (16) it follows that which are subjected to the given initial conditions y 1 (x 1 ) and y 2 (x 1 ). However, z := [y 1 (x 2 ), y 2 (x 2 ), . . . , y 1 (x m ), y 2 (x m )] T is an unknown vector. If z is available, we can apply the fourth-order Runge-Kutta method (RK4), as shown in the Appendix, to integrate the ODEs in Eq. (20) to obtain y(x) = y 1 (x), and then u(x) is obtained from Eq. (15) by inserting y(x). We depict the iterative boundary shape function method (BSFM) for finding u(x) in Eqs. (1)-(3): (i) Derive s 1 (x), s 2 (x), give y 1 (x 1 ), y 2 (x 1 ), z 0 , , and x = (x mx 1 )/N with N given.
(ii) For k = 0, 1, 2, . . . , applying RK4 to integrate the following ODEs with N 1 steps to if the stopping criterion r k := z k+1 -z k < is satisfied, then we stop the iteration; otherwise, for the next iteration go to step (ii). When y(x) = y 1 (x) is solved for, u(x) is obtained from Eq. (15) by inserting y(x).
For the details of the algorithm, we use the following four-point (x 1 = 0 < x 2 < x 3 < x 4 = 1) boundary conditions as a demonstrative example to find s 1 (x) and s 2 (x): Upon letting it follows that which can be arranged into (a + bx 2 ) + 2(a + bx 3 ) -6(a + bx 1 ) = 6, 2(a + bx 2 ) + 5(a + bx 3 ) -10(a + bx 4 ) = 0, and further changed to Thus, we can derive The variable transformation is then the transformed ODE is Starting from the given initial conditions y(x 1 ) = y (x 1 ) = 0, we can apply RK4 as shown in the Appendix to integrate the above ODE. With the initial guesses y(x 2 ) = y(x 3 ) = y(x 4 ) = 0, integrating with N 1 steps to x 2 , we can obtain the new value y(x 2 ); then N 2 steps to x 3 gives the new value y(x 3 ), and N steps to x 4 provides the new value y( ; z), we integrate it again. The process is continued, until the old values and the new values of y(x 2 ), y(x 3 ), y(x 4 ) are very close to satisfy the specified convergence criterion.
Remark 1 Liu [28] has pointed out the drawback of the shooting method, which assumes some unknown initial conditions u(0) and u (0) for Eq. (1) to convert the BVP into an IVP. It often requires many iterations to match the targets defined by the multipoint boundary conditions (2) and (3) through trial and error. In general, it is very difficult to find the exact values u(0) and u (0) for the nonlinear BVP with nonseparable and multipoint boundary conditions. Strictly speaking, the IVP used in the shooting method is not an exact one because its initial conditions are unknown. The current IVP being obtained exactly by using the variable transformation from u(x) to y(x) is different from the IVP that appeared in the shooting method in two aspects: the governing equation is Eq. (16) instead of Eq. (1), and the initial conditions y(0) and y (0) are given arbitrarily, not unknown values.

Numerical tests
In order to investigate the stability of the BSFM, the data b 1 and b 2 in Eqs. (2) and (3) are polluted by noise aŝ where s is the intensity of noise and R(i) are random numbers between [-1, 1]. Hence, sometimes we useb i , instead of b i , in the computations.
For the following parameters y 1 (1) = y 2 (1) = 0, z 0 = (0, 0) T , N = 200, and = 10 -10 , the iterative algorithm BSFM converges after 80 iterations as shown in Fig. 1(a). From Fig. 1(b), the numerical u(x) almost coincides with the exact one, with the maximum error (ME) being 7.16 × 10 -9 . Although the nonlinear nonseparable three-point BVP is difficult to be  treated by numerical methods, the accuracy of this problem is good, which is much better than that computed in [29] by about five orders. In order to test the influence of the noise on the numerical solution, in Table 1 we compare the the ME and iterations number (IN) for different noise levels. Upon comparing with the maximum value 17 of u(x), these MEs are acceptable.
We can observe in Table 1 that the IN is not sensitive to the noise level s. Instead of RK4, we have employed the fourth-order group preserving scheme [30] to integrate the resulting IVP with nonzero initial conditions y 1 (1) = 1, y 2 (1) = 0. With the same initial conditions, the BSFM converges with 81 iterations and the ME is 7.15 × 10 -9 . For the same parameter values, the fourth-order group preserving scheme converges with  Fig. 2(a). As observed in Fig. 2(b), the solution u(x) obtained almost coincides with the exact one, with the ME being 3.04 × 10 -10 . Although the nonseparable three-point BVP is more complex, the accuracy is much better than that computed in [29] by about eight orders.
The exact solution u(x) is still given in Eq. (23).
In order to test the influence of the noise on the numerical solution and the effect of large spatial range, in Table 2 we compare the ME and iterations number (IN) for different noise levels.
The fictitious time integration method (FTIM) was first developed by Liu and Atluri [31] to solve the following nonlinear algebraic equations: where u 1 , . . . , u n are unknown variables, and F i are given functions. After introducing the fictitious time τ , Eq. (28) is recast by Liu and Atluri [31] as a system of ODEs: They employed the forward Euler scheme to integrate the above ODEs, until the steady solution of u 1 , . . . , u n was obtained: According to Liu [29], FTIM for Eqs. (26) and (27) is given by where u i are the nodal values of u at the points x i = (i -1)/(n -1), and k = 1/(2 τ ) + 1. In Fig. 4(b), we compare the numerical solution u(x) with that computed by Liu [29] using FTIM, which are very close. In FTIM, we must choose some suitable values of ν i , τ , Liu [19] also applied the two-stage Lie-group shooting method (TSLGSM) to solve this problem, whose result is close to that obtained from FTIM and BSFM, and we do not plot it in Fig. 4(b). As shown in [19], one needs to solve four nonlinear algebraic equations derived from the Lie-group shooting equations to determine six unknown variables. The process of the TSLGSM is complex and is hard to be extended to an m-point BVP with m > 3.
For the following parameters y 1 (0) = -1, y 2 (0) = -1, z 0 = (0, 0, 0) T , N = 400, and = 10 -10 , the BSFM converges after 16 iterations as shown in Fig. 5(a). In Fig. 5(b), we com-  pare the numerical solution u(x) with the exact one in Eq. (33), whose ME is 2.63 × 10 -9 . The accuracy is very good, although the presented problem is nonlinear and is subjected to the nonseparable four-point boundary conditions. In order to test the influence of the noise on the numerical solution, in Table 3 we compare the the ME and iterations number (IN) for different noise levels. Upon comparing with the maximum value 16 of u(x), these MEs are acceptable.
In Table 5, we compare the ME and iterations number (IN) for different noise levels.
Remark 2 As shown in Examples 1, 3-7, where exact solutions are available, the accuracy obtained by BSFM is quite well, on the order of 10 -9 and 10 -10 . Because we have exactly transformed the multipoint BVP to the corresponding IVP, and integrated it by using RK4, the accuracy is on the order of ( x) 4 . For example, with x = 0.01, we have the error bound of 10 -8 . Therefore, we can estimate the ME by where M 0 is some positive constant. If the input data are noised by s, from Tables 1-3 and 5, we can observe that ME ≤ s.
We also compare Examples 3 and 5 with that obtained from Abd-Elhameed et al. [27]. For Example 3, Abd-Elhameed et al. [27] can obtain the exact solution with the ME being zero. On the other hand, BSFM led to the ME being 9.38 × 10 -8 . For Example 5, the accuracies obtained from BSFM and Abd-Elhameed et al. [27] are competitive. The wavelets collocation method with Chebyshev polynomials as the bases [27] led to residual algebraic equations to be solved to determine the coefficients. For some cases, the accuracy is very high. As shown by Example 3 in [27], the wavelets collocation method can also be applied to solve the singular nonlinear BVP as x(1x)u = 6 cosh x + 2 + xx 2 + sinh x sinh x -6u -2uu 2 , 0<x < 1, and with high accuracy as shown in Table 5 there. However, BSFM cannot treat this problem due to the left-hand side being zero at x = 0 and x = 1 when we apply RK4 to integrate the resultant IVP. The BSFM without needing to solve algebraic equations is an alternative candidate to solve the multipoint BVP efficiently.

Conclusions
In the paper, the boundary shape function was derived, which exactly and automatically satisfies the prescribed multipoint boundary conditions. It is of utmost importance that we can design the numerical method to exactly match the given multipoint boundary conditions. According to the new idea of boundary shape function, we have developed an iterative numerical algorithm used in solutions of the second-order nonlinear multipoint BVPs. The main contributions are introducing the boundary shape function, deriving a variable transformation, and then transforming the nonlinear BVP to the initial value problem (IVP). The resulting iterative algorithm resorting on the boundary shape function method (BSFM) is convergent very fast to a solution, and automatically satisfies the prescribed multipoint boundary conditions. Numerical examples confirmed that the BSFM is highly accurate and efficient. Even for some problems with large interval and subjected to the noise imposed on the boundary data, the presented new method is still workable to provide quite accurate solutions. The current idea of boundary shape function has been extended to multidimensional boundary value problems, for example, 2D problem [33,34] and 3D problem [35]. There, the higher-dimensional homogenization functions are constructed in a similar manner as the construction of BSF from the free function and with simple shape functions.