A Nonlinear Shooting Method and Its Application to Nonlinear Rayleigh-Bénard Convection

The simple shooting method is revisited in order to solve nonlinear two-point BVP numerically. The BVP of the type X(z) = f(z,X(z)), X(z) ∈ Rn, n > 1, z ∈ [0, 1], is considered where m components of X(z) are known at one of the boundaries and (n − m) components of X(z) are specified at the other boundary. The map f is assumed to be smooth and satisfies the Lipschitz condition. The two-point BVP is transformed into a system of nonlinear algebraic equations in several variables which, is solved numerically using the Newton method. Unlike the one-dimensional case, the Newton method does not always have quadratic convergence in general. However, we prove that the rate of convergence of the Newton iterative scheme associated with the BVPs of present type is at least quadratic. This indeed justifies and generalizes the shooting method of Ha (2001) to the BVPs arising in the higher order nonlinear ODEs.With at least quadratic convergence of Newton’s method, an explicit application in solving nonlinear Rayleigh-Bénard convection in a horizontal fluid layer heated from the below is discussed where rapid convergence in nonlinear shooting essentially plays an important role.


Introduction
Let X : [0, 1] → R  ,  > 1, be a vector valued function defined by X() = ( 1 (),  2 (), . . .,   ()) where each map   : [0, 1] → R, 1 ≤  ≤ , is smooth over the interval (0, 1).We consider the following vector differential equation satisfied by X(): Throughout, the map  : [0, 1] × R  → R  is assumed to be smooth and satisfying the Lipschitz condition on a closed rectangle  := [0, 1] ×  ⊂ [0, 1] × R  with a Lipschitz constant  > 0 s.t.     (, X ()) −  (, Y ())     ≤  ‖X () − Y ()‖ , (2) for all X, Y ∈  and  ∈ [0, 1]; furthermore, the components of X() are assumed to satisfy  initial conditions at first boundary given by and (−)-conditions at the other boundary which are given by ( (1) , . . .,  (−) ) (1) = ( 1 , . . .,  − ) ∈ R − , where  is a permutation of the symmetric group   .Equations ( 1)-( 4) lead to a two-point boundary value problem whose solution is not known a priori at either of the boundaries.Such BVPs are a common object of study in mathematics, physics, engineering, stochastic analysis, and optimization.In general, it is not possible to solve these BVPs analytically, and one needs to look for their numerical solutions in order to unfold the inherent dynamics.To do so, the vector differential equation may be reduced to a set of nonlinear algebraic equations by approximating the solution X() with a (finite) Galerkin expansion in terms of a suitably chosen set of orthogonal functions already satisfying the conditions of the BVP in hand.The resulting set of the nonlinear algebraic equations is solved numerically for the unknown Galerkin coefficients using iterative methods such as the Newton-Raphson scheme.Despite its general applicability, the Galerkin method becomes handy as the number of Galerkin coefficients becomes large.
A delicate, easy, and direct numerical approach to solve two-point BVPs (see ( 1)-( 4)) is using conventional shooting methods which are quite general and applicable to a wide variety of ODEs.Depending upon the type of two-point BVP, people have developed shooting methods (see for details [1][2][3][4][5][6][7][8][9]).Some of these works utilize Adomian decomposition method to obtain a series solution of the initial value problem involved.
Shooting method appears to be a straightforward way to solve nonlinear BVPs but many difficulties arise while actually carrying out numerical calculations using it.For instance, it requires an intelligent guess for the missing initial conditions.Here, the problem of convergence of the Newton iterates can be very troublesome, the major problem being the location of the root within the radius of convergence and that too is not known a priori.The degree of difficulty increases as the dimension  increases and shooting method practically becomes useless.Due to this drawback, the conventional simple shooting methods are not preferred to obtain numerical solution of nonlinear BVPs.
To overcome difficulty of convergence of the simple shooting methods, more advanced parallel shooting is used where the interval of the solution is subdivided into small parts and the BVP is solved using simple shooting over each part with a continuity condition at the nodal points.The parallel shooting is advantageous over simple shooting as the former is capable of handling a larger class of two-point BVPs.However, it may require a huge number of parameter updates at each iteration causing large computation times [1,4,5].Ha [7] has discussed the simple shooting method for nonlinear two-point boundary value problems and observed the rapid convergence in his numerical experiments.In the last decade Liu [10] has developed a noble approach of Lie-group preserving schemes for integrating the nonlinear system represented by (1), using homogeneous coordinates in the space of Minköwski type.Based on this approach, Liu [11] has developed an efficient Lie-group shooting method for second order nonlinear BVPs and successfully applied it to solve the boundary-layer flow problems for power law fluids.An attempt towards tackling higher order nonlinear twopoint BVPs involving a modification of Adomian decomposition has been made by Wazwaz [6].
In spite of the wide literature and convergence analysis, it seems that the use of the simple shooting method has been neglected over the years by researchers to obtain numerical solutions of the higher order nonlinear two-point BVPs.The goal of the present discussion is to address the importance and capability of the simple shooting method to easily handle numerical solutions of the higher order two-point BVPs arising in other areas of research.Moreover, we prove that under certain conditions on the BVP in hand, the Newton iterates derived from the shooting method always have at least quadratic convergence.This result then ensures the rapid convergence in the nonlinear shooting.It also gives a possible explanation for the behavior of convergence on the initial guess as observed by Ha [7].In establishing these results, we only need the following well-known results in the study of a system of ODEs (see Birkhoff and Rota [12]).
Theorem 1 (existence).Let (, X()) be continuous and satisfying Lipschitz condition as in (2) on interval | − | ≤  for all X ∈ .Then for any constant vector c, the differential equation (1) has a solution on | − | ≤  with initial condition X() = c.Theorem 2 (uniqueness).If (, X()) satisfies Lipschitz condition on , then there is at most one solution of (1) with initial condition X() = c.Theorem 3 (continuity).If X() and Y() are any two solutions of (1) where (, X()) is continuous and satisfies Lipschitz condition of (2), then The next theorem is the statement of the Newton-Kantorovich theorem that establishes convergence of the Newton-Raphson method in the higher dimensions.For a quick review, the reader is referred to the work of Ortega [13], Tapia [14], Rall [15], and Gragg and Tapia [16].
The rest of the paper is organized as follows.The simple shooting method is discussed in Section 2 where the convergence Theorem 6 is proved.Based on Theorem 6, the underlying numerical integration scheme is explained to solve the underlying BVP in Section 2.1 and is successfully implemented in Section 3 to verify its validity via solving a nonlinear two-point BVP of order eighteen arising in the hydrodynamic studies.
Theorem 6.The hypothesis of Theorem 5 is satisfied by the sequence {(c  )} of functions defined by on the closed rectangle  ⊂ R − containing point c in its interior where The sequence {c  } is defined recursively by the following: Proof.Since the Jacobian matrix X (, c) :=  X/c(, c) satisfies the linear matrix differential equation where f is the vector with ( − ) components of  depending upon the permutation , it follows from the existence and uniqueness theorem for the system of ODEs that X ( . ( Note that  > 0 as (/X)(, X(, c)) is nonsingular and hence not zero identically.Consequently, with an application of Theorem 3, we arrive at the following: showing that {  (c  )},  = 0, 1, . .., is a constant sequence.It follows that   satisfies the Lipschitz condition on  0 ∋ c with any real number  > 0 as a Lipschitz constant.We choose  > 0 small enough such that for ‖  (c) −1 ‖ ≤  ̸ = 0 and ‖  (c) −1 (c)‖ ≤  ̸ = 0, the following strict inequality holds: This completes the proof.Proof.The convergence of the Newton iterates follows from the Newton-Kantorovich Theorems 5 and 6 with  =  and convergence of the solution X(, c) to a solution X(, c * ) as c converges to c * , follows from the corollary to Theorem 3.
Remark 8.It turns out that the function defined by the sequence {(c  )} as aforementioned satisfies the Lipschitz condition with a Lipschitz constant   and hence is dominated by the sequence {c  }.To establish this, let Since  is given to satisfy the hypothesis of Theorem 3 with Lipschitz constant , it follows that for a fixed  and each  =  −1 , the following holds for the two solutions X(, c +1 ) and X(, c  ) of ( 8): and correspondingly, the following is satisfied: From ( 18) we obtain the desired assertion for the function (c  ),  = 0, 1 . .., as follows: where we have made use of the fact that the sequence {  (c  )} is constant.It also follows from the estimate (19 Remark 9.It is also known that if   () satisfies Lipschitz condition then the underlying Newton iterates have at least quadratic convergence to a root of () = 0 (see Rall [15]).
2.1.Numerical Scheme.The associated numerical scheme involves computation of numerical solution X(, c * ) of the BVP defined by (1)-(4) which at ( + 1)th iteration involves simultaneous numerical integration of (8) from  = 0 to  = 1 to obtain (c  ) and yet another numerical integration of the IVPs defined by the linear system in (10) for each   from  = 0 to  = 1 to obtain X  (1, c  ).Since the underlying numerical integrations are to be carried out on interval [0, 1], we define another equivalent vector differential equation as follows: for each  = 1, . . ., ( − ) . (20) These are ( − ) IVPs, each of order 2 and can be integrated numerically using the fourth order Runge-Kutta method over the interval [0, 1], to obtain ( X X/  ) (1,   ) for each   =   1 , . . .,   − s.t.c  = (  1 , . . .,   − ),  = 0, 1, . . .at each Newton iteration from which it is easy to extract X(1, c  ) and the Jacobian matrix X (1, c  ) given by Now each Newton iteration is performed to obtain c +1 by solving the following system of ( − ) linear algebraic equations in  +1 1 , . . .,  +1 − given by To avoid any numerical singularity, the linear system in ( 22) is solved using an iterative method such as the Gausselimination method which is easy to implement in computer programs with a control over accuracy.We have developed a numerical code in MATLAB for simultaneously solving (20) and ( 22), recursively, till the scheme converges and meets the desired tolerance.To check the validity of the method, we consider the following example.

Example 10. Consider the nonlinear BVP (
, with the boundary conditions  1 (1) = 2,  1 (2) = 5/2.The exact solution on [1,2] is  1 () =  + 1/.We choose step size ℎ = 0.05 to integrate the underlying IVP with 20 points of evaluation using the fourth order Runge-Kutta method with a tolerance of 10 −5 .Table 1 gives the description of the convergence analysis of the solution for several values of the trial initial guess for c 0 = (2,  0 1 )  .It is clear from the table that the initial guess  0 1 = 0 is appropriate for which the numerical scheme converges at the very first iteration to zero within an error of the order of 10 −12 .If we take the initial guess  0 1 away from 0 then a little larger number of iterations is required to achieve convergence of the numerical scheme within an error at most of the order of 10 −6 which corresponds to the value  0 1 = 10.
This example has been considered by Ha [7] where no suitable solution was found via his numerical experiment with the same BVP for  0 1 ≥ 1.We have not encountered this situation except for possibly  0 1 > 10 when the initial guess does not fall into the the disc of radius  * * (see Theorem 5).Also for  0 1 = 0.6975, Ha achieves convergence after 9982 iterations; however, we have found that it takes only 11 iterations to obtain the solution correct within an error of at most 5.3561 − 008.Note that the system of ODEs considered here is locally Lipschitz with Lipschitz constant   ≥ 40 on the closed rectangular strip of width 1 defined by  := {( 1 ,  2 ) ∈ R 2 | ‖  1 − 2 ‖ ≤ 0.5}; therefore, by Theorem 6, it follows that the present numerical scheme will always converge for each initial guess taken from  on the line  1 = 2, but with large  0 1 , the speed of convergence will slow down as indicated in Table 1.
The numerical solutions  1 and  2 are drawn as ⋇ in Figure 1 and the actual solution corresponds to the solid line in the plots.It is clear from Figure 1 that the numerical solution converges to the actual solution for all initial guesses given in Table 1.
We remark that our numerical scheme converges with all initial guesses for BVPs originating from locally Lipschitz systems of ODEs.In such systems, we need to evaluate the Jacobian only once for all iterations which facilitates considerable amount of lowering down the cost of numerical calculations.

The Rayleigh-Bénard Convection
We verify the validity of Theorem 6 by implementing the shooting method to solve a typical higher order nonlinear system of ODEs arising in the research studies involving nonlinear normal-mode stability analysis of the Navier-Stokes equations for the standard Rayleigh-Bénard convection occurring in an initially quiescent fluid layer resting between two parallel rigid horizontal planes  = 0 and  =  (see Figure 2).
The fluid layer is maintained at different constant temperatures  1 and  2 at its two boundaries, respectively.The fluid conduction state is known to prevail till a critical temperature gradient across the fluid layer is exceeded when the instability appears in the form of parallel rolls, squares, or hexagon pattern (for details see Chandrasekhar [17] and Drazin and Reid [18]).At any time  and any point (, , ) in the fluid layer, u is the fluid velocity,  is the fluid temperature,  is the fluid pressure, and  is the fluid density, and then the flow is described by the following equations [17]: where  0 is the fluid density at some reference fluid temperature   , ] is the fluid viscosity, g = (0, 0, −) is the acceleration due to gravity,   is the diffusivity constant,   is the specific heat of the fluid at constant volume, and  := (/)| =  .
The system (23) admits a closed form steady state of the form where the subscript  denotes the equilibrium state.We discuss the stability of the solution described by (24) by superimposing perturbations on it in the following form: where we have utilized the thickness  of the fluid layer as the length scale,   / as the velocity scale, and  1 −  2 as the temperature scale.Assuming the perturbations given by (25) as the solution of the governing equations, the perturbations satisfy the following set of partial differential equations where 0 ≤  ≤ 1,  > 0; the dimensionless numbers Pr := ]/  and Ra := ( 1 −  2 ) 3 /  ] are called the Prandtl number and the Rayleigh number, respectively.The boundary conditions for the velocity field and the temperature field are given by Rigid-Boundaries: u| =0,1 = 0, Free-Boundaries: u| =0,1 = 0, To solve the nonlinear system of (26), we consider the solution in the following form: and we obtain the nonlinear system of ODEs as in ( 30) and (34) for a fixed value of (ℓ, ) = (3, 0) for rolls solutions and (ℓ, ) = (1, 1) for square patterns, respectively, for the case of rigid boundaries.For a fixed (ℓ, )-mode such that ℓ+ = , we assume the following relationship between the horizontal wave numbers  1 and  2 : Consequently, for  = 0 the choice for  2 is arbitrary while for ℓ = 0,  1 =  2 .Therefore, for a choice  1 =  2 , the two problems for (ℓ, 0) and (0, ℓ) modes are essentially the same.For ℓ ̸ = 0 and  ̸ = 0, the problem for (, ℓ)-mode is equivalent to the corresponding problem for (ℓ, )-mode under the change  1 → ( √ (ℓ + ) 2 − ℓ 2 /) 1 .The solution corresponding to  = 0 is called the rolls solution, and for ℓ =  the solution is called the square pattern.
The BVP defined by ( 30) is solved numerically using the proposed numerical scheme (see ( 20) and ( 22)) with the help of MATLAB programming.The time lapsed for obtaining the missing boundary conditions in each run is found to be between 2 and 5 seconds, which shows that the numerical convergence is quite rapid.Only a few Newton iterations are required for convergence of the present numerical scheme with an initial guess of c = 0 ∈  = {x ∈ R The corresponding normalized solution curves satisfying the given boundary conditions are illustrated in Figure 3 which straightway verifies the correctness of the present numerical scheme.Here in Figure 3,  11 ,  21 , and  31 are the fluid velocity profiles corresponding to the linear mode, second order nonlinear mode, and the third order nonlinear mode, respectively.Similarly,  13 ,  23 , and  33 are the corresponding profiles for the temperature field perturbations inside the fluid.

Square Pattern.
In this subsection, we obtain the stationary square pattern in the Rayleigh-Bénard convection in the horizontal fluid layer with rigid boundaries which corresponds to ℓ = 1 and  = 1.Here: the underlying nonlinear dynamical system is given by the following twelve ODEs: for  = 0, 1 where   := ( + 1) and the twelve boundary conditions are given by for each  = 0, 1.The fixed parametric values are taken as the following ( 0 ,  1 , Pr,  1 , Ra) = (85.80,0, 10, 1.603, 15874) .In this case, the missing boundary conditions for  = 0 using an initial guess of c = 0 in the rectangle  defined by  := {x ∈ R 6 |‖ x − c ‖ < 2200} are found to be the following: Here the Newton iterates are also found to converge rapidly.The nonlinear square pattern exhibited by the Rayleigh-Bénard convection as obtained using the previous numerical calculations is shown in Figures 4 and 5.The division of the fluid layer in the form of square pattern of the streamlines and the isotherms in the horizontal plane is clear from Figure 4.

Concluding Remarks
The simple shooting method for solving nonlinear higher order BVPs is revisited to establish the rapid convergence of the underlying Newton iterates.We prove at least quadratic convergence of the underlying Newton iterates arising in the numerical solution of the two-point BVPs by the conventional shooting method.Only the smoothness and nonsingularity of the underlying dynamical system and the Lipschitz condition, are assumed.The Newton method in the higher dimensions usually has linear convergence.Based on the present analysis, we have proposed a numerical scheme for the numerical solution of the higher order two-point BVPs  and used it successfully to solve hydrodynamic problem of the Rayleigh-Bénard convection in a horizontal fluid layer heated from the below.Convergence of the proposed numerical scheme is found to be rapid.Owing to the simplicity and high accuracy of the proposed shooting method, we recommend its use for handling a wide variety of nonlinear BVPs arising in different disciplines.
Finally, we would like to remark that the rapid convergence of the simple shooting method (Theorem 6) in the BVPs arising from the nonsingular Lipschitz systems of ODEs establishes the importance of nonlinear shooting specially while dealing with the research studies involving numerical solutions of higher order systems.The present analysis does not consider the two-point BVPs with mostgeneral boundary conditions.This is another direction for the future studies in the present context.

Figure 1 :
Figure 1: The numerical solution ⋇ versus the actual solution of the BVP of Example 10.

Figure 2 :
Figure 2: The schematic diagram for Rayleigh-Bénard convection in an infinite horizontal fluid layer confined between the planes  = 0 and  =  which are maintained at constant temperatures  1 and  2 , respectively.

Table 1 :
Dependence of convergence of the numerical scheme on the initial guess.