# A Third-Orderp-Laplacian Boundary Value Problem Solved by an SL(3,ℝ) Lie-Group Shooting Method

Authors: Chein-Shan Liu
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The boundary layer problem for power-law fluid can be recast to a third-orderp-Laplacian boundary value problem (BVP). In this paper, we transform the third-order p-Laplacian into a new system which exhibits a Lie-symmetry SL(3,ℝ). Then, the closure property of the Lie-group is used to derive a linear transformation between the boundary values at two ends of a spatial interval. Hence, we can iteratively solve the missing left boundary conditions, which are determined by matching the right boundary conditions through a finer tuning of r∈[0,1]. The present SL(3,ℝ) Lie-group shooting method is easily implemented and is efficient to tackle the multiple solutions of the third-order p-Laplacian. When the missing left boundary values can be determined accurately, we can apply the fourth-order Runge-Kutta (RK4) method to obtain a quite accurate numerical solution of the p-Laplacian.


---

## Body

## 1. Introduction

The power-law fluids have been called the Ostwald-de Waele fluids which have been well examined in the past several decades, because the constitutive equation of such a fluid not only gives a good expression for a large portion of the non-Newtonian fluids but also encompasses a Newtonian fluid as well. The theoretical boundary layer theory for power-law fluids was first investigated by Schowalter [1], and then Acrivos et al. [2] obtained a similarity solution. The experimental results that a significant drag reduction can be achieved by injecting fluid into the boundary layer motivated the investigations of the non-Newtonian boundary layer flows with injection or suction at the surface. Flows with suction or injection through a porous wall are of practical interest for cooling, delaying transition to turbulence, and the prevention of separation in an adverse pressure gradient.The drag force inside the shear layer is a consequence of pressure distribution on the surface. Realizing the nature of this force by a mathematical modelling to predict the drag force and the associated behavior of fluid flow has been the focus of considerable research. The reason for the interest in the analysis of the boundary layer flows along solid surfaces is the possibility of applying the theory to the efficient design of supersonic and hypersonic flights. Besides, the mathematical model considered in the present research has importance in studying many problems of engineering, meteorology, and oceanography, for example, Howell et al. [3], Nachman and Callegari [4], Ozisik [5], Schlichting [6], Shu and Wilks [7], and Zheng and Zhang [8].We assume that the moving flat plate is semi-infinite with a porous surface and that the plate is moving at a constant speedUw in the direction parallel to an oncoming flow with a constant speed U∞. By the assumption of incompressibility and the conservation of momentum, the laminar flow satisfies
(1)∂U∂X+∂V∂Y=0,U∂U∂X+V∂U∂Y=1ρ∂τXY∂Y.
In the above, X and Y are the coordinates attached to the plate in the horizontal and perpendicular directions, and U and V are, respectively, the velocity components of the flow in the X and Y directions. The fluid density ρ is assumed to be a constant.The shear stress is governed by a power law(2)τXY=K|∂U∂Y|N-1∂U∂Y,
where K>0 is a constant and the power N>0 reflects the discrepancy to the Newtonian fluids (with N=1). The case with N<1 is the power law of pseudoplastic fluids, and N>1 is the dilatant fluids. The corresponding boundary conditions are given by
(3)U(X,0)=Uw,U(X,+∞)=U∞,V(X,0)=Vw(X)=V0X-N/(N+1).After introducing a similarity variable and a stream function(4)η=BXβY,ϕ(X,Y)=AXσf(η)
with
(5)σ=1N+1,β=-σ,B=(ρU∞2-N(N+1)K)1/(N+1),A=U∞B,
we can obtain
(6)(|f′′(η)|N-1f′′(η))′+f(η)f′′(η)=0,
which is subjected to the following boundary conditions:
(7)f(0)=-C0,f′(0)=ξ,f′(+∞)=1.
In the above, ξ=Uw/U∞ is the velocity ratio. When ξ<0, we have a reverse flow attached near the boundary. When 0<ξ<1, the speed of the oncoming fluid is larger than that of the plate. When ξ>1, the speed of the moving plate is faster than the speed of the oncoming fluid. The term C0=(N+1)BV0/U∞ is a constant related to the situation of suction if it is negative or injection if it is positive.Previously, the author and his coworkers have developed the Lie-group shooting method based on the Lorentz-group [9, 10] to solve the boundary layer equations [11–14]. Liu et al. [15] have found the multiple solutions of boundary layer problem by an enhanced fictitious time integration method for the nonlinear algebraic equations discretized from the governing equation. In this paper, we propose a more simple and powerful Lie-group shooting method directly based on the three-dimensional special linear group SL(3,ℝ) to solve the p(=N+1)-Laplacian boundary layer equations with multiple solutions. Liu [16] used the SL(2,ℝ) Lie-group shooting method to solve the eigenvalue problem of the second-order Sturm-Liouville equation. Recently, Liu [17] has successfully applied the SL(3,ℝ) Lie-group shooting method to effectively solve the Falkner-Skan boundary layer problem. The present study is an extension of these researches.The operator(|f′′(η)|N-1f′′(η))′ in (6) is one sort of the third-order p(=N+1)-Laplacian [18]. Besides the non-Newtonian fluid [19–21], the p-Laplacian operator arises in some different physical and mathematical modelling of the combustion theory [22], population dynamics [23, 24], and also the Monge-Kantorovich partial differential equation [25]. There had been many papers for studying the second-order p-Laplacian about its existence of positive solutions. Only a few papers are devoted to the numerical solution of the third-order p-Laplacian.This paper is organized as follows. In Section2, we consider (6) as a special case of the third-order p-Laplacian, and then by a translation of the variable to a new variable with a positive value, we can transform the p-Laplacian into a new system, which exhibits a Lie-symmetry of SL(3,ℝ). We introduce some mathematical requirements of the Lie-group formulation of the resulting ODEs and construct a Lie-group shooting method based on the Lie-group SL(3,ℝ). The boundary conditions that the present SL(3,ℝ) Lie-group shooting method can be applied to are discussed. In Section 3, we test the performance of the newly developed SL(3,ℝ) Lie-group shooting method for several examples, and in Section 4, the multiple solutions of the boundary layer problem of power-law fluids are investigated. Finally, we draw conclusions in Section 5.


