Forward-Backward Alternating Parallel Shooting Method for Multi-layer Boundary Value Problems

Multi-layer boundary value problems have received a great deal of attention in the past few years. This is due to the fact that they model many engineering applications. Examples of applications include uid ow though multi-layer porous media such as ground water and oil reservoirs. In this work, we present a new method for solving multi-layer boundary value problems. The method is based on an e cient adaption of the classical shooting method, where a boundary value problem is solved by means of solving a sequence of initial value problems. We propose, an alternating forward-backward shooting strategy that reduces computational cost. Illustration of the method is presented on application to uid ow through multi-layer porous media. The examples presented suggested that the method is reliable and accurate.


Introduction
Multi-layer boundary value problems (MLBVPs) have received a great deal of attention in the past few years. Many engineering problems are modelled by multi-layer boundary problems, where the system is described by a dierent dierential equation on dierent subsets of the physical domain. For example, MLBVPs are used as mathematical models for underground uid ow through porous channels with dierent characteristics. It is well known that the classical shooting method [1] is one of the most powerful methods used to solve single boundary value problems.
The classical shooting method is based on successively solving a sequence of initial value problems (IVPs) whose initial conditions are iteratively updated. The iterations continuously change (update) the assumed initial conditions until the solution of the initial-value problem satises the boundary conditions. Many modications to the classical shooting method either to improve the rate of convergence or to t a certain class of boundary value problems have been made. Examples of these modications include the modied simple shooting method [2], the multiple shooting method [3,4], and the parallel shooting method [5]. All of these modications deal with using the shooting method to solve a single boundary value problem over a nite domain. In this work, we consider an adaptation of the shooting method to solve multi-layer boundary value problems, which we refer to as multi-layer parallel shooting method (MLPSM). The advantage of this method is that it is suitable for parallel computing implementation. This kind of problem has been considered in [16], for a two layer problem, and in [17] for a multi-layer setting, both using a nite-dierence approach.
As an application of the suggested multilayer parallel shooting method, we will investigate the uid mechanics of multi-layer ows through porous media. The problem of the two layers ows was investigated by several authors [10]- [13]. In [10], the authors considered the interface region between two dierent porous media. They employed a special formulation of the shooting method to study the characteristics of the uid ow. In [11], a nite dierence approach was employed for the two channel problem. The multi-channel case was treated in [16,17] using nite dierences. In the present work, we extend the earlier works and formulate a multi-layer parallel shooting method for the solution of the multi-layer case using an alternating forward-backward shooting strategy that reduces the computational cost.
The paper is organised as follows. In section 2, the mathematical formulation of the problem is presented. The proposed method is presented in section 3. In section 4, application of the uid ow through multi-layer porous media is presented. Concluded remarks are given in section 5.

Mathematical Formulation of the Problem
We consider a second order multi-layer boundary value problem consisting of the following set of second order boundary value problems with outer boundary conditions and interior boundary conditions at the interface nodes The functions f i , i = 1, ..., n, are assumed to be C 2 . The interior nodes x i , i = 1, . . . , n − 1, divide the overall domain [x 0 , x n ], not necessarily uniformly, into subdomains [x i−1 , x i ], as shown in Fig. 1. The set of interior boundary conditions, equation (3), are imposed such that the overall solution y(x) dened over the whole domain [x 0 , x n ], is smooth across the nodes x i . This requirement is important in many physical problems such as in uid ow in multi-layer porous channels [11,12]. It is important to mention that the interior boundary conditions (3) make the problem more restrictive in the sense that the dierent solutions in neighbouring subdomains must agree smoothly at the nodes x i . Therefore, it is important that any numerical scheme should produce a numerical solution satisfying (3) as accurate as possible. As mentioned earlier, nite dierence methods can be, and in fact has been used [17], to solve equation (1) subject to (2) and (3). However, the accuracy in satisfying (3) dependents on the mesh size, whether it is a uniform mesh size across [x 0 , x n ] or dierent mesh sizes h i for each subdomain [x i−1 , x i ], as was done in [17].
This limitation of nite dierences method and the good performance of the classical shooting method prompt us to consider an adaptation of the classical shooting method to solve equations (1) subject to (2) and (3). The next section describes the proposed method. It is well-known that the shooting method is an accurate method for solving boundary value problems via solving a sequence of initial value problems. The classical shooting mechanism is performed in a forward way as shown in Fig. 2. Using, this classical mechanism for problem (1) subject to (2) and (3) requires the introduction of 2n − 1 parameters (1 at x 0 and 2 at each x i , i = 1, 2, . . . , n − 1). However, as we will see in the section the proposed strategy requires the introduction of n parameters, which reduces the computational complexity by half.

