About this Journal Submit a Manuscript Table of Contents
ISRN Mathematical Physics
Volume 2013 (2013), Article ID 650208, 9 pages
Research Article

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

Department of Mathematics, Guru Nanak Dev University, Amritsar 143005, India

Received 5 June 2013; Accepted 3 July 2013

Academic Editors: D. Dürr and A. Qadir

Copyright © 2013 Jitender Singh. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


The simple shooting method is revisited in order to solve nonlinear two-point BVP numerically. The BVP of the type is considered where components of are known at one of the boundaries and components of are specified at the other boundary. The map 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.

1. Introduction

Let ,  , be a vector valued function defined by , where each map , , is smooth over the interval . We consider the following vector differential equation satisfied by : Throughout, the map is assumed to be smooth and satisfying the Lipschitz condition on a closed rectangle with a Lipschitz constant s.t. for all   and ; furthermore, the components of are assumed to satisfy initial conditions at first boundary given by and -conditions at the other boundary which are given by 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 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 [19]). 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 two-point 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 be continuous and satisfying Lipschitz condition as in (2) on interval for all . Then for any constant vector , the differential equation (1) has a solution on with initial condition .

Theorem 2 (uniqueness). If satisfies Lipschitz condition on , then there is at most one solution of (1) with initial condition .

Theorem 3 (continuity). If and are any two solutions of (1) where is continuous and satisfies Lipschitz condition of (2), then

Corollary 4. Let be a solution of (1) with and hypothesis of Theorem 3 being satisfied. Let be defined on , . Then is a continuous function of and . Further, if , then uniformly for .

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].

Theorem 5 (Newton-Kantorovich theorem). Let and be Banach spaces and such that on an open convex set , the map is Fréchèt differentiable and satisfies For some assume that is defined on all of and that , where and . Define , , and . Then the Newton iterates are well-defined and lie in and converge to a solution of which is unique in . Moreover if , then the order of convergence is at least quadratic.

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.

2. The Simple Shooting Method

We introduce a vector , , independent of the variable s.t. and the smooth function satisfies (1)–(4); that is

Assume that and are defined on the closed rectangle and are continuously differentiable on such that is nonsingular there. Observe that at point ,  ,

If then we see that , where is an element of the standard basis of the vector space . It follows that for each component of , the derivative satisfies initial value problems (IVPs) one for each . Under these considerations, we have the following main theorem.

Theorem 6. The hypothesis of Theorem 5 is satisfied by the sequence of functions defined by
on the closed rectangle containing point in its interior where s.t. ;      , ;     . The sequence is defined recursively by the following:

Proof. Since the Jacobian matrix satisfies the linear matrix differential equation
where is the vector with components of depending upon the permutation , it follows from the existence and uniqueness theorem for the system of ODEs that is defined and invertible in an open rectangle containing the point . Also, since is a solution of a linear vector differential equation in satisfies the Lipschitz condition with a Lipschitz constant
Note that as is nonsingular and hence not zero identically. Consequently, with an application of Theorem 3, we arrive at the following:
showing that , is a constant sequence. It follows that satisfies the Lipschitz condition on with any real number as a Lipschitz constant. We choose small enough such that for and , the following strict inequality holds: This completes the proof.

Corollary 7. The sequence of Newton iterates , defined by (12) of the Theorem 6, is well-defined and lies in the set and converges to unique solution of the equation on . Furthermore, as .

Proof. The convergence of the Newton iterates follows from the Newton-Kantorovich Theorems 5 and 6 with and convergence of the solution to a solution as converges to , follows from the corollary to Theorem 3.

Remark 8. It turns out that the function defined by the sequence as aforementioned satisfies the Lipschitz condition with a Lipschitz constant and hence is dominated by the sequence . To establish this, let ,  ,  , be a partition of . Since is given to satisfy the hypothesis of Theorem 3 with Lipschitz constant , it follows that for a fixed and each , the following holds for the two solutions and of (8):
and correspondingly, the following is satisfied:
From (18) we obtain the desired assertion for the function , as follows:
where we have made use of the fact that the sequence is constant. It also follows from the estimate (19) that for .

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 (see Rall [15]).

2.1. Numerical Scheme

The associated numerical scheme involves computation of numerical solution of the BVP defined by (1)–(4) which at th iteration involves simultaneous numerical integration of (8) from to to obtain and yet another numerical integration of the IVPs defined by the linear system in (10) for each from to to obtain . Since the underlying numerical integrations are to be carried out on interval , we define another equivalent vector differential equation as follows:

These are IVPs, each of order and can be integrated numerically using the fourth order Runge-Kutta method over the interval , to obtain for each s.t. ,   at each Newton iteration from which it is easy to extract and the Jacobian matrix given by

Now each Newton iteration is performed to obtain by solving the following system of linear algebraic equations in given by

To avoid any numerical singularity, the linear system in (22) is solved using an iterative method such as the Gauss-elimination 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 ,  . The exact solution on is . We choose step size to integrate the underlying IVP with 20 points of evaluation using the fourth order Runge-Kutta method with a tolerance of . Table 1 gives the description of the convergence analysis of the solution for several values of the trial initial guess for . It is clear from the table that the initial guess is appropriate for which the numerical scheme converges at the very first iteration to zero within an error of the order of . If we take the initial guess away from 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 which corresponds to the value .

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

This example has been considered by Ha [7] where no suitable solution was found via his numerical experiment with the same BVP for . We have not encountered this situation except for possibly when the initial guess does not fall into the the disc of radius (see Theorem 5). Also for , 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 . Note that the system of ODEs considered here is locally Lipschitz with Lipschitz constant on the closed rectangular strip of width defined by ; therefore, by Theorem 6, it follows that the present numerical scheme will always converge for each initial guess taken from on the line , but with large , the speed of convergence will slow down as indicated in Table 1.

