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

Abstract

The simple shooting method is revisited in order to solve nonlinear two-point BVP numerically. The BVP of the type 𝐗^′(𝑧)=𝑓(𝑧,𝐗(𝑧)),𝐗(𝑧)∈ℝ^𝑛,𝑛>1,𝑧∈[0,1], 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 𝐗∶[0,1]→ℝ^𝑛,  𝑛>1, be a vector valued function defined by 𝐗(𝑧)=(𝑋₁(𝑧), 𝑋₂(𝑧),…,𝑋_𝑛(𝑧)) where each map 𝑋_𝑖∶[0,1]→ℝ, 1≤𝑖≤𝑛, is smooth over the interval (0,1). We consider the following vector differential equation satisfied by 𝐗(𝑧): 𝑑𝐗(𝑧)((𝑧[].𝑑𝑧=𝑓𝑧,𝐗)),𝑧∈0,1 (1) Throughout, the map 𝑓∶[0,1]×ℝ^𝑛→ℝ^𝑛 is assumed to be smooth and satisfying the Lipschitz condition on a closed rectangle 𝑅∶=[0,1]×𝑈⊂[0,1]×ℝ^𝑛 with a Lipschitz constant 𝐾>0 s.t. 𝑓(𝑧,𝐗(𝑧))−𝑓(𝑧,𝐘(𝑧))≤𝐾‖𝐗(𝑧)−𝐘(𝑧)‖, (2) for all  𝐗,𝐘∈𝑈 and 𝑧∈[0,1]; furthermore, the components of 𝐗(𝑧) are assumed to satisfy 𝑚 initial conditions at first boundary given by 𝑋₁,𝑋₂,…,𝑋_𝑚(0)=𝟎∈ℝ^𝑚, (3) and (𝑛−𝑚)-conditions at the other boundary which are given by 𝑋_𝜎(1),…,𝑋_{𝜎(𝑛−𝑚)}(1)=𝜇₁,…,𝜇_𝑛−𝑚∈ℝ^𝑛−𝑚, (4) 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 [1–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 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 ‖𝐗(𝑎+ℎ)−𝐘(𝑎+ℎ)‖≤𝑒^𝐾|ℎ|‖𝐗(𝑎)−𝐘(𝑎)‖. (5)

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 𝐹^′(𝐱)−𝐹^′(𝐲)≤𝐿‖𝐱−𝐲‖,∀𝐱,𝐲,∈𝑈₀. (6) For some 𝐱₀∈𝑈₀ assume that 𝐹^′(𝐱₀)⁻¹ is defined on all of 𝐵 and that ℎ∶=𝐿𝛽𝜂≤1/2, where ‖𝐹^′(𝐱₀)⁻¹‖≤𝛽 and ‖𝐹^′(𝐱₀)⁻¹𝐹(𝐱₀)‖≤𝜂. Define 𝑡^∗∶=(1−√1−2ℎ)/𝛽𝐿, 𝑡^∗∗∶=(1+√1−2ℎ)/𝛽𝐿, and 𝑆∶={𝐱∈𝑈₀∣‖𝐱−𝐱₀‖≤𝑡^∗}. Then the Newton iterates 𝐱_𝑘+1=𝐱_𝑘−𝐹^′𝐱_𝑘⁻¹𝐹𝐱_𝑘,𝑘=0,1,2… (7) are well-defined and lie in 𝑆 and converge to a solution 𝐱^∗ of 𝐹(𝐱^∗)=0 which is unique in 𝑈₀∩{𝐱∣‖𝐱−𝐱₀‖<𝑡^∗∗}. Moreover if ℎ<1/2, 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 𝑑𝐗()((;𝐗()(0)()(1).𝑑𝑧𝑧,𝐜=𝑓𝑧,𝐗𝑧,𝐜))0,𝐜=𝐗;𝐗1,𝐜=𝐗 (8)

Assume that 𝐗(𝑧,𝐜) and 𝑓 are defined on the closed rectangle 𝑅∶=[0,1]×𝑈⊂[0,1]×ℝ^𝑛−𝑚 and are continuously differentiable on 𝑅−𝜕𝑅 such that 𝜕𝑓/𝜕𝐗 is nonsingular there. Observe that at point (𝑧,𝑐_𝑖),  𝑖=1,2,…𝑛−𝑚, 𝑑𝑑𝑧𝑑𝐗𝑑𝑐_𝑖=𝑑𝑑𝑐_𝑖𝑑𝐗=𝑑𝑑𝑧𝑑𝑐_𝑖𝑓𝑧,𝐗𝑧,𝑐_𝑖=𝜕𝑓𝜕𝐗𝑑𝐗𝑑𝑐_𝑖. (9)

If 𝐗(0,𝐜)=(0,0,…,0_𝑚-times,𝑐₁,𝑐₂,…,𝑐_𝑛−𝑚) then we see that 𝑑𝐗/𝑑𝑐_𝑖=(0,0,…,0,0,0,…0,1_}(𝑚+𝑖),0,…,0)=𝐞_𝑚+𝑖, 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) 𝑑𝑑𝑧𝑑𝐗=𝑑𝑐𝜕𝑓𝜕𝐗𝑑𝐗,𝑑𝑐𝑑𝐗(0)𝑑𝑐=𝐞_𝑘, (10) one for each 𝑘=𝑚+1,…,𝑛. 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 𝐹𝐜^𝑘=𝐗1,𝐜^𝑘−𝐗(1),𝑘=0,1,2,…, (11)
on the closed rectangle 𝑈⊂ℝ^𝑛−𝑚 containing point 𝐜 in its interior where 𝐗(𝑧,𝐜^𝑘)∈[0,1]×𝑈 s.t. 𝐗(0)=0;  𝐗(1)  =  (𝑋_𝜎(1),𝑋_𝜎(2),…,𝑋_{𝜎(𝑛−𝑚)})(1), 𝐗(0,𝐜^𝑘)=𝐜^𝑘; 𝐗(1,𝐜^𝑘)  =  (𝑋_𝜎(1),𝑋_𝜎(2),…,𝑋_{𝜎(𝑛−𝑚)})(1,𝐜^𝑘). The sequence {𝐜^𝑘} is defined recursively by the following: 𝐜⁰=𝐜,𝐜^𝑘+1=𝐜^𝑘−𝐹^′𝐜^𝑘−1𝐹𝐜^𝑘,𝑘=1,2,…. (12)