## 2. AnSL(3,ℝ) Lie-Group Shooting Method

Equation (6) is a special case of the following p-Laplacian:
(8)(ϕp(y′′(x)))′=F(x,y(x),y′(x),y′′(x)),x∈(x0,xf),
of which the boundary conditions have many types to be discussed below. In the above, ϕp(s)=|s|p-2s,p>1,ϕp-1=ϕq, where 1/p+1/q=1.In order to develop anSL(3,ℝ) Lie-group shooting method, we suppose that y>-∞, such that there exists a constant k0 rendering
(9)u(x)=y(x)+k0>0,∀x∈[x0,xf].
Then, (8) becomes
(10)(ϕp(u′′(x)))′=F(x,u(x)-k0,u′(x),u′′(x))=:H(x,u(x),u′(x),u′′(x)).


### 2.1. A Group-Preserving Scheme

Let(11)v(x)=ϕp(u′′(x)),
and thus by ϕp-1=ϕq, we have
(12)u′′(x)=ϕq(v(x)).
At the same time, from (10) and (11), it follows that
(13)v′(x)=H(x,u(x),u′(x),ϕq(v(x)))=:f(x,u(x),u′(x),v(x)).Letw(x)=u′(x), and we have w′(x)=ϕq(v(x)) by (12). Now, let
(14)u1=u,u2=w=u′,u3=v.
The above equations u′(x)=w(x),w′(x)=ϕq(v(x)), and (13) can be written as a system of three first-order ordinary differential equations (ODEs) by
(15)ddx[u1(x)u2(x)u3(x)]=[010ϕq(u3(x))u100f(x,u1(x),u2(x),u3(x))u100]×[u1(x)u2(x)u3(x)],
which is well defined because of u1=u>0 by (9).The differential equations system (15) is highly nonlinear due to the appearance of f(x,u1,u2,u3)/u1 and ϕq(u3(x))/u1 in the coefficient matrix; however, it allows a Lie-symmetry SL(3,ℝ):
(16)ddxG(x)=A(x)G(x),G(x0)=I3,(17)A(x)=[010ϕqu100fu100],(18)detG(x)=1,
because of trA=0. Here, for saving notations, we use f=f(x,u1,u2,u3) and ϕq=ϕq(u3). The above G is a fundamental matrix of (15).Accordingly, we can develop a group-preserving scheme (GPS) to solve (15)
(19)Uk+1=G(k)Uk,
where Uk denotes the numerical value of U=(u1,u2,u3)T at the discrete space xk and G(k)∈SL(3,ℝ). This Lie-symmetry is known as the three-dimensional real-valued special linear group, denoted by SL(3,ℝ).


### 2.2. A Generalized Midpoint Rule

Applying the GPS in (19) to (15) with an initial condition U(x0)=U0, we can find U(x). Assuming that the stepsize used in the GPS is Δx=ℓ/K, where ℓ=xf-x0, we can calculate the value of U at x=xf by
(20)Uf=GK(Δx)⋯G1(Δx)U0.Now, we prove the following closure property of the Lie-groupSL(3,ℝ):
(21)G1(x),G2(x)∈SL(3,ℝ)⇒G2(x)G1(x)∈SL(3,ℝ).
By the assumptions of G1(x)∈SL(3,ℝ) and G2(x)∈SL(3,ℝ), we have detG1(x)=1 and detG2(x)=1. Then, by using the following result:
(22)det[G2(x)G1(x)]=detG2(x)detG1(x),
it is straightforward to verify that det[G2(x)G1(x)]=1, which means that G2(x)G1(x)∈SL(3,ℝ). Thus, we have proven (21).Because eachGi,i=1,…,K in (20) is an element of the Lie-group SL(3,ℝ), and by the above closure property of the Lie-group SL(3,ℝ),GK(Δx)⋯G1(Δx) is also a Lie-group element of SL(3,ℝ), denoted by G. Hence, we have
(23)Uf=GU0.
This is a one-step Lie-group transformation from U0 to Uf, acting by G∈SL(3,ℝ).However, it is very hard to obtain an exact solution ofG because the differential equations system is highly nonlinear. Before the derivation of a suitable form for G, let us recall the mean value theorem for a continuous function f(x), which is defined in an interval of x∈[a,b]. The mean value theorem asserts that there exists at least one c∈[a,b], such that the following equality holds
(24)∫ab‍f(x)dx=f(c)[b-a],
where the value of c depends on the function f(x). In terms of the weighting factor r∈[0,1], we can write c=ra+(1-r)b. Therefore, it means that there exists at least one r∈[0,1], such that (24) is satisfied. The above theorem enables us to evaluate the value of the integral in (24) by an area of a rectangle with a width b-a times a height f(c), where f(c) is calculated by a mid-point rule with a suitable c∈[a,b].BecauseG is a solution of (16), we can formally write it by an exponential mapping
(25)G(x)=exp[∫x0x‍A(ξ)dξ].
When A is not a constant matrix, in general we do not have a closed-form solution of G(x). However, motivated by the above mean value theorem and to be a reasonable approximation, we can calculate G in (23) by a generalized mid-point rule, which is obtained from an exponential mapping of A by taking the values of the variables in A at a suitable mid-point: x^=rx0+(1-r)xf, where r∈[0,1] is an unknown constant to be determined by the shooting method. So we can compute this G by
(26)G(r)=exp[(xf-x0)A^],
which is corresponding to a constant matrix
(27)A^=:[010b00a00]=[010ϕ^qu^100f^u^100],
where a=f^/u^1 and b=ϕ^q/u^1 are supposed to be constant.ThisSL(3,ℝ) Lie-group element generated from such a constant matrix A^∈sl(3,ℝ) has a closed-form solution. If b>0 and c=a/b, then we have
(28)G(r)=[G11G12G13G21G22G23G31G32G33]=[cosh(ℓb)1bsinh(ℓb)0bsinh(ℓb)cosh(ℓb)0csinh(ℓb)cb[cosh(ℓb)-1]1],
where ℓ=xf-x0 denotes the length of the interval [x0,xf].Ifb<0 and c=a/-b, then we have
(29)G(r)=[G11G12G13G21G22G23G31G32G33]=[cos(ℓ-b)1-bsin(ℓ-b)0--bsin(ℓ-b)cos(ℓ-b)0csin(ℓ-b)c-b[1-cos(ℓ-b)]1].
In the above, we have taken
(30)x^=rx0+(1-r)xf,u^1=ru10+(1-r)u1f,u^2=ru20+(1-r)u2f,u^3=ru30+(1-r)u3f,f^:=f(x^,u^1,u^2,u^3),ϕ^q:=ϕq(u^3),a=f^u^1,b=ϕ^qu^1,c=a|b|.For the special case ofb=0, we can derive
(31)G(r)=[G11G12G13G21G22G23G31G32G33]=[1ℓ0010aℓaℓ221].


