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

Authors: Jitender Singh
Journal: ISRN Mathematical Physics
 (2013)

---

## Abstract

The simple shooting method is revisited in order to solve nonlinear two-point BVP numerically. The BVP of the typeX'(z)=f(z,X(z)),X(z)∈ℝn,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.


---

## Body

## 1. Introduction

LetX:[0,1]→ℝn,  n>1, be a vector valued function defined by X(z)=(X1(z), X2(z),…,Xn(z)) where each map Xi:[0,1]→ℝ, 1≤i≤n, is smooth over the interval (0,1). We consider the following vector differential equation satisfied by X(z):
(1)ddzX(z)=f(z,X(z)),z∈[0,1].
Throughout, the map f:[0,1]×ℝn→ℝn is assumed to be smooth and satisfying the Lipschitz condition on a closed rectangle R:=[0,1]×U-⊂[0,1]×ℝn with a Lipschitz constant K>0 s.t.                      (2)∥f(z,X(z))-f(z,Y(z))∥≤K∥X(z)-Y(z)∥,
for all  X,Y∈U- and z∈[0,1]; furthermore, the components of X(z) are assumed to satisfy m initial conditions at first boundary given by
(3)(X1,X2,…,Xm)(0)=0∈ℝm,
and (n-m)-conditions at the other boundary which are given by
(4)(Xσ(1),…,Xσ(n-m))(1)=(μ1,…,μn-m)∈ℝn-m,
where σ is a permutation of the symmetric group Sn.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(z) 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 knowna priori. The degree of difficulty increases as the dimension n 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 advancedparallel 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).
Letf(z,X(z)) be continuous and satisfying Lipschitz condition as in (2) on interval |z-a|≤T for all X∈U-. Then for any constant vector c, the differential equation (1) has a solution on |z-a|≤T with initial condition X(a)=c.Theorem 2 (uniqueness).
Iff(z,X(z)) satisfies Lipschitz condition on U-, then there is at most one solution of (1) with initial condition X(a)=c.Theorem 3 (continuity).
IfX(z) and Y(z) are any two solutions of (1) where f(z,X(z)) is continuous and satisfies Lipschitz condition of (2), then
(5)∥X(a+h)-Y(a+h)∥≤eK|h|∥X(a)-Y(a)∥.Corollary 4.
LetX(z,c) be a solution of (1) with X(z,c)=c and hypothesis of Theorem 3 being satisfied. Let X(z,c) be defined on ∥c-c0∥≤K, |z-a|≤T. Then X(z,c) is a continuous function of z and c. Further, if c→c0, then X(z,c)→X(z,c0) uniformly for |z-a|≤T.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).
LetA and B be Banach spaces and F:U⊂A→B such that on an open convex set U0⊂U, the map F is Fréchèt differentiable and satisfies
(6)∥F′(x)-F′(y)∥≤L∥x-y∥,∀x,y,∈U0.
For some x0∈U0 assume that F′(x0)-1 is defined on all of B and that h:=Lβη≤1/2, where ∥F′(x0)-1∥≤β and ∥F′(x0)-1F(x0)∥≤η. Define t*:=(1-1-2h)/βL, t**:=(1+1-2h)/βL, and S:={x∈U0∣∥x-x0∥≤t*}. Then the Newton iterates
(7)xk+1=xk-F′(xk)-1F(xk),k=0,1,2…
are well-defined and lie in S and converge to a solution x* of F(x*)=0 which is unique in U0∩{x∣∥x-x0∥<t**}. Moreover if h<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 Section2 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 vectorc∈ℝn-m, n>m, independent of the variable z s.t. c:=(c1,c2,…,cn-m) and the smooth function X(z,c) satisfies (1)–(4); that is
(8)ddzX(z,c)=f(z,X(z,c));X(0,c)=X(0);X(1,c)=X(1).Assume thatX(z,c) and f are defined on the closed rectangle R:=[0,1]×U-⊂[0,1]×ℝn-m and are continuously differentiable on R-∂R such that ∂f/∂X is nonsingular there. Observe that at point (z,ci),  i=1,2,…n-m,
(9)ddz(dXdci)=ddcidXdz=ddcif(z,X(z,ci))=∂f∂XdXdci.IfX(0,c)=(0,0,…,0︸m-times,c1,c2,…,cn-m) then we see that dX/dci=(0,0,…,0,0,0,…0,1︸(m+i),0,…,0)=em+i, where em+i is an element of the standard basis of the vector space ℝn(ℝ). It follows that for each component ci of c, the derivative dX/dci satisfies (n-m) initial value problems (IVPs)
(10)ddz(dXdc)=∂f∂X(dXdc),(dXdc)(0)=ek,
one for each k=m+1,…,n. Under these considerations, we have the following main theorem.Theorem 6.
The hypothesis of Theorem5 is satisfied by the sequence {F(ck)} of functions defined by
(11)F(ck)=X~(1,ck)-X~(1),k=0,1,2,…,
on the closed rectangleU-⊂ℝn-m containing point c in its interior where X~(z,ck)∈[0,1]×U- s.t. X~(0)=0;  X~(1)  =  (Xσ(1),Xσ(2),…,Xσ(n-m))(1),  X~(0,ck)=ck; X~(1,ck)  =  (Xσ(1),Xσ(2),…,Xσ(n-m))(1,ck). The sequence {ck} is defined recursively by the following:
(12)c0=c,ck+1=ck-(F′(ck))-1F(ck),k=1,2,….Proof.
Since the Jacobian matrixX~′(z,c):=∂X~/∂c(z,c) satisfies the linear matrix differential equation
(13)ddzX~′(z,c)=f~′(z,X~(z,c))X~′(z,c),X~′(0,c)=(δij)(n-m)×(n-m),
wheref~ is the vector with (n-m) components of f depending upon the permutation σ, it follows from the existence and uniqueness theorem for the system of ODEs that X~′(1,c)=F′(c) is defined and invertible in an open rectangle U0⊂U- containing the point c. Also, since X~′(z,c) is a solution of a linear vector differential equation in z,f′(z,X~(z,c)) satisfies the Lipschitz condition with a Lipschitz constant
(14)ω:=∑i,j‍supt∈(0,1)∥Jij(t)∥;Jij(t):=[∂f∂X(t,X(t,c))](i,j).
Note thatω>0 as (∂f/∂X)(t,X(t,c)) is nonsingular and hence not zero identically. Consequently, with an application of Theorem 3, we arrive at the following:
(15)∥F′(ck+1)-F′(ck)∥=∥X~′(1,ck+1)-X~′(1,ck)∥≤eω∥X~′(0,ck+1)-X~′(0,ck)∥=0,
showing that{F′(ck)},k=0,1,…, is a constant sequence. It follows that F′ satisfies the Lipschitz condition on U0∋c with any real number ε>0 as a Lipschitz constant. We choose ε>0 small enough such that for ∥F′(c)-1∥≤β≠0 and ∥F′(c)-1F(c)∥≤η≠0, the following strict inequality holds:
(16)h:=εβη<12.
This completes the proof.Corollary 7.
The sequence of Newton iterates{ck}, k=0,1,… defined by (12) of the Theorem 6, is well-defined and lies in the set {x∣∥x-c∥<t*}⊂U0 and converges to unique solution c* of the equation F(x)=0 on U0∩{x∣∥x-c∥<t**}. Furthermore, X~(z,c)→X~(z,c*) as c→c*.Proof.
The convergence of the Newton iterates follows from the Newton-Kantorovich Theorems5 and 6 with L=ε and convergence of the solution X~(z,c) to a solution X~(z,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{F(ck)} as aforementioned satisfies the Lipschitz condition with a Lipschitz constant eK and hence is dominated by the sequence {ck}. To establish this, let 0=z0<z1<⋯<zN-1<zN=1,  zi=i/N,  1≤i≤N, be a partition of [0,1]. Since f is given to satisfy the hypothesis of Theorem 3 with Lipschitz constant K, it follows that for a fixed k and each a=zi-1, the following holds for the two solutions X(z,ck+1) and X(z,ck) of (8):
(17)∥X(zi,ck+1)-X(zi,ck)∥≤eK(zi-zi-1)∥X(zi-1,ck+1)-X(zi-1,ck)∥⋯≤eK(zi-z0)∥X(z0,ck+1)-X(z0,ck)∥,
and correspondingly, the following is satisfied:(18)∥X~(zi,ck+1)-X~(zi,ck)∥≤eK(zi-zi-1)∥X~(zi-1,ck+1)-X~(zi-1,ck)∥⋯≤eK(zi-z0)∥X~(z0,ck+1)-X~(z0,ck)∥.
From (18) we obtain the desired assertion for the function F(ck),k=0,1…, as follows:
(19)∥F(ck+1)-F(ck)∥=∥X~(1,ck+1)-X~(1,ck)∥≤eK∥X~(0,ck+1)-X~(0,ck)∥=eK∥ck+1-ck∥=eK∥(F′(c))-1F(ck)∥,
where we have made use of the fact that the sequence{F′(ck)} is constant. It also follows from the estimate (19) that F(ck)→F(c*)→0 for ck→c*.Remark 9.
It is also known that ifF′(x) satisfies Lipschitz condition then the underlying Newton iterates have at least quadratic convergence to a root of F(x)=0 (see Rall [15]).


### 2.1. Numerical Scheme

The associated numerical scheme involves computation of numerical solutionX(z,c*) of the BVP defined by (1)–(4) which at (k+1)th iteration involves simultaneous numerical integration of (8) from z=0 to z=1 to obtain F(ck) and yet another numerical integration of the IVPs defined by the linear system in (10) for each ci from z=0 to z=1 to obtain X′(1,ck). Since the underlying numerical integrations are to be carried out on interval [0,1], we define another equivalent vector differential equation as follows:
(20)ddz(XdXdck)(z,ck)=(f(z,X(z,ck))∂f∂XdXdck(z,ck));(XdXdck)(0,ck)=em+i,foreachi=1,…,(n-m).These are(n-m) IVPs, each of order 2n and can be integrated numerically using the fourth order Runge-Kutta method over the interval [0,1], to obtain (XdX/dck)(1,ck) for each ck=c1k,…,cn-mk s.t. ck=(c1k,…,cn-mk),  k=0,1,… at each Newton iteration from which it is easy to extract X~(1,ck) and the Jacobian matrix X~′(1,ck) given by
(21)X~′(1,ck):=(dX~dc1k,dX~dc2k,…,dX~dcn-mk)(1,ck).Now each Newton iteration is performed to obtainck+1 by solving the following system of (n-m) linear algebraic equations in c1k+1,…,cn-mk+1 given by
(22)X~′(1,ck)ck+1=X~′(1,ck)ck-{X~(1,ck)-X~(1)}.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(X˙1X˙2)=(X22X13-6X1-2z3),  z∈[1,2], with the boundary conditions X1(1)=2,  X1(2)=5/2. The exact solution on [1,2] is X1(z)=z+1/z. We choose step size h=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 c0=(2,c10)′. It is clear from the table that the initial guess c10=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 c10 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 c10=10.Table 1
Dependence of convergence of the numerical scheme on the initial guess.
c
1
0
k
c
1
k
max
z
∈
[
0,1
]
{
∥
X
~
(
z
,
c
k
)
-
X
~
(
z
)
∥
}
−1.0
18
1.4551
e
-
007
1.7920
e
-
006
0.0
1
0.0000
e
-
000
5.8100
e
-
012
0.1
5
1.2995
e
-
011
1.6582
e
-
010
0.5
12
1.5494
e
-
007
1.9081
e
-
006
1.0
14
1.8527
e
-
008
2.2818
e
-
007
5.0
33
1.8362
e
-
009
2.2619
e
-
008
10.0
30
7.8688
e
-
007
9.6907
e
-
006This example has been considered by Ha [7] where no suitable solution was found via his numerical experiment with the same BVP for c10≥1. We have not encountered this situation except for possibly c10>10 when the initial guess does not fall into the the disc of radius t** (see Theorem 5). Also for c10=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.3561e-008. Note that the system of ODEs considered here is locally Lipschitz with Lipschitz constant K≥40 on the closed rectangular strip of width 1 defined by R:={(X1,X2)∈ℝ2∣∥X1-2∥≤0.5}; therefore, by Theorem 6, it follows that the present numerical scheme will always converge for each initial guess taken from R on the line X1=2, but with large c10, the speed of convergence will slow down as indicated in Table 1.The numerical solutionsX1 and X2 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.The numerical solution⋇ versus the actual solution of the BVP of Example 10.
(a)
(b)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.


## 2.1. Numerical Scheme

The associated numerical scheme involves computation of numerical solutionX(z,c*) of the BVP defined by (1)–(4) which at (k+1)th iteration involves simultaneous numerical integration of (8) from z=0 to z=1 to obtain F(ck) and yet another numerical integration of the IVPs defined by the linear system in (10) for each ci from z=0 to z=1 to obtain X′(1,ck). Since the underlying numerical integrations are to be carried out on interval [0,1], we define another equivalent vector differential equation as follows:
(20)ddz(XdXdck)(z,ck)=(f(z,X(z,ck))∂f∂XdXdck(z,ck));(XdXdck)(0,ck)=em+i,foreachi=1,…,(n-m).These are(n-m) IVPs, each of order 2n and can be integrated numerically using the fourth order Runge-Kutta method over the interval [0,1], to obtain (XdX/dck)(1,ck) for each ck=c1k,…,cn-mk s.t. ck=(c1k,…,cn-mk),  k=0,1,… at each Newton iteration from which it is easy to extract X~(1,ck) and the Jacobian matrix X~′(1,ck) given by
(21)X~′(1,ck):=(dX~dc1k,dX~dc2k,…,dX~dcn-mk)(1,ck).Now each Newton iteration is performed to obtainck+1 by solving the following system of (n-m) linear algebraic equations in c1k+1,…,cn-mk+1 given by
(22)X~′(1,ck)ck+1=X~′(1,ck)ck-{X~(1,ck)-X~(1)}.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(X˙1X˙2)=(X22X13-6X1-2z3),  z∈[1,2], with the boundary conditions X1(1)=2,  X1(2)=5/2. The exact solution on [1,2] is X1(z)=z+1/z. We choose step size h=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 c0=(2,c10)′. It is clear from the table that the initial guess c10=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 c10 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 c10=10.Table 1
Dependence of convergence of the numerical scheme on the initial guess.
c
1
0
k
c
1
k
max
z
∈
[
0,1
]
{
∥
X
~
(
z
,
c
k
)
-
X
~
(
z
)
∥
}
−1.0
18
1.4551
e
-
007
1.7920
e
-
006
0.0
1
0.0000
e
-
000
5.8100
e
-
012
0.1
5
1.2995
e
-
011
1.6582
e
-
010
0.5
12
1.5494
e
-
007
1.9081
e
-
006
1.0
14
1.8527
e
-
008
2.2818
e
-
007
5.0
33
1.8362
e
-
009
2.2619
e
-
008
10.0
30
7.8688
e
-
007
9.6907
e
-
006This example has been considered by Ha [7] where no suitable solution was found via his numerical experiment with the same BVP for c10≥1. We have not encountered this situation except for possibly c10>10 when the initial guess does not fall into the the disc of radius t** (see Theorem 5). Also for c10=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.3561e-008. Note that the system of ODEs considered here is locally Lipschitz with Lipschitz constant K≥40 on the closed rectangular strip of width 1 defined by R:={(X1,X2)∈ℝ2∣∥X1-2∥≤0.5}; therefore, by Theorem 6, it follows that the present numerical scheme will always converge for each initial guess taken from R on the line X1=2, but with large c10, the speed of convergence will slow down as indicated in Table 1.The numerical solutionsX1 and X2 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.The numerical solution⋇ versus the actual solution of the BVP of Example 10.
(a)
(b)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 Theorem6 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 z=0 and z=d (see Figure 2).Figure 2
The schematic diagram for Rayleigh-Bénard convection in an infinite horizontal fluid layer confined between the planesz=0 and z=d which are maintained at constant temperatures T1 and T2, respectively.The fluid layer is maintained at different constant temperaturesT1 and T2 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 t and any point P(x,y,z) in the fluid layer, u is the fluid velocity, T is the fluid temperature, p is the fluid pressure, and ρ is the fluid density, and then the flow is described by the following equations [17]:
(23)ρ0∂u∂t+ρ0u·∇u=-∇p+ρ0ν∇2u+ρg,∂T∂t+u·∇T=kTρ0CV∇2T,∇·u=0,ρ=ρ0{1-α(T-Ta)},
where ρ0 is the fluid density at some reference fluid temperature Ta,  ν is the fluid viscosity, g=(0,0,-g) is the acceleration due to gravity, kT is the diffusivity constant, CV is the specific heat of the fluid at constant volume, and α:=(∂ρ/∂T)|T=Ta.The system (23) admits a closed form steady state of the form
(24)ue=(0,0,0),Te=T1-(T1-T2)zd,dpedz=-ρ0g(1-α(Te-Ta)),
where the subscript e denotes the equilibrium state. We discuss the stability of the solution described by (24) by superimposing perturbations on it in the following form:
(25)u=kTd(u,v,w),T=Te+(T1-T2)θ,p=pe+ρ0kT2d4P,
where we have utilized the thickness d of the fluid layer as the length scale, kT/d as the velocity scale, and T1-T2 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
(26)(∂∂t-Pr∇2+u·∇)(∂u∂y-∂v∂x)+∂(w,u)∂(y,z)+∂(w,v)∂(z,x)=0,(∂∂t-Pr∇2+u·∇)(∂v∂z-∂w∂y)+∂(u,v)∂(z,x)+∂(u,w)∂(x,y)=-PrRa∂θ∂y,(∂∂t-Pr∇2+u·∇)(∂w∂x-∂u∂z)+∂(v,w)∂(x,y)+∂(v,u)∂(y,z)=PrRa∂θ∂x,(∂∂t-Pr∇2+u·∇)θ=-DTw,∂u∂x+∂v∂y+∂w∂z=0,
where 0≤z≤1,t>0; the dimensionless numbers Pr:=ν/kT and Ra:=(T1-T2)αd3g/kTν are called the Prandtl number and the Rayleigh number, respectively. The boundary conditions for the velocity field and the temperature field are given by
(27)Rigid-Boundaries:u|z=0,1=0,Dw|z=0,1=0,θ|z=0,1=0,Free-Boundaries:u|z=0,1=0,D2w|z=0,1=0,θ|z=0,1=0.To solve the nonlinear system of (26), we consider the solution in the following form:
(28)(uvwθ)(x,y,z,t)=∑ℓ=0∞‍∑m=0∞‍(uℓ,mvℓ,mwℓ,mθℓ,m)(z)exp{ωℓt+ιℓk1x+ιmk2y},
and we obtain the nonlinear system of ODEs as in (30) and (34) for a fixed value of (ℓ,m)=(3,0) for rolls solutions and (ℓ,m)=(1,1) for square patterns, respectively, for the case of rigid boundaries. For a fixed (ℓ,m)-mode such that ℓ+m=L, we assume the following relationship between the horizontal wave numbers k1 and k2:
(29)k12ℓ2+k22m2=(ℓ+m)2k12.
Consequently, for m=0 the choice for k2 is arbitrary while for ℓ=0, k1=k2. Therefore, for a choice k1=k2, the two problems for (ℓ,0) and (0,ℓ) modes are essentially the same. For ℓ≠0 and m≠0, the problem for (m,ℓ)-mode is equivalent to the corresponding problem for (ℓ,m)-mode under the change k1→((ℓ+m)2-ℓ2/m)k1. The solution corresponding to m=0 is called the rolls solution, and for ℓ=m the solution is called the square pattern.


### 3.1. Rolls Solution

The underlying hydrodynamic system consists of6ℓ ODEs for each normal mode ℓ=1,2,… and m=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: (30)(X˙ℓ1X˙ℓ2X˙ℓ3X˙ℓ4X˙ℓ5X˙ℓ6)=(Xℓ2Xℓ4Xℓ5Xℓ6-Xℓ1+(kℓ2+ωℓ)Xℓ3+∑p=1ℓ-1‍{(1-ℓp)Xp2X(ℓ-p)3+Xp1X(ℓ-p)5}-(kℓ4+ωℓkℓ2Pr)Xℓ1+(2kℓ2+ωℓPr)Xℓ4+kℓ2RaXℓ3+ℓPr∑p=1ℓ-1‍Nℓ,p),where
(31)(Xℓ1,Xℓ2,Xℓ3,Xℓ4,Xℓ5,Xℓ6)=(wℓ,0,Dwℓ,0,θℓ,0,D2wℓ,0,Dθℓ,0,D3wℓ,0),Nℓ,p:={-Xp2X(ℓ-p)4/p+Xp1X(ℓ-p)6/ℓ-p}, ℓ=1,2,3, X˙ℓj=d/dzXℓj, for all 1≤j≤6; kℓ2=ℓ2k12; 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
(32)(Xℓ1,Xℓ2,Xℓ3)(0)=(0,0,0)=(Xℓ1,Xℓ2,Xℓ3)(1),∀ℓ=1,2,3.The vector differential equation (30) involves unknown constants ωℓ,k1, and Ra for a given value of Pr. As a check we take (ω1,ω2,ω3,Pr,k1,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 c=0∈R={x∈ℝ9∣∥x-c∥<2}. The missing boundary conditions for z=0 are found to be given by
(33)(X14X15X16X24X25X26X34X35X36)(0,c*)=c*=(+0.023007339918880+0.023171813739103-1.000000000000000-0.000017193958663-0.000000054567985+0.000192804606564+0.000000053314000+0.000000000027238-0.000000563466641).
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, X11, X21, and X31 are the fluid velocity profiles corresponding to the linear mode, second order nonlinear mode, and the third order nonlinear mode, respectively. Similarly, X13, X23, and X33 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,k1,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 m=1. Here: the underlying nonlinear dynamical system is given by the following twelve ODEs:
(34)(X˙1m1X˙1m2X˙1m3X˙1m4X˙1m5X˙1m6)=(X1m2X1m4X1m5X1m6Lm1+2mX1m1X1m5Lm2+Nm)
for m=0,1 where km:=(m+1)k1; Lm1:=-X1m1+(km2+ωm)X1m3,
(35)Lm2:=-(km4+ωmkm2Pr)X1m1+(2km2+ωmPr)X1m4+km2RaX1m3,Nm:=2mPr(X101X106+X102X104-km2X101X102),
and the twelve boundary conditions are given by
(36)(X1m1,X1m2X1m3)(0)=(0,0,0)=(X1m1,X1m2X1m3)(1),
for each m=0,1. The fixed parametric values are taken as the following
(37)(ω0,ω1,Pr,k1,Ra)=(85.80,0,10,1.603,15874).In this case, the missing boundary conditions forz=0 using an initial guess of c=0 in the rectangle R defined by R:={x∈ℝ6∣∥x-c∥<2200} are found to be the following:
(38)(X104X105X106X114X115X116)(0,c*)=c*=(+79.20602598022311+0.045512918963491-500.0000000000000-129.3119890195200-0.119928096846000+1357.197011463896).Here the Newton iterates are also found to converge rapidly. Thenonlinear 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.Square pattern in the Rayleigh-Bénard convection for(ω0,ω1,Pr,k1,Ra)=(85.80,0,10,1.603,15874).
(a)
(b)Square pattern (side view inxz-plane) in the Rayleigh-Bénard convection for (ω0,ω1,Pr,k1,Ra)=(85.80,0,10,1.603,15874).
(a)
(b)


## 3.1. Rolls Solution

The underlying hydrodynamic system consists of6ℓ ODEs for each normal mode ℓ=1,2,… and m=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: (30)(X˙ℓ1X˙ℓ2X˙ℓ3X˙ℓ4X˙ℓ5X˙ℓ6)=(Xℓ2Xℓ4Xℓ5Xℓ6-Xℓ1+(kℓ2+ωℓ)Xℓ3+∑p=1ℓ-1‍{(1-ℓp)Xp2X(ℓ-p)3+Xp1X(ℓ-p)5}-(kℓ4+ωℓkℓ2Pr)Xℓ1+(2kℓ2+ωℓPr)Xℓ4+kℓ2RaXℓ3+ℓPr∑p=1ℓ-1‍Nℓ,p),where
(31)(Xℓ1,Xℓ2,Xℓ3,Xℓ4,Xℓ5,Xℓ6)=(wℓ,0,Dwℓ,0,θℓ,0,D2wℓ,0,Dθℓ,0,D3wℓ,0),Nℓ,p:={-Xp2X(ℓ-p)4/p+Xp1X(ℓ-p)6/ℓ-p}, ℓ=1,2,3, X˙ℓj=d/dzXℓj, for all 1≤j≤6; kℓ2=ℓ2k12; 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
(32)(Xℓ1,Xℓ2,Xℓ3)(0)=(0,0,0)=(Xℓ1,Xℓ2,Xℓ3)(1),∀ℓ=1,2,3.The vector differential equation (30) involves unknown constants ωℓ,k1, and Ra for a given value of Pr. As a check we take (ω1,ω2,ω3,Pr,k1,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 c=0∈R={x∈ℝ9∣∥x-c∥<2}. The missing boundary conditions for z=0 are found to be given by
(33)(X14X15X16X24X25X26X34X35X36)(0,c*)=c*=(+0.023007339918880+0.023171813739103-1.000000000000000-0.000017193958663-0.000000054567985+0.000192804606564+0.000000053314000+0.000000000027238-0.000000563466641).
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, X11, X21, and X31 are the fluid velocity profiles corresponding to the linear mode, second order nonlinear mode, and the third order nonlinear mode, respectively. Similarly, X13, X23, and X33 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,k1,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 m=1. Here: the underlying nonlinear dynamical system is given by the following twelve ODEs:
(34)(X˙1m1X˙1m2X˙1m3X˙1m4X˙1m5X˙1m6)=(X1m2X1m4X1m5X1m6Lm1+2mX1m1X1m5Lm2+Nm)
for m=0,1 where km:=(m+1)k1; Lm1:=-X1m1+(km2+ωm)X1m3,
(35)Lm2:=-(km4+ωmkm2Pr)X1m1+(2km2+ωmPr)X1m4+km2RaX1m3,Nm:=2mPr(X101X106+X102X104-km2X101X102),
and the twelve boundary conditions are given by
(36)(X1m1,X1m2X1m3)(0)=(0,0,0)=(X1m1,X1m2X1m3)(1),
for each m=0,1. The fixed parametric values are taken as the following
(37)(ω0,ω1,Pr,k1,Ra)=(85.80,0,10,1.603,15874).In this case, the missing boundary conditions forz=0 using an initial guess of c=0 in the rectangle R defined by R:={x∈ℝ6∣∥x-c∥<2200} are found to be the following:
(38)(X104X105X106X114X115X116)(0,c*)=c*=(+79.20602598022311+0.045512918963491-500.0000000000000-129.3119890195200-0.119928096846000+1357.197011463896).Here the Newton iterates are also found to converge rapidly. Thenonlinear 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.Square pattern in the Rayleigh-Bénard convection for(ω0,ω1,Pr,k1,Ra)=(85.80,0,10,1.603,15874).
(a)
(b)Square pattern (side view inxz-plane) in the Rayleigh-Bénard convection for (ω0,ω1,Pr,k1,Ra)=(85.80,0,10,1.603,15874).
(a)
(b)


## 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 (Theorem6) 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.



---

*Source: 650208-2013-08-22.xml*