The numerical solutions and 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.

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

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.

3. 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 and (see Figure 2).

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

The fluid layer is maintained at different constant temperatures and 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, 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 is the fluid density at some reference fluid temperature ,   is the fluid viscosity, 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 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 ; the dimensionless numbers and are called the Prandtl number and the Rayleigh number, respectively. The boundary conditions for the velocity field and the temperature field are given by

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 for rolls solutions and 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 and : Consequently, for the choice for is arbitrary while for , . Therefore, for a choice , the two problems for and modes are essentially the same. For and , the problem for mode is equivalent to the corresponding problem for mode under the change . The solution corresponding to is called the rolls solution, and for the solution is called the square pattern.

3.1. Rolls Solution

The underlying hydrodynamic system consists of ODEs for each normal mode and . The linearized system corresponds to the lowest normal mode . Higher modes include the nonlinear effects. For a typical value the dynamical system is of order eighteen and is described by the following set of ODEs: where , , , for all ; ; for the terms under summation sign in the R.H.S. of (30) are absent. This system satisfies nine conditions at each of the rigid boundaries

The vector differential equation (30) involves unknown constants , and for a given value of . As a check we take .

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 . The missing boundary conditions for are found to be given by 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, , , and are the fluid velocity profiles corresponding to the linear mode, second order nonlinear mode, and the third order nonlinear mode, respectively. Similarly, , , and are the corresponding profiles for the temperature field perturbations inside the fluid.

Figure 3: A third order nonlinear solution (normalized) of (30) for parametric values .
3.2. 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 and . Here: the underlying nonlinear dynamical system is given by the following twelve ODEs: for where ; , and the twelve boundary conditions are given by for each . The fixed parametric values are taken as the following

In this case, the missing boundary conditions for using an initial guess of in the rectangle defined by 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.

Figure 4: Square pattern in the Rayleigh-Bénard convection for .
Figure 5: Square pattern (side view in plane) in the Rayleigh-Bénard convection for .

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 most-general boundary conditions. This is another direction for the future studies in the present context.


  1. H. B. Keller, Numerical Solution of Two Point Boundary Value problems, SIAM, Philadelphia, Pa, USA, 1976. View at Zentralblatt MATH · View at MathSciNet
  2. A. Granas, R. B. Guenther, and J. W. Lee, “The shooting method for the numerical solution of a class of nonlinear boundary value problems,” SIAM Journal on Numerical Analysis, vol. 16, no. 5, pp. 828–836, 1979. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  3. R. M. M. Mattheij and G. W. M. Staarink, “On optimal shooting intervals,” Mathematics of Computation, vol. 42, no. 165, pp. 25–40, 1984. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  4. J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, vol. 12 of Texts in Applied Mathematics, Springer, New York, NY, USA, 2nd edition, 1993. View at MathSciNet
  5. M. E. Kramer and R. M. M. Mattheij, “Application of global methods in parallel shooting,” SIAM Journal on Numerical Analysis, vol. 30, no. 6, pp. 1723–1739, 1993. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  6. A.-M. Wazwaz, “Approximate solutions to boundary value problems of higher order by the modified decomposition method,” Computers & Mathematics with Applications, vol. 40, no. 6-7, pp. 679–691, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  7. S. N. Ha, “A nonlinear shooting method for two-point boundary value problems,” Computers & Mathematics with Applications, vol. 42, no. 10-11, pp. 1411–1420, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  8. A. M. Wazwaz, “A reliable algorithm for obtaining positive solutions for nonlinear boundary value problems,” Computers & Mathematics with Applications, vol. 41, no. 10-11, pp. 1237–1244, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  9. B. S. Attili and M. I. Syam, “Efficient shooting method for solving two point boundary value problems,” Chaos, Solitons and Fractals, vol. 35, no. 5, pp. 895–903, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  10. C.-S. Liu, “Cone of non-linear dynamical system and group preserving schemes,” International Journal of Non-Linear Mechanics, vol. 36, no. 7, pp. 1047–1068, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  11. C. S. Liu, “The Lie-group shooting method for boundary-layer problms with suction/injection/reverse flow conditions for power-law fluids,” International Journal of Non-Linear Mechanics, vol. 46, pp. 1001–1008, 2011.
  12. G. Birkhoff and G.-C. Rota, Ordinary Differential Equations, John Wiley & Sons, New York, NY, USA, 3rd edition, 1978. View at MathSciNet
  13. J. M. Ortega, “The Newton-Kantorovich Theorem,” The American Mathematical Monthly, vol. 75, pp. 658–660, 1968. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  14. R. A. Tapia, “Classroom notes: the Kantorovich theorem for Newton's method,” The American Mathematical Monthly, vol. 78, no. 4, pp. 389–392, 1971. View at Publisher · View at Google Scholar · View at MathSciNet
  15. L. B. Rall, “A note on the convergence of Newton's method,” SIAM Journal on Numerical Analysis, vol. 11, pp. 34–36, 1974. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  16. W. B. Gragg and R. A. Tapia, “Optimal error bounds for the Newton-Kantorovich theorem,” SIAM Journal on Numerical Analysis, vol. 11, pp. 10–13, 1974. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  17. S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, The International Series of Monographs on Physics, Clarendon Press, Oxford, UK, 1961. View at MathSciNet
  18. P. G. Drazin and W. H. Reid, Hydrodynamic Stability, Cambridge Mathematical Library, Cambridge University Press, Cambridge, UK, 2nd edition, 2004. View at Publisher · View at Google Scholar · View at MathSciNet