### 2.3. Specification of Boundary Conditions

For the third-order ODEs, there are several different type boundary conditions. In this section, we study this problem that under what type boundary conditions the presentSL(3,ℝ) Lie-group shooting method is applicable.Let(32)U0=[u10u20u30],Uf=[u1fu2fu3f]
denote, respectively, the left-end and right-end boundary values of U=(u1,u2,u3)T. For linear type boundary conditions (separable or nonseparable), we can describe the boundary conditions by the following equation:
(33)B0U0+BfUf=b,
where both B0 and Bf are 3×3 matrices and b∈ℝ3 is a constant vector, which might be zero.Inserting (23) into (33), we have
(34)[B0+BfG]U0=b,
such that for a nonempty solution of U0, we require
(35)det[B0+BfG]≠0.
It means that the matrix B0+BfG must be invertible.In order to demonstrate the above idea about the specification of the boundary conditions, of which the presentSL(3,ℝ) Lie-group shooting method is applicable, let us take p=N+1 and F=-uu′′ and (6) under the following boundary conditions:
(36)u(0)=-C0,u′(0)=ξ,u′(η∞)=1,thatis,u1(0)=-C0,u2(0)=ξ,u2(η∞)=1,
where we have replaced ∞ by a finite number η∞. In terms of (33), we have
(37)B0=[100010000],Bf=[000000010],b=[-C0ξ1].
Then, by (28) or (29), we can obtain
(38)B0+BfG=[100010000]+[000000010][G11G120G21G220G31G321]=[100010G21G220].
Thus, the above matrix B0+BfG is not invertible, because of det[B0+BfG]=0. So we can conclude that for the power-law fluid under the boundary conditions (36), the present method is not applicable. In Section 4, we will give another type approach.Physically, we can specify that at a large distance from the boundary layer the shear stress is quite small. Then, instead of (36), we can specify
(39)u(0)=-C0,u′(0)=ξ,u′′(η∞)=s0,thatis,u1(0)=-C0,u2(0)=ξ,u3(η∞)=ϕN+1(s0),
where s0 is a small number. In terms of (33), we have
(40)B0=[100010000],Bf=[000000001],b=[-C0ξϕN+1(s0)].
Then, by (28) or (29), we can obtain
(41)B0+BfG=[100010000]+[000000001][G11G120G21G220G31G321]=[100010G31G321],
which is invertible, due to det[B0+BfG]=1. Then, the present method is applicable. Below, we discuss the Lie-group shooting solution for the boundary layer problem of the power-law fluid under the boundary conditions (39).


### 2.4. AnSL(3,ℝ) Lie-Group Shooting Method