The Proposed Method
The proposed method to solve problem (1)(3) consists of a multi-layer shooting method which can be performed in parallel. The shooting strategy we propose is in a forward-backward alternating fashion as illustrated in Fig. 3. The derivations are the same for both n even and odd with a slight modication in the Jacobian matrix. We, therefore, consider in detail the case n even. According to the shooting strategy suggested in Fig. 3, we introduce n parameters λ i , i = 1, . . . , n, such that y (x n ) = λ n .
Then we solve, in parallel, the following set of initial-value problems: 1.
The parameters λ i , i = 1, . . . , n, are to be determined such that the interior boundary conditions (3) are satised to within a prescribed tolerance. Let y 1 (x; λ 1 ), y 2i (x; λ 2i , λ 2i+1 ), y 2i+1 (x; λ 2i , λ 2i+1 ), i = 1, ..., n/2−1, and y n (x; λ n ), be the solutions to (7)(10), respectively. Imposing conditions (3) at the oddly-indexed nodes x 2i−1 , i = 1, 2, . . . , n/2 − 1, (eq. (3)) is already satised at x 2i ), we obtain the set of equations The above set of equations can be regarded as a nonlinear homogeneous system for the unknown parameters λ i . Introducing the notation Λ = (λ 1 , λ 2 , . . . , λ n ), it can be written in vector form as .., F n (Λ)] T : R n −→ R n , is a vector valued function of the parameter vector Λ and F i (·), i = 1, . . . , n, are scalar functions given by equations (11)(16). A fast and ecient method for solving the above system is the well-known multidimensional Newton's method. Let 2 , . . . , λ (k) n ) T be the values of the parameters at the k th iteration. The classical multidimensional Newton's method calculates new values Λ (k+1) according to where is the inverse of the n × n Jacobian matrix J evaluated at Λ (k) . The entries of the Jacobian matrix are given by It is worth mentioning that there are modied versions of Newton's method which converge cubically. In the one-dimensional case, cubically convergent modied Newton's schemes were derived in [6,8,9] and references therein. A generalization to multivariate case was later given in [7]. The scheme given in [7] is as follows The Jacobian matrix J turn out to be pentadiagonal and has the form: for i = 3, 5, 7, . . . , n − 3, and J n−1,: = 0 n−3 where the notation J k,: stands for the kth row of J and 0 k stands for a row of k zeros. The calculation of the entries of the Jacobian matrix requires the calculation of with the appropriate x j and λ k for a given y i . Specically, we require for i = 1, 2, . . . , n/2 − 1 and k = 2i, 2i + 1, and The quantities in (21) and (22) are obtained from the solutions of an other set of initial value problems as explained below. From equation (1), dierentiating with respect to λ, we get ∂y ∂λ = ∂f ∂y ∂y ∂λ + ∂f ∂y ∂y ∂λ .
Let z j,k = ∂y j ∂λ k . Assuming we can interchange the order of dierentiation, we nd that z j,k satises Therefore, we have the following I.V.Ps.
Once the solutions z j,k of the above I.V.Ps (23) with (24)(29) are obtained, the dierent Jacobian entries (21) and (22) are given by The proposed algorithm is summarized as follows.
It is important to mention here that at each step solving the IVPs (7)(10) can be done in parallel. Similarly, solving the IVPs (23) with (24)(29) can also be done in parallel. This is an important advantage of the proposed method over nite dierences methods.