Proof. Since the Jacobian matrix 𝐗^′(𝑧,𝐜)∶=𝜕𝐗/𝜕𝐜(𝑧,𝐜) satisfies the linear matrix differential equation 𝑑𝐗𝑑𝑧^′(𝑧,𝐜)=𝑓^′𝑧,𝐗(𝑧,𝐜)𝐗^′(),𝐗𝑧,𝐜_′(0,𝐜)=𝛿_{𝑖𝑗(𝑛−𝑚)×(𝑛−𝑚)}, (13)
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 𝐗^′(1,𝐜)=𝐹^′(𝐜) 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 𝜔∶=_𝑖,𝑗sup_𝑡∈(0,1)𝐽_𝑖𝑗(𝑡);𝐽_𝑖𝑗(𝑡)∶=𝜕𝑓((𝜕𝐗𝑡,𝐗𝑡,𝐜))_(𝑖,𝑗). (14)
Note that 𝜔>0 as (𝜕𝑓/𝜕𝐗)(𝑡,𝐗(𝑡,𝐜)) is nonsingular and hence not zero identically. Consequently, with an application of Theorem 3, we arrive at the following: 𝐹^′𝐜^𝑘+1−𝐹^′𝐜^𝑘=𝐗^′1,𝐜^𝑘+1−𝐗^′1,𝐜^𝑘≤𝑒^𝜔𝐗^′0,𝐜^𝑘+1−𝐗^′0,𝐜^𝑘=0, (15)
showing that {𝐹^′(𝐜^𝑘)},𝑘=0,1,…, is a constant sequence. It follows that 𝐹^′ satisfies the Lipschitz condition on 𝑈₀∋𝐜 with any real number 𝜀>0 as a Lipschitz constant. We choose 𝜀>0 small enough such that for ‖𝐹^′(𝐜)⁻¹‖≤𝛽≠0 and ‖𝐹^′(𝐜)⁻¹𝐹(𝐜)‖≤𝜂≠0, the following strict inequality holds: 1ℎ∶=𝜀𝛽𝜂<2. (16) This completes the proof.