In order to demonstrate the application of theSL(3,ℝ) Lie-group shooting method to find the missing left boundary conditions, as a representative case, let us take p=N+1 and F=-yy′′, and then we can recover to (6)
(42)(|y′′(x)|N-1y′′(x))′=-y(x)y′′(x),x∈(0,η∞),
which is subjected to the following boundary conditions:
(43)y(0)=-C0,y′(0)=ξ,y′′(η∞)=s0=0.001,
where C0 and ξ are given constants, and we use a large value η∞, say η∞=6, to replace the last boundary condition in (6), and also y′(η∞)=1 is changed to y′′(η∞)=s0 as just mentioned above.The stepping technique developed for solving the initial value problem (IVP) requires both the initial conditions ofu1=u,u2=u′, and u3=v for the third-order ODEs. Starting from the initial values of u1,u2, and u3, we can numerically integrate the following IVP step by step from x=x0 to x=xf:
(44)ddx[u1(x)u2(x)u3(x)]=[010ϕq(u3(x))u100f(x,u1(x),u2(x),u3(x))u100]×[u1(x)u2(x)u3(x)],(45)u1(0)=u10,(46)u2(0)=u20,(47)u3(0)=u30,
where some unknown initial values are to be found by the SL(3,ℝ) Lie-group shooting method.In (43), u10=-C0 and u20=ξ are given, but u30 is an unknown constant to be determined such that we can satisfy the target equation of u3f=ϕN+1(s0). Starting from an initial guess of u1f,u3f, and u30, we can solve the unknown initial value u30 by the following iterative processes:
(48)u30=u3f-G31u10-G32u20G33,u1f=G11u10+G12u20+G13u30,u2f=G21u10+G22u20+G23u30,
which are obtained from (23). Inserting the initially guessed values of u1f,u3f,  and u30 and the given values of u10=-C0,u20=ξ, and u3f=ϕN+1(s0) into (28) or (29) with a specified r∈[0,1], we can evaluate Gij,i,j=1,2,3, and then by (48), we can generate the new values of u30,u1f, and u2f, until they are convergent. If the new values of u30,u1f, and u2f converge to satisfy the following convergence criterion:
(49)([u2f(k+1)-u2f(k)]2[u30(k+1)-u30(k)]2+[u1f(k+1)-u1f(k)]2+[u2f(k+1)-u2f(k)]2)1/2≤ε,
then the iterations stop. Here, u30(k+1) and u30(k) denote, respectively, the (k+1)th and the kth iteration values of u30. They are defined similarly for u1f and u2f.For a trialr, we can calculate u30 from the above equations by a few iterations and then numerically integrate (44) by the fourth-order Runge-Kutta method (RK4) from 0 to η∞ and compare the end value of u3f with the exact one u2(η∞)=ϕN+1(s0), which is a target equation to be matched. Indeed, we need to find the root of the equation u3f-ϕN+1(s0)=0, where u3f is a numerically integrated result, depending on r. It can be done in practice by adjusting the value of r to a point such that the curve of mismatching error is intersected with the zero line at that point.


## 2.1. A Group-Preserving Scheme

Let(11)v(x)=ϕp(u′′(x)),
and thus by ϕp-1=ϕq, we have
(12)u′′(x)=ϕq(v(x)).
At the same time, from (10) and (11), it follows that
(13)v′(x)=H(x,u(x),u′(x),ϕq(v(x)))=:f(x,u(x),u′(x),v(x)).Letw(x)=u′(x), and we have w′(x)=ϕq(v(x)) by (12). Now, let
(14)u1=u,u2=w=u′,u3=v.
The above equations u′(x)=w(x),w′(x)=ϕq(v(x)), and (13) can be written as a system of three first-order ordinary differential equations (ODEs) by
(15)ddx[u1(x)u2(x)u3(x)]=[010ϕq(u3(x))u100f(x,u1(x),u2(x),u3(x))u100]×[u1(x)u2(x)u3(x)],
which is well defined because of u1=u>0 by (9).The differential equations system (15) is highly nonlinear due to the appearance of f(x,u1,u2,u3)/u1 and ϕq(u3(x))/u1 in the coefficient matrix; however, it allows a Lie-symmetry SL(3,ℝ):
(16)ddxG(x)=A(x)G(x),G(x0)=I3,(17)A(x)=[010ϕqu100fu100],(18)detG(x)=1,
because of trA=0. Here, for saving notations, we use f=f(x,u1,u2,u3) and ϕq=ϕq(u3). The above G is a fundamental matrix of (15).Accordingly, we can develop a group-preserving scheme (GPS) to solve (15)
(19)Uk+1=G(k)Uk,
where Uk denotes the numerical value of U=(u1,u2,u3)T at the discrete space xk and G(k)∈SL(3,ℝ). This Lie-symmetry is known as the three-dimensional real-valued special linear group, denoted by SL(3,ℝ).


## 2.2. A Generalized Midpoint Rule

Applying the GPS in (19) to (15) with an initial condition U(x0)=U0, we can find U(x). Assuming that the stepsize used in the GPS is Δx=ℓ/K, where ℓ=xf-x0, we can calculate the value of U at x=xf by
(20)Uf=GK(Δx)⋯G1(Δx)U0.Now, we prove the following closure property of the Lie-groupSL(3,ℝ):
(21)G1(x),G2(x)∈SL(3,ℝ)⇒G2(x)G1(x)∈SL(3,ℝ).
By the assumptions of G1(x)∈SL(3,ℝ) and G2(x)∈SL(3,ℝ), we have detG1(x)=1 and detG2(x)=1. Then, by using the following result:
(22)det[G2(x)G1(x)]=detG2(x)detG1(x),
it is straightforward to verify that det[G2(x)G1(x)]=1, which means that G2(x)G1(x)∈SL(3,ℝ). Thus, we have proven (21).Because eachGi,i=1,…,K in (20) is an element of the Lie-group SL(3,ℝ), and by the above closure property of the Lie-group SL(3,ℝ),GK(Δx)⋯G1(Δx) is also a Lie-group element of SL(3,ℝ), denoted by G. Hence, we have
(23)Uf=GU0.
This is a one-step Lie-group transformation from U0 to Uf, acting by G∈SL(3,ℝ).However, it is very hard to obtain an exact solution ofG because the differential equations system is highly nonlinear. Before the derivation of a suitable form for G, let us recall the mean value theorem for a continuous function f(x), which is defined in an interval of x∈[a,b]. The mean value theorem asserts that there exists at least one c∈[a,b], such that the following equality holds
(24)∫ab‍f(x)dx=f(c)[b-a],
where the value of c depends on the function f(x). In terms of the weighting factor r∈[0,1], we can write c=ra+(1-r)b. Therefore, it means that there exists at least one r∈[0,1], such that (24) is satisfied. The above theorem enables us to evaluate the value of the integral in (24) by an area of a rectangle with a width b-a times a height f(c), where f(c) is calculated by a mid-point rule with a suitable c∈[a,b].BecauseG is a solution of (16), we can formally write it by an exponential mapping
(25)G(x)=exp[∫x0x‍A(ξ)dξ].
When A is not a constant matrix, in general we do not have a closed-form solution of G(x). However, motivated by the above mean value theorem and to be a reasonable approximation, we can calculate G in (23) by a generalized mid-point rule, which is obtained from an exponential mapping of A by taking the values of the variables in A at a suitable mid-point: x^=rx0+(1-r)xf, where r∈[0,1] is an unknown constant to be determined by the shooting method. So we can compute this G by
(26)G(r)=exp[(xf-x0)A^],
which is corresponding to a constant matrix
(27)A^=:[010b00a00]=[010ϕ^qu^100f^u^100],
where a=f^/u^1 and b=ϕ^q/u^1 are supposed to be constant.ThisSL(3,ℝ) Lie-group element generated from such a constant matrix A^∈sl(3,ℝ) has a closed-form solution. If b>0 and c=a/b, then we have
(28)G(r)=[G11G12G13G21G22G23G31G32G33]=[cosh(ℓb)1bsinh(ℓb)0bsinh(ℓb)cosh(ℓb)0csinh(ℓb)cb[cosh(ℓb)-1]1],
where ℓ=xf-x0 denotes the length of the interval [x0,xf].Ifb<0 and c=a/-b, then we have
(29)G(r)=[G11G12G13G21G22G23G31G32G33]=[cos(ℓ-b)1-bsin(ℓ-b)0--bsin(ℓ-b)cos(ℓ-b)0csin(ℓ-b)c-b[1-cos(ℓ-b)]1].
In the above, we have taken
(30)x^=rx0+(1-r)xf,u^1=ru10+(1-r)u1f,u^2=ru20+(1-r)u2f,u^3=ru30+(1-r)u3f,f^:=f(x^,u^1,u^2,u^3),ϕ^q:=ϕq(u^3),a=f^u^1,b=ϕ^qu^1,c=a|b|.For the special case ofb=0, we can derive
(31)G(r)=[G11G12G13G21G22G23G31G32G33]=[1ℓ0010aℓaℓ221].