Numerical Examples
In this section we consider two cooked examples with exact solutions in order to test the accuracy of the proposed method. For the numerical simulations, we adopted the stopping criteria F ≤ ε = 10 −6 . Example 4.1. As a rst example, we consider the following second-order BVP with the boundary conditions y(0) = 2 and y(π) = 2.
The exact solution is We applied the proposed algorithm to this example and convergence to within the prescribed accuracy was achieved after 6 iterations. Fig. 5, displays the results for the 6 iterations. Table 2 displays the values of y(x ± i ) and y (x ± i ) versus iterations which can be seen to converge to those of the exact solution.   The previous two example demonstrate that the proposed method is a reliable, fast and accurate for the numerical simulation of multi-layer problems in (1)- (3). As an application of the proposed method to a physical problem, we shall, in the next section, use it to numerically resolve the velocity prole of uid ow through multi-layer porous media.

Application to Fluid Flow Through Multi-layer Porous Media
In this section, we present an application of the algorithm presented in Section 3 to the problem of determining the velocity prole of uid ow through multi-layer porous media. The media is composed of a number of porous layers. The top and the bottom layers are bounded above and below, respectively, by solid walls, see Fig. 6 which depicts a multi-layer porous media conguration. In each porous layer, the ow is assumed to be governed by the Darcy-Forchheimer-Brinkman (DFB) model [12] d 2 u dy 2 = Re C + where u(y), a ≤ y ≤ b, is the velocity of the uid, and the various physical parameters in Eq. (33) are dened as follows.
1. Re = ρU ∞ L/µ is the Reynolds number, ρ is the uid density, U ∞ is the free stream characteristic velocity, µ is the uid viscosity, and L is the channel characteristic length. 2. k is the permeability of the porous channel. 3. C d is the form drag coecient. 4. C < 0 is a dimensionless pressure gradient.
When C d = 0, we have what is known as the Darcy-Lapwook-Brinkman (DLB) [12] model It is worth mentioning that there have been a lot of work on uid ow through porous media [10]- [17]. The book by Chen et al. [15] presents an excellent account on computational methods for multiphase ows in porous media. The method used in [17] is based on nite dierence approximation of the derivatives.
In this work, we aim to show that the proposed method is suitable for solving for the velocity prole of uid ow through any number of porous-layer conguration. We will use the two examples considered in [17] of the three-and ve-layer congurations. In our simulations, we x the following parameter values. The Reynolds number Re = 10, C = −10 and C d = 0.55, and without loss of generality, we assume that overall domain is −1 ≤ y ≤ 1. Since the top (at y = 1) and bottom (at y = −1) layers are bounded by solid impermeable walls, we assume a non-slip condition at y = ±1, that is u(±1) = 0. At the interface nodes x i , we assume that the ow and shear stress are continuous, that is, . Therefore, here, what distinguish between the dierent porous layers is the permeability parameter k.
We consider the cases of a three-porous media conguration and ve-layer porous media conguration, as considered in [17], where the ow in all three layers is governed by the DFB model (Eq. (33)). For the three-layer congurations, we use the following permeability values: k b = 0.001, k m = 0.01, k t = 100, where k b , k m , and k t are the permeabilities of the bottom, middle, and top layers, respectively. The interface nodes are x 1 = −1/2 and x 2 = 1/2. For the ve-layer conguration, we use the following permeability values: k 1 = 0.1, k 2 = 1, k 3 = 0.001, k 4 = 10, k 5 = 0.001, where k 1 is for the bottom layer (−1 ≤ y ≤ x 1 ) and k 5 for the top layer (x 4 ≤ y ≤ 1). The interface nodes for this case are x 1 = −0.6, x 2 = 0, x 3 = 0.4, x 4 = 0.6. The uid velocity proles for both cases are displayed in Fig. 7 which show an agreement with the results found in [17].

Conclusion
In this paper we have presented a parallel shooting method to solve multi-layer boundary values problems. The method relies on the adaption of the classical shooting method for single boundary value problems to multi-layer boundary value problems. We have adopted an alternating forward-backward shooting strategy which decreases the size of the problem, compared to a classical forward shooting strategy. A good advantage of this method is that it can be implemented in parallel, unlike the nite-dierence method, hence reducing the computation time, especially when the number of layers exceeds 3. A real application of the method is in the resolution of the velocity prole of uid ow through multi-layer porous media. The numerical simulations suggest that the proposed method is reliable and accurate.