Corollary 7. The sequence of Newton iterates {𝐜^𝑘}, 𝑘=0,1,… defined by (12) of the Theorem 6, is well-defined and lies in the set {𝐱∣‖𝐱−𝐜‖<𝑡^∗}⊂𝑈₀ and converges to unique solution 𝑐^∗ of the equation 𝐹(𝐱)=0 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 0=𝑧₀<𝑧₁<⋯<𝑧_𝑁−1<𝑧_𝑁=1,  𝑧_𝑖=𝑖/𝑁,  1≤𝑖≤𝑁, be a partition of [0,1]. 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 𝐗(𝑧,𝐜^𝑘+1) and 𝐗(𝑧,𝐜^𝑘) of (8): 𝐗𝑧_𝑖,𝐜^𝑘+1−𝐗𝑧_𝑖,𝐜^𝑘≤𝑒^{𝐾(𝑧𝑖−𝑧𝑖−1)}𝐗𝑧_𝑖−1,𝐜^𝑘+1−𝐗𝑧_𝑖−1,𝐜^𝑘⋯≤𝑒^{𝐾(𝑧𝑖−𝑧0)}𝐗𝑧₀,𝐜^𝑘+1−𝐗𝑧₀,𝐜^𝑘, (17)
and correspondingly, the following is satisfied: 𝐗𝑧_𝑖,𝐜^𝑘+1−𝐗𝑧_𝑖,𝐜^𝑘≤𝑒^{𝐾(𝑧𝑖−𝑧𝑖−1)}𝐗𝑧_𝑖−1,𝐜^𝑘+1−𝐗𝑧_𝑖−1,𝐜^𝑘⋯≤𝑒^{𝐾(𝑧𝑖−𝑧0)}𝐗𝑧₀,𝐜^𝑘+1−𝐗𝑧₀,𝐜^𝑘. (18)
From (18) we obtain the desired assertion for the function 𝐹(𝐜^𝑘),𝑘=0,1…, as follows: 𝐹𝐜^𝑘+1−𝐹𝐜^𝑘=𝐗1,𝐜^𝑘+1−𝐗1,𝐜^𝑘≤𝑒^𝐾𝐗0,𝐜^𝑘+1−𝐗0,𝐜^𝑘=𝑒^𝐾𝐜^𝑘+1−𝐜^𝑘=𝑒^𝐾𝐹^′(𝐜)⁻¹𝐹𝐜^𝑘, (19)
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 𝐹(𝑥)=0 (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 (𝑘+1)th iteration involves simultaneous numerical integration of (8) from 𝑧=0 to 𝑧=1 to obtain 𝐹(𝐜^𝑘) and yet another numerical integration of the IVPs defined by the linear system in (10) for each 𝑐_𝑖 from 𝑧=0 to 𝑧=1 to obtain 𝐗^′(1,𝐜^𝑘). Since the underlying numerical integrations are to be carried out on interval [0,1], we define another equivalent vector differential equation as follows: 𝑑𝐗𝑑𝑧𝑑𝐗𝑑𝑐^𝑘𝑧,𝑐^𝑘=𝑓𝑧,𝐗𝑧,𝑐^𝑘𝜕𝑓𝜕𝐗𝑑𝐗𝑑𝑐^𝑘𝑧,𝑐^𝑘;𝐗𝑑𝐗𝑑𝑐^𝑘0,𝑐^𝑘=𝐞_𝑚+𝑖,foreach𝑖=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 ^𝐗_{𝑑𝐗/𝑑𝑐𝑘}(1,𝑐^𝑘) for each 𝑐^𝑘=𝑐^𝑘₁,…,𝑐^𝑘_𝑛−𝑚 s.t. 𝐜^𝑘=(𝑐^𝑘₁,…,𝑐^𝑘_𝑛−𝑚),  𝑘=0,1,… at each Newton iteration from which it is easy to extract 𝐗(1,𝐜^𝑘) and the Jacobian matrix 𝐗^′(1,𝐜^𝑘) given by 𝐗^′1,𝐜^𝑘𝑑𝐗∶=𝑑𝑐^𝑘₁,𝑑𝐗𝑑𝑐^𝑘₂𝑑𝐗,…,𝑑𝑐^𝑘_𝑛−𝑚1,𝐜^𝑘. (21)

Now each Newton iteration is performed to obtain 𝐜^𝑘+1 by solving the following system of (𝑛−𝑚) linear algebraic equations in 𝑐₁^𝑘+1,…,𝑐^𝑘+1_𝑛−𝑚 given by 𝐗^′1,𝐜^𝑘𝐜^𝑘+1=𝐗^′1,𝐜^𝑘𝐜^𝑘−𝐗1,𝐜^𝑘−𝐗(1). (22)

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 ^̇𝑋1_̇𝑋2=^𝑋2_{2𝑋31−6𝑋1−2𝑧3},  𝑧∈[1,2], with the boundary conditions 𝑋₁(1)=2,  𝑋₁(2)=5/2. The exact solution on [1,2] is 𝑋₁(𝑧)=𝑧+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⁻⁵. Table 1 gives the description of the convergence analysis of the solution for several values of the trial initial guess for 𝐜⁰=(2,𝑐⁰₁)^′. It is clear from the table that the initial guess 𝑐⁰₁=0 is appropriate for which the numerical scheme converges at the very first iteration to zero within an error of the order of 10⁻¹². If we take the initial guess 𝑐⁰₁ 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⁻⁶ which corresponds to the value 𝑐⁰₁=10.

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 𝑐⁰₁≥1. We have not encountered this situation except for possibly 𝑐⁰₁>10 when the initial guess does not fall into the the disc of radius 𝑡^∗∗ (see Theorem 5). Also for 𝑐⁰₁=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 𝑅∶={(𝑋₁,𝑋₂)∈ℝ²∣‖𝑋₁−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 𝑋₁=2, 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.

(a)

(b)

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 𝑧=0 and 𝑧=𝑑 (see 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.

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]: 𝜌₀𝜕𝐮𝜕𝑡+𝜌₀𝐮⋅∇𝐮=−∇𝑝+𝜌₀𝜈∇²𝐮+𝜌𝐠,𝜕𝑇𝑘𝜕𝑡+𝐮⋅∇𝑇=^𝑇𝜌₀𝐶_𝑉∇²𝑇,∇⋅𝐮=0,𝜌=𝜌₀1−𝛼𝑇−𝑇_𝑎, (23) where 𝜌₀ is the fluid density at some reference fluid temperature 𝑇_𝑎,  𝜈 is the fluid viscosity, 𝐠=(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 𝐮_𝑒=(0,0,0),𝑇_𝑒=𝑇₁−𝑇₁−𝑇₂𝑧𝑑,𝑑𝑝_𝑒𝑑𝑧=−𝜌₀𝑔1−𝛼𝑇_𝑒−𝑇_𝑎, (24) 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: 𝑘𝐮=^𝑇𝑑(𝑢,𝑣,𝑤),𝑇=𝑇_𝑒+𝑇₁−𝑇₂𝜃,𝑝=𝑝_𝑒+𝜌₀𝑘²_𝑇𝑑⁴𝑃, (25) 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 𝜕𝜕𝑡−Pr∇²+𝐮⋅∇𝜕𝑢−𝜕𝑦𝜕𝑣+𝜕()𝜕𝑥𝑤,𝑢𝜕+𝜕()𝑦,𝑧𝑤,𝑣𝜕()𝜕𝑧,𝑥=0,𝜕𝑡−Pr∇²+𝐮⋅∇𝜕𝑣−𝜕𝑧𝜕𝑤+𝜕()𝜕𝑦𝑢,𝑣𝜕()+𝜕()𝑧,𝑥𝑢,𝑤𝜕𝑥,𝑦=−PrRa𝜕𝜃,𝜕𝜕𝑦𝜕𝑡−Pr∇²+𝐮⋅∇𝜕𝑤−𝜕𝑥𝜕𝑢+𝜕()𝜕𝑧𝑣,𝑤𝜕+𝜕()𝑥,𝑦𝑣,𝑢𝜕𝑦,𝑧=PrRa𝜕𝜃,𝜕𝜕𝑥𝜕𝑡−Pr∇²+𝐮⋅∇𝜃=−𝐷𝑇𝑤,𝜕𝑢+𝜕𝑥𝜕𝑣+𝜕𝑦𝜕𝑤𝜕𝑧=0, (26) where 0≤𝑧≤1,𝑡>0; the dimensionless numbers Pr∶=𝜈/𝑘_𝑇 and Ra∶=(𝑇₁−𝑇₂)𝛼𝑑³𝑔/𝑘_𝑇𝜈 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:𝐮|_𝑧=0,1=0,D𝑤|_𝑧=0,1=0,𝜃|_𝑧=0,1=0,Free-Boundaries:𝐮|_𝑧=0,1D=0,₂𝑤|_𝑧=0,1=0,𝜃|_𝑧=0,1=0. (27)

To solve the nonlinear system of (26), we consider the solution in the following form: 𝑢𝑣𝑤𝜃=𝑥,𝑦,𝑧,𝑡_{∞∞ℓ=0𝑚=0}𝑢_ℓ,𝑚𝑣_ℓ,𝑚𝑤_ℓ,𝑚𝜃_ℓ,𝑚(𝑧)exp𝜔_ℓ𝑡+𝜄ℓ𝑘₁𝑥+𝜄𝑚𝑘₂𝑦, (28) 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 𝑘₁ and 𝑘₂: 𝑘²₁ℓ²+𝑘²₂𝑚²=(ℓ+𝑚)²𝑘²₁. (29) Consequently, for 𝑚=0 the choice for 𝑘₂ is arbitrary while for ℓ=0, 𝑘₁=𝑘₂. Therefore, for a choice 𝑘₁=𝑘₂, 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 𝑘₁→((ℓ+𝑚)²−ℓ²/𝑚)𝑘₁. The solution corresponding to 𝑚=0 is called the rolls solution, and for ℓ=𝑚 the solution is called the square pattern.

3.1. Rolls Solution

The underlying hydrodynamic system consists of 6ℓ ODEs for each normal mode ℓ=1,2,… and 𝑚=0. The linearized system corresponds to the lowest normal mode ℓ=1. Higher modes include the nonlinear effects. For a typical value ℓ=3 the dynamical system is of order eighteen and is described by the following set of ODEs: 𝑋̇_ℓ1𝑋̇_ℓ2𝑋̇_ℓ3𝑋̇_ℓ4𝑋̇_ℓ5𝑋̇_ℓ6=𝑋_ℓ2𝑋_ℓ4𝑋_ℓ5𝑋_ℓ6−𝑋_ℓ1+𝑘²_ℓ+𝜔_ℓ𝑋_ℓ3+∑^ℓ−1_𝑝=1ℓ1−𝑝𝑋_𝑝2𝑋_{(ℓ−𝑝)3}+𝑋_𝑝1𝑋_{(ℓ−𝑝)5}−𝑘⁴_ℓ+𝜔_ℓ𝑘²_ℓ𝑋Pr^ℓ1+2𝑘²_ℓ+𝜔_ℓ𝑋Pr^ℓ4+𝑘²_ℓRa𝑋_ℓ3+ℓPr^ℓ−1_𝑝∑₌₁𝑁_ℓ,𝑝, (30)where 𝑋_ℓ1,𝑋_ℓ2,𝑋_ℓ3,𝑋_ℓ4,𝑋_ℓ5,𝑋_ℓ6=𝑤_ℓ,0,𝐷𝑤_ℓ,0,𝜃_ℓ,0,𝐷²𝑤_ℓ,0,𝐷𝜃_ℓ,0,𝐷³𝑤_ℓ,0, (31)𝑁_ℓ,𝑝∶={−𝑋_𝑝2𝑋_{(ℓ−𝑝)4}/𝑝+𝑋_𝑝1𝑋_{(ℓ−𝑝)6}/ℓ−𝑝}, ℓ=1,2,3, 𝑋̇_ℓ𝑗=𝑑/𝑑𝑧𝑋_ℓ𝑗, for all 1≤𝑗≤6; 𝑘²_ℓ=ℓ²𝑘²₁; for ℓ=1 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 𝑋_ℓ1,𝑋_ℓ2,𝑋_ℓ3(0)=()=𝑋0,0,0_ℓ1,𝑋_ℓ2,𝑋_ℓ3(1),∀ℓ=1,2,3. (32)

The vector differential equation (30) involves unknown constants 𝜔_ℓ,𝑘₁, and Ra for a given value of Pr. As a check we take (𝜔₁,𝜔₂,𝜔₃,Pr,𝑘₁,Ra)=(−1000,0,0,10,1.603,15874).

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 𝐜=0∈𝑅={𝐱∈ℝ⁹∣‖𝐱−𝐜‖<2}. The missing boundary conditions for 𝑧=0 are found to be given by 𝑋₁₄𝑋₁₅𝑋₁₆𝑋₂₄𝑋₂₅𝑋₂₆𝑋₃₄𝑋₃₅𝑋₃₆0,𝐜^∗=𝐜^∗=+0.023007339918880+0.023171813739103−1.000000000000000−0.000017193958663−0.000000054567985+0.000192804606564+0.000000053314000+0.000000000027238−0.000000563466641. (33) 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 (𝜔1,𝜔2,𝜔3,Pr,𝑘1,Ra)=(−1000,0,0,10,1.603,15874).

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 ℓ=1 and 𝑚=1. Here: the underlying nonlinear dynamical system is given by the following twelve ODEs: 𝑋̇_1𝑚1𝑋̇_1𝑚2𝑋̇_1𝑚3𝑋̇_1𝑚4𝑋̇_1𝑚5𝑋̇_1𝑚6=𝑋_1𝑚2𝑋_1𝑚4𝑋_1𝑚5𝑋_1𝑚6𝐿¹_𝑚+2𝑚𝑋_1𝑚1𝑋_1𝑚5𝐿²_𝑚+𝑁_𝑚 (34) for 𝑚=0,1 where 𝑘_𝑚∶=(𝑚+1)𝑘₁; 𝐿¹_𝑚∶=−𝑋_1𝑚1+(𝑘²_𝑚+𝜔_𝑚)𝑋_1𝑚3, 𝐿²_𝑚∶=−𝑘⁴_𝑚+𝜔_𝑚𝑘²_𝑚𝑋Pr^1𝑚1+2𝑘²_𝑚+𝜔_𝑚𝑋Pr^1𝑚4+𝑘²_𝑚Ra𝑋_1𝑚3,𝑁_𝑚∶=2𝑚𝑋Pr¹⁰¹𝑋₁₀₆+𝑋₁₀₂𝑋₁₀₄−𝑘²_𝑚𝑋₁₀₁𝑋₁₀₂, (35) and the twelve boundary conditions are given by 𝑋_1𝑚1,𝑋_1𝑚2𝑋_1𝑚3(0)=()=𝑋0,0,0_1𝑚1,𝑋_1𝑚2𝑋_1𝑚3(1), (36) for each 𝑚=0,1. The fixed parametric values are taken as the following 𝜔₀,𝜔₁,Pr,𝑘₁,Ra=(85.80,0,10,1.603,15874). (37)

In this case, the missing boundary conditions for 𝑧=0 using an initial guess of 𝐜=0 in the rectangle 𝑅 defined by 𝑅∶={𝐱∈ℝ⁶∣‖𝐱−𝐜‖<2200} are found to be the following: 𝑋₁₀₄𝑋₁₀₅𝑋₁₀₆𝑋₁₁₄𝑋₁₁₅𝑋₁₁₆0,𝐜^∗=𝐜^∗=+79.20602598022311+0.045512918963491−500.0000000000000−129.3119890195200−0.119928096846000+1357.197011463896. (38)

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.

(a)

(b)

Figure 4 

Square pattern in the Rayleigh-Bénard convection for (𝜔0,𝜔1,Pr,𝑘1,Ra)=(85.80,0,10,1.603,15874).

(a)

(b)

Figure 5 

Square pattern (side view in 𝑥𝑧-plane) in the Rayleigh-Bénard convection for (𝜔0,𝜔1,Pr,𝑘1,Ra)=(85.80,0,10,1.603,15874).

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.