## 2.3. Specification of Boundary Conditions

For the third-order ODEs, there are several different type boundary conditions. In this section, we study this problem that under what type boundary conditions the presentSL(3,ℝ) Lie-group shooting method is applicable.Let(32)U0=[u10u20u30],Uf=[u1fu2fu3f]
denote, respectively, the left-end and right-end boundary values of U=(u1,u2,u3)T. For linear type boundary conditions (separable or nonseparable), we can describe the boundary conditions by the following equation:
(33)B0U0+BfUf=b,
where both B0 and Bf are 3×3 matrices and b∈ℝ3 is a constant vector, which might be zero.Inserting (23) into (33), we have
(34)[B0+BfG]U0=b,
such that for a nonempty solution of U0, we require
(35)det[B0+BfG]≠0.
It means that the matrix B0+BfG must be invertible.In order to demonstrate the above idea about the specification of the boundary conditions, of which the presentSL(3,ℝ) Lie-group shooting method is applicable, let us take p=N+1 and F=-uu′′ and (6) under the following boundary conditions:
(36)u(0)=-C0,u′(0)=ξ,u′(η∞)=1,thatis,u1(0)=-C0,u2(0)=ξ,u2(η∞)=1,
where we have replaced ∞ by a finite number η∞. In terms of (33), we have
(37)B0=[100010000],Bf=[000000010],b=[-C0ξ1].
Then, by (28) or (29), we can obtain
(38)B0+BfG=[100010000]+[000000010][G11G120G21G220G31G321]=[100010G21G220].
Thus, the above matrix B0+BfG is not invertible, because of det[B0+BfG]=0. So we can conclude that for the power-law fluid under the boundary conditions (36), the present method is not applicable. In Section 4, we will give another type approach.Physically, we can specify that at a large distance from the boundary layer the shear stress is quite small. Then, instead of (36), we can specify
(39)u(0)=-C0,u′(0)=ξ,u′′(η∞)=s0,thatis,u1(0)=-C0,u2(0)=ξ,u3(η∞)=ϕN+1(s0),
where s0 is a small number. In terms of (33), we have
(40)B0=[100010000],Bf=[000000001],b=[-C0ξϕN+1(s0)].
Then, by (28) or (29), we can obtain
(41)B0+BfG=[100010000]+[000000001][G11G120G21G220G31G321]=[100010G31G321],
which is invertible, due to det[B0+BfG]=1. Then, the present method is applicable. Below, we discuss the Lie-group shooting solution for the boundary layer problem of the power-law fluid under the boundary conditions (39).


## 2.4. AnSL(3,ℝ) Lie-Group Shooting Method

In order to demonstrate the application of theSL(3,ℝ) Lie-group shooting method to find the missing left boundary conditions, as a representative case, let us take p=N+1 and F=-yy′′, and then we can recover to (6)
(42)(|y′′(x)|N-1y′′(x))′=-y(x)y′′(x),x∈(0,η∞),
which is subjected to the following boundary conditions:
(43)y(0)=-C0,y′(0)=ξ,y′′(η∞)=s0=0.001,
where C0 and ξ are given constants, and we use a large value η∞, say η∞=6, to replace the last boundary condition in (6), and also y′(η∞)=1 is changed to y′′(η∞)=s0 as just mentioned above.The stepping technique developed for solving the initial value problem (IVP) requires both the initial conditions ofu1=u,u2=u′, and u3=v for the third-order ODEs. Starting from the initial values of u1,u2, and u3, we can numerically integrate the following IVP step by step from x=x0 to x=xf:
(44)ddx[u1(x)u2(x)u3(x)]=[010ϕq(u3(x))u100f(x,u1(x),u2(x),u3(x))u100]×[u1(x)u2(x)u3(x)],(45)u1(0)=u10,(46)u2(0)=u20,(47)u3(0)=u30,
where some unknown initial values are to be found by the SL(3,ℝ) Lie-group shooting method.In (43), u10=-C0 and u20=ξ are given, but u30 is an unknown constant to be determined such that we can satisfy the target equation of u3f=ϕN+1(s0). Starting from an initial guess of u1f,u3f, and u30, we can solve the unknown initial value u30 by the following iterative processes:
(48)u30=u3f-G31u10-G32u20G33,u1f=G11u10+G12u20+G13u30,u2f=G21u10+G22u20+G23u30,
which are obtained from (23). Inserting the initially guessed values of u1f,u3f,  and u30 and the given values of u10=-C0,u20=ξ, and u3f=ϕN+1(s0) into (28) or (29) with a specified r∈[0,1], we can evaluate Gij,i,j=1,2,3, and then by (48), we can generate the new values of u30,u1f, and u2f, until they are convergent. If the new values of u30,u1f, and u2f converge to satisfy the following convergence criterion:
(49)([u2f(k+1)-u2f(k)]2[u30(k+1)-u30(k)]2+[u1f(k+1)-u1f(k)]2+[u2f(k+1)-u2f(k)]2)1/2≤ε,
then the iterations stop. Here, u30(k+1) and u30(k) denote, respectively, the (k+1)th and the kth iteration values of u30. They are defined similarly for u1f and u2f.For a trialr, we can calculate u30 from the above equations by a few iterations and then numerically integrate (44) by the fourth-order Runge-Kutta method (RK4) from 0 to η∞ and compare the end value of u3f with the exact one u2(η∞)=ϕN+1(s0), which is a target equation to be matched. Indeed, we need to find the root of the equation u3f-ϕN+1(s0)=0, where u3f is a numerically integrated result, depending on r. It can be done in practice by adjusting the value of r to a point such that the curve of mismatching error is intersected with the zero line at that point.


## 3. Numerical Examples

Example 1.
First, we consider the followingp-Laplacian:
(50)(ϕp(y′′(x)))′=h(x)+y(x),x∈(0,1),y(0)=0,y′(1)=π,y′′(0)=0.
We assume that the closed-form solution is y(x)=-sin(πx). Hence,
(51)h(x)=π3(p-1)[π2sin(πx)]p-2cos(πx)+sin(πx).
We can use the following equations to iteratively solve the unknown initial value ofu20:
(52)u20=u2f-G21u10-G23u30G22,u1f=G11u10+G12u20+G13u30,u3f=G31u10+G32u20+G33u30.
We takek0=3.5 for p=3 and k0=4 for p=3.5 and use the SL(3,ℝ) Lie-group shooting method developed in Section 2.4 by adjusting the value of u20. If the target equation u′(1)=π is satisfied, then we obtain the numerical solution.
The convergence criterion isε=10-8. Although under this stringent convergence criterion the iteration process to find u20 is convergent very fast as shown in Figure 1(a), where for r∈[0.58,0.59] the iteration numbers are between 16 and 23 for the case p=3. In Figure 1(a), we also plot the mismatching errors with respect to r in a range [0.58,0.59] for p=3, while in a range [0.54,0.55] for p=3.5. Both have an intersection point with the zero line. Then, through a finer tuning of the value to r=0.58135439224 for the case p=3 and to r=0.547436304748 for the case p=3.5, we can match the right-end boundary condition very precisely with an error in the order 10-11. The numerical solutions of y and y′ are, respectively, plotted in Figures 1(b) and 1(c), which are almost coincident with the closed-form solutions. Therefore, we plot the numerical errors, which are the absolute differences between exact solutions and numerical solutions, in Figure 2 for p=3 and p=3.5. It can be seen that the numerical results are quite accurate.Example1 solved by the SL(3,ℝ) method, (a) showing the error of mismatching and iterations number and comparing numerical solutions and exact solutions (b) of y and (c) y′.
(a)
(b)
(c)Example1 with p=3 and p=3.5, showing the numerical errors of (a) y, (b) y′, and (c) y′′.
(a)
(b)
(c)Example 2.
Then, we consider a different case of (50) by
(53)(ϕp(y′′(x)))′=h(x)-y2(x),x∈(0,1),y(0)=0,y(1)=0,y′′(0)=0.
Similarly, the closed-form solution is y(x)=-sin(πx). Hence,
(54)h(x)=π3(p-1)[π2sin(πx)]p-2cos(πx)+sin2(πx).
We can use the following equations to iteratively solve the unknown initial value ofu20:
(55)u20=u1f-G11u10-G13u30G12,u2f=G21u10+G22u20+G23u30,u3f=G31u10+G32u20+G33u30.
We takek0=3.5 and p=3. If the target equation u(1)=k0 is satisfied, then we obtain the numerical solution. In Figure 3(a), we plot the mismatching error and the number of iterations with respect to r in a range [0.575,0.585]. Then through a finer tuning of the value to r=0.5813782916, we can match the right-end boundary condition very precisely with an error in the order 10-10. Upon comparing the numerical solutions with the closed-form solutions, the numerical errors of y and y′ are, respectively, plotted in Figures 3(b) and 3(c). It can be seen that for y the accuracy is in the order of 10-9, while that for y′ the accuracy is in the order of 10-6.Example2 solved by the SL(3,ℝ) method, (a) showing the error of mismatching and iterations number and showing the numerical errors (b) of y and (c) y′.
(a)
(b)
(c)Example 3.
We consider the same equation (53) but under the following boundary conditions:
(56)y(0)=0,y′′(0)=1,y′′(1)=1.
When the closed-form solution is given by
(57)y(x)=x22-sin(πx),
the term h(x) is given by
(58)h(x)=π3(p-1)[1+π2sin(πx)]p-2cos(πx)+[x22-sin(πx)]2.
For the above boundary conditions, we can use the following equations to iteratively solve the unknown initial value ofu20:
(59)u20=u3f-G31u10-G33u30G32,u1f=G11u10+G12u20+G13u30,u2f=G21u10+G22u20+G23u30.
We takek0=3 and p=1.5. If the target equation u′′(1)=k0+1 is satisfied, then we obtain the numerical solution. When we plot the mismatching error with respect to r in a range [0.4675,0.4755] in Figure 4, we find that there exist two intersection points at r=0.467958888 and r=0.47500045, which means that there exist two solutions. In (57), we only give one exact solution, but we do not have another solution as given in a closed-form.
Then, through a finer tuning of the value ofr, we can match the right-end boundary condition very precisely with the error in the order 10-8, and the first numerical solution is obtained with r=0.467958888, while the second numerical solution is obtained with r=0.47500045. Upon comparing the first numerical solution with the closed-form solution in (57), the numerical errors of y,y′, and y′′ are, respectively, plotted in Figures 5(a), 5(b), and 5(c). It can be seen that all the accuracies are in the order of 10-8. In Figure 6, we compare the first numerical solution and the second numerical solution with the exact one. It can be seen that when the first numerical solution is almost coincident with the exact solution, the second numerical solution is obviously different from the first numerical solution. We can also observe that the second numerical solution satisfies the boundary conditions in (56) very precisely.Figure 4
Example3 solved by the SL(3,ℝ) method, showing the error of mismatching. There are two roots.Example3 showing the errors of the first numerical solution by comparing with the exact solution.
(a)
(b)
(c)Example3 comparing the first numerical solution and the second numerical solution with the exact solution.
(a)
(b)
(c)Example 4.
Then, we consider a more general boundary conditions for the followingp-Laplacian:
(60)(ϕp(y′′(x)))′=F(x,y(x)),x∈(0,π),a11y(0)+a12y′(0)=B1,a21y(π)+a22y′(π)=B2,y′′(0)=0.
Similarly, we consider a translation with u(x)=y(x)+k0>0, such that we have
(61)(ϕp(u′′(x)))′=f(x,u(x))≔F(x,u(x)-k0),x∈(0,π),a11u(0)+a12u′(0)=b1:=a11k0+B1,a21u(π)+a22u′(π)=a21k0+B2,u′′(0)=0.
In terms of (33), we can write
(62)(B0+BfG)U0=([a11a120000001]+[000a21a220000][G11G120G21G220G31G321])×[u10u20u30]=[a11a120a21G11+a22G21a21G12+a22G220001][u10u20u30]=[b1b20].
Thus, from (62) and (23), we can solve
(63)u10=b1(a21G12+a22G22)-b2a12a11(a21G12+a22G22)-a12(a21G11+a22G21),u20=b2a11-b1(a21G11+a22G21)a11(a21G12+a22G22)-a12(a21G11+a22G21),u1f=G11u10+G12u20+G13u30,u2f=G21u10+G22u20+G23u30,u3f=G31u10+G32u20+G33u30.
The above five equations can be used to iteratively solve the five unknowns of u10,u20,u1f,u2f, and u3f. We note that u30=0 by the last boundary condition in (61).
As a demonstrative case, we takeF(x,y(x))=h(x)-y(x), that is, f(x,u(x))=h(x)-[u(x)-k0],p=3, and y(x)=excosx to be the exact solution, where h(x) can be computed from (60) by inserting the above y(x). The boundary conditions are given by
(64)y(0)+y′(0)=2,y(π)+y′(π)=-2eπ,y′′(0)=0,thatis,u(0)+u′(0)=k0+2,u(π)+u′(π)=k0-2eπ,u′′(0)=0.
We takek0=25. If the target equation u1f+u2f-k0+2eπ=0 is satisfied, then we obtain the numerical solution. When we plot the mismatching error with respect to r in a finer range [0.7102,0.7103] in Figure 7, we find that there exists one intersection point at r=0.710283575975. We can match the right-end boundary condition very precisely with an error being -3.694×10-10. Because there are many equations to be solved iteratively, the number of iterations as shown in Figure 7(a) is between 45 and 48, which is higher than the previous three examples. In Figure 7(b), we compare the numerical solution of y(x) with the exact solution y(x)=excosx, whose numerical error as shown in Figure 7(c) is quite accurate in the order of 10-6.
Alternatively, we consider a nonlinear perturbation of the above example under the same boundary conditions but with(65)F(x,y(x))=h(x)-y(x)+12y2(x),h(x)=-2(p-1)ex(cosx+sinx)|2exsinx|p-2,
where we also fix p=3. For this problem, we do not have a closed-form solution. However, we take k0=26, and by taking r=0.6863523 and r=0.799, we can obtain two numerical solutions as shown in Figure 8. The numerical solution as shown by the dashed line is quite unstable. For the purpose of comparison, we also plot the numerical solutions obtained in the last example in (64) by the dashed-dotted lines in Figure 8. It can be seen that the solid lines are somewhat perturbed from the ones of the dashed-dotted lines, but the unstable ones are quite different from the above two solutions.Example4 solved by the SL(3,ℝ) method, (a) showing the error of mismatching and iterations number, (b) comparing the numerical and exact solutions, and (c) showing the numerical error of y.
(a)
(b)
(c)A nonlinear perturbation of Example4, displaying two solutions of (a) y, (b) y′, and (c) y′′.
(a)
(b)
(c)


## 4. Power-Law Fluids

In this section, we consider the boundary layer problems of power-law fluid in (6). We use the SL(3,ℝ) Lie-group shooting method developed in Section 2.3 by adjusting the value of u30. However, before that we need to treat the difficulty mentioned in Section 2.3 if we use the boundary conditions in (36).First, we need to point out that for this boundary layer problem we only consider that the functionf(x) is convex, that is, f′′(x)≥0. So the term b defined in (30) is positive. Then, from (28), (23), and (32), we have
(66)u1f=cosh(ℓb)u10+1bsinh(ℓb)u20,(67)u2f=bsinh(ℓb)u10+cosh(ℓb)u20,(68)u3f=csinh(ℓb)u10+cb[cosh(ℓb)-1]u20+u30,
where ℓ=η∞.Letz:=ℓb. From (67), we can derive a scalar equation to solve z:
(69)u20ℓcoshz+u10zsinhz-ℓu2f=0.
Because u10=k0-C0,u20=ξ, and u2f=1 are given from the boundary conditions, and ℓ=η∞ is selected, we can apply the Newton method to solve the above equation, whose solution is denoted by z0. Hence, we have
(70)b=(z0ℓ)2,
by the definition of z=ℓb. Furthermore, by (30), we have
(71)ru30+(1-r)u3f=ϕp(bu^1).
From (68) and (71), we can solve
(72)u30=ϕp(bu^1)-(1-r)×[csinh(ℓb)u10+cb[cosh(ℓb)-1]u20],u3f=ϕp(bu^1)+r[csinh(ℓb)u10+cb[cosh(ℓb)-1]u20],
where c=-(u^1-k0)ϕq(u^3)/(bu^1). When b is solved, u1f is determined by (66); hence, (72) can be used iteratively to solve the two unknowns of u30 and u3f. We can satisfy the target equation u2f=1 by selecting the best value of r.Example 5.
We fixN=0.8,C0=0,ξ=0.2,η∞=10, and k0=1, and the convergence criterion is ε=10-5. Although under this stringent convergence criterion the iteration process to find u30 is convergent very fast as shown in Figure 9(a), where for r∈[0.65,0.7] the iteration numbers are all to be four. In Figure 9(a), we plot the mismatching error with respect to r in the same range. It can be seen that the mismatching error curve is intersected with the zero line at a point near to 0.68. Then, through a finer tuning of the value to r=0.689064, we can match the right-end boundary condition very precisely with an error in the order of 10-7. The unknown initial value of u3(0)=v(0)=0.4998652 (or f′′(0)=0.4203065) is obtained. The numerical results of f,f′, and f′′ are, respectively, plotted in Figures 9(b)–9(d).The boundary layer problem of a power-law fluid withN=0.8 solved by the SL(3,ℝ) method, (a) showing the error of mismatching and iterations number and showing the numerical solutions of (b) f, (c) f′, and (d) f′′.
(a)
(b)
(c)
(d)Example 6.
We fixN=1.2,C0=0.2,ξ=-0.2,η∞=15, and k0=1, and the convergence criterion is ε=10-5. As shown in Figure 10(a), where for r∈[0.2,0.8] the iteration numbers are all to be four, the mismatching error curve is intersected with the zero line at two points a and b. Then, through a finer tuning of the values to r=0.2408982 and r=0.7191693, we can match the right-end boundary condition very precisely with the errors in the order of 10-8, and thus we obtain two numerical solutions as compared in Figures 10(b)–10(d). For the first solution, the unknown initial value of u3(0)=v(0)=0.18159807 (or f′′(0)=0.24131965) is obtained, while for the second solution the unknown initial value of u3(0)=v(0)=0.024415987 (or f′′(0)=0.04533098) is obtained. The numerical results of v(0)=0.18159807 and v(0)=0.024415987 are very close to that obtained by Liu [14].The boundary layer problem of power-law fluid withN=1.2 solved by the SL(3,ℝ) method, (a) showing the error of mismatching and iterations number and showing two different numerical solutions of (b) f, (c) f′, and (d) f′′.
(a)
(b)
(c)
(d)Example 7.
We fixN=0.3,C0=0.2,ξ=-0.2,η∞=10, and k0=0.5. As shown in Figure 11(a), the mismatching error curve is intersected with the zero line at three points a, b, and c. Then, through a finer tuning of the values to r=0.081054,r=0.6812354, and r=0.96791285, we can obtain three corresponding numerical solutions as compared in Figures 11(b)–11(d). For the first solution, the unknown initial value of f′′(0)=0.150134848 is obtained, and for the second solution the unknown initial value is f′′(0)=0.0337905, while that for the third solution the unknown initial value is f′′(0)=0.01283368. For the last solution, f′′ grows rapidly after x=8.The boundary layer problem of a power-law fluid withN=0.3 solved by the SL(3,ℝ) method, (a) showing the error of mismatching and iterations number and showing three different numerical solutions of (b) f, (c) f′, and (d) f′′.
(a)
(b)
(c)
(d)


## 5. Conclusions

In the present paper, we have offered a rather accurate and simple method with only a few iterations to find the unknown left boundary conditions by applying theSL(3,ℝ) Lie-group shooting method to the third-order p-Laplacian boundary value problems. Also, as an application, we have solved the boundary layer problems of power-law fluids by the present method. The SL(3,ℝ) Lie-group shooting method allows us to express the missing left-end boundary conditions by the closed-form functions of r∈[0,1], where the best r is determined iteratively by matching the right-end boundary conditions. Because the iterations to find the missing left-end boundary conditions are convergent very fast, the Lie-group shooting method based on SL(3,ℝ) is quite computationally efficient. The new method was effective to find the multiple solutions, although for the highly nonlinear case with multiple unknown left boundary conditions.



---

*Source: 497863-2013-03-28.xml*