Analysis and numerical approximation of singular boundary value problems with the p -Laplacian in fluid mechanics

Journal of Computational and Applied Mathematics 262 (2014) 87–104

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Analysis and numerical approximation of singular boundary value problems with the p-Laplacian in fluid mechanics G.Yu. Kulikov a,∗ , P.M. Lima a , M.L. Morgado b a

CEMAT, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal

b

CM-UTAD and Dept Matemática, Univ. Trás-os-Montes e Alto Douro, Apartado 1013, 5000-311 Vila Real, Portugal

article

abstract

info

Article history: Received 13 October 2012 Received in revised form 13 June 2013

This paper studies a generalization of the Cahn–Hilliard continuum model for multi-phase fluids where the classical Laplacian has been replaced by a degenerate one (i.e., the socalled p-Laplacian). The solution’s asymptotic behavior is analyzed at two singular points; namely, at the origin and at infinity. An efficient technique for treating such singular boundary value problems is presented, and results of numerical integration are discussed and compared with earlier computed data. © 2013 Elsevier B.V. All rights reserved.

MSC: 65L05 34B16 Keywords: Singular boundary value problems Nonlinear ordinary differential equations Degenerate Laplacian Shooting method Nested implicit Runge–Kutta formulas with global error control

1. Introduction Mathematical models of phase transitions have been the topic of analytical and computational investigations for a long time. Particular attention has been given to the study of the interface between phases, for example, gaseous bubbles in a liquid or liquid droplets in a gas. According to [1], the Cahn–Hilliard theory [2] is the simplest continuum model for multiphase fluids. Denoting the mass density of a nonhomogeneous fluid by ρ , the above theory assumes that the volume free energy E of the fluid is the sum of the volume energy W (ρ) and a term that takes into account the nonhomogeneity of the fluid; i.e, E (ρ, grad(ρ)) = W (ρ) +

γ 2

(grad(ρ))2

(1)

where γ is a known constant. In the most common models the function W (ρ) is a two-well function with two local minima, which are further referred to as ρl and ρv . Because of the shape of W (ρ), the fluid tends to divide into two phases with the densities ρ = ρl and ρ = ρv . Also, the term γ /2(grad(ρ))2 tends to reduce the variation of the field ρ , turning the interface into a thin layer and endowing it with an energy, the surface tension (see [3,1]). In contrast to what happens in other fields of Physics, in the theory of fluid mixtures, authors often use units such that the fluid density may take either negative or positive values (typically it is positive in one of the phases and negative in the other). Having expressed the

∗

Corresponding author. Tel.: +351 217743097. E-mail addresses: [email protected], [email protected] (G.Yu. Kulikov), [email protected] (P.M. Lima), [email protected] (M.L. Morgado).

0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.09.071

88

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

volume free energy in the form (1), the density profile ρ(r ) can be obtained in the stationary case by means of minimization of the functional J (ρ) =

 

γ

Ω

2

W (ρ(r )) +

 (grad(ρ(r )))2 dr

(2)

where Ω ⊂ RN . This minimization leads to the partial differential equation (PDE)

γ 1ρ = µ(ρ) − µ0

(3)

where µ(ρ) := dW (ρ)/dρ is called the chemical potential and µ0 is a constant. Using spherical symmetry, Eq. (3) can be reduced to the second order ordinary differential equation (ODE) r 1−N r N −1 ρ ′ (r )



′

= f (ρ(r )) ,

r > 0,

(4)

where f is a given function (typically, this is a cubic polynomial with three real zeros); N > 1 is the space dimension. Formula (4) is referred to as the density profile equation. One searches for a monotone solution of (4) that satisfies the following boundary conditions:

ρ ′ (0) = 0,

lim ρ(r ) = ξ

(5)

r →∞

where ξ is a given real number (whose physical meaning is the liquid density). The differential problem (4), (5) is an important topic of research, both from the point of view of qualitative analysis and numerical approximation (see, for instance, [4–9]). In the present paper we deal with a more general model of phase transitions where the volume free energy is presented, instead of (1), by E (ρ) = Wp (ρ) + γ (grad(ρ))p /p with p > 1. Here, Wp (ρ) is also a double potential function, but its form depends on the specific value of p. The introduction of this parameter allows the formation of interfaces of more general form, compared to the classical Cahn–Hilliard theory, to be studied. A mathematical model of such a sort has been analyzed in [10]. There, the authors consider the following energy functional: Jp (ρ) =

  Ω

Wp (ρ(r )) +

γ p

(grad(ρ(r )))

p



dr .

(6)

Note that energy functionals of the form (6) have also been explored in [11,12]. The quantity (grad(ρ))p on the right-hand side of (6) can be regarded as a generalized form of kinetic energy, which reduces to the classical one when p = 2. Having replaced the classical energy functional (2) with the more general formula (6), and choosing a system of units such that µ0 = 0 in (3), one arrives at the PDE

γ 1p ρ = Wp′ (ρ) (7)   where 1p ρ := div |grad(ρ)|p−2 grad(ρ) denotes the so-called degenerate p-Laplacian, which reduces to the classical Laplacian 1ρ when p = 2. By means of a variable separation, the differential problem (7) is again converted to the radial equation

′

r 1−N r N −1 |ρ ′ (r )|p−2 ρ ′ (r )



= fp (ρ(r )) ,

(8)

which is an analogue of the density profile Eq. (4), but for the case of the p-Laplacian. The function fp (ρ) on the right-hand side of (8) is a generalization of the cubic polynomial f (ρ) in (4). It is no longer a polynomial but, as in the classical case p = 2, it has three real zeros physically corresponding to critical values of the fluid density (i.e. to lower and upper bounds of the density and its average density at the interface). As for the classical Laplacian, we study monotone solutions of (8) that satisfy the boundary conditions (5). It is proved in [13] that this boundary value problem (BVP) possesses a unique solution when the right-hand side is a function given by the formula fp (ρ) = −2pλ2 (ξ − ρ)p−1 ρ(ρ + 1)p−1

(9)

where λ is a constant. We stress that this result covers only the case of p ≥ 2, and existence of the solution to the BVP (8), (5) when p ∈ (1, 2) remains an open issue. Nevertheless, our computations in Section 5 show that this differential problem is solvable for p = 1.5 at least. Eventually, we search for a strictly monotone solution of the ODE

′

r 1−N r N −1 |ρ ′ (r )|p−2 ρ ′ (r )



= fp (ρ(r )) ,

(10)

which satisfies the boundary values

ρ ′ (0) = 0,

lim ρ(r ) = ξ .

r →∞

(11)

It is worthwhile to mention that similar BVPs with the p-Laplacian have also been studied in [14,15]. However, we have to point out that those differential equations are considered only on finite intervals. Most importantly, the studied BVPs are

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

89

solvable when their right-hand side is an odd function (see the cited papers). In our case the problem (10), (11) possesses a unique solution only if ξ < 1 in the right-hand side of (9), that is when this function is not odd (see the proof in [13]). Since we deal with numerical solution of the BVP (10), (11), when p > 1, it is convenient to discuss here some specific properties of this problem, which create additional difficulties for application of numerical methods. As in the case p = 2, this problem has singularities at r = 0 and r = ∞. However, there is an additional difficulty when p ̸= 2. The secondorder (and any higher) derivatives of the solution are unbounded at the origin and we have a singularity of the second kind (see details in Section 2). The latter implies a reduction in the convergence order of some classical numerical methods (e.g. collocation methods, which are utilized in [4]). This is the main reason why we ‘‘transfer’’ the boundary condition at the origin in our approach (that is, we replace it by an equivalent condition in a neighboring point). Moreover, since another singularity exists at infinity it makes sense to split the original BVP into two auxiliary problems for two subintervals [r0 , R] and [R, r∞ ] where r0 is close to the origin, r∞ is sufficiently large and R is a certain positive number (all details of this decomposition are explained in Section 4). In this way, each auxiliary problem has only one singular point, which is treated by means of asymptotic approximation of the solution about it. The described approach has been used successfully in the case p = 2 (see [4,6]). There is a large variety of well-known computational methods and available BVP solvers for ODEs. We do not attempt to give a comprehensive overview, but only mention some important contributions. We start with the Fortran code COLSYS and its more recent version COLNEW published in [16,17], respectively. Another series of BVP solvers (the Matlab codes bvp4c, bvp5c, bvp6c) are described in [18–20]. However, the latter codes are not specifically designed for singular problems. An important contribution to a numerical analysis of singular BVPs is given in [21]. The code sbvp and the bvpsuite series presented in [22–24] are based on numerical methods that preserve their convergence properties in the presence of singularities. However, even the cited codes may experience difficulties when applied to certain particular cases of BVP (10), (11), as, for instance, observed in [13]. A possible reason for such difficulties is the problems’ ill-conditioning and/or non-smooth solutions. An advantage of the numerical technique introduced here is that it avoids solving large systems of equations, where strong ill-conditioning may arise. 2. Asymptotic behavior of the solution at the origin In this section, we explore the behavior of the solution ρ(r ) in a neighborhood of the singular point r = 0. For that, we consider an initial value problem (IVP) of the form tx′ = Ax + φ(t , x) + g (t ),

0 < t ≤ t0 ,

lim x(t ) = 0

(12)

t →0

where A is a square constant matrix of order n with eigenvalues satisfying the condition Re λ(A) ≤ 0, the vector-function φ(t , x) is continuous in the domain (t , x) ∈ [0, t0 ] × {x ∈ Rn : ∥x∥∞ ≤ a} where ∥ · ∥∞ is the conventional sup-norm in Rn , φ(t , 0) = 0, and g (t ) is a continuous n-value vector-function in the interval [0, t0 ] such that g (0) = 0. In relation to (12), the following result has been proved in [25, Theorem 5]. Theorem 1. Suppose that the right-hand side of the ODE (12) is sufficiently smooth and Re λ(A) ≤ 0, g (t ) and φ(t , x) can be represented in the form of the convergent series g (t ) =

∞ 

g (m) t m

and φ(t , x) =



fl

(m) l m

xt ,

t ≤ t0 ,

(13)

∥l∥∞ ≥1,m≥0,∥l∥∞ +m≥2

m=1 l

where xl denotes the product x11 · · · xlnn with lj , j = 1, . . . , n, being non-negative integers, for all (t , x) ∈ [0, t0 ] × {x ∈ Rn : ∥x∥∞ ≤ a}. Then the following assertions hold. 1. For any sufficiently small t0 , a solution θ (t ) to the IVP (12) exists and can be represented in the form of the convergent series

θ (t ) =

∞ 

θ (m) t m ,

t ≤ t0 ,

(14)

m=1

where the coefficients θ (m) are found by the formal substitution of all the above expansions into (12); moreover, the function θ (t ) is the unique solution to this problem in the class of analytic functions. 2. Next, the solution presented by formula (14) is unique among all sufficiently smooth vector-functions x(t ) with the property ∥x(t )∥∞ ≤ Ct ν , where C is a constant and 0 < ν ≤ 1, as t → 0. 3. In addition, the solution (14) is unique in the class of all sufficiently smooth vector-functions x(t ) when the inhomogeneous (0) term φ(t , x) on the left-hand side of Eq. (12) is given by the second formula in (13) with the coefficients fl = 0 for all l such that ∥l∥∞ ≥ 1. Theorem 1 proves the existence and uniqueness of the solution to the singular IVP (12). Unfortunately, this result is not applicable yet to the case of the BVP (10), (11) when p ̸= 2 since expansions (13) admit only integer powers, but our Eq. (10) implies non-integer ones because of the right-hand side given by formula (9). Therefore, we have to extend Theorem 1 to this case first.

90

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

Corollary 1. Suppose that the right-hand side of the ODE (12) is sufficiently smooth and Re λ(A) ≤ 0, g (t ) and φ(t , x) can be represented in the form of the convergent series g (t ) =

∞ 

g (m) t mη



and φ(t , x) =

fl

(m) l mη

xt

,

t ≤ t0 ,

∥l∥∞ ≥1,m≥0,∥l∥∞ +m≥2

m=1 l

where xl denotes the product x11 · · · xlnn with lj , j = 1, . . . , n, being non-negative integers and η is any positive real number, for all (t , x) ∈ [0, t0 ] × {x ∈ Rn : ∥x∥∞ ≤ a}. Then the assertions 1–3 of Theorem 1 hold with t m being replaced by t mη in (14). Proof. This result is easily proved by the variable substitution τ = t η . The latter transforms the IVP (12) to the following form:

˜ + 1 f˜ (τ , x) + 1 g˜ (τ ), τ x˙ = Ax η η

0 < τ ≤ τ0 ,

lim x(τ ) = 0

(15)

τ →0

where A˜ := A/η, f˜ (τ , x) := φ(τ 1/η , x) and g˜ (τ ) := g (τ 1/η ). Clearly, the functions g˜ (τ ) and f˜ (τ , x) can be expanded in series (m) of the forms (13), respectively, but with coefficients g˜ (m) and f˜l . We conclude that the IVP (15) satisfies all the conditions mentioned in the formulation of Theorem 1. We apply now the latter result to prove Corollary 1 for any sufficiently small τ0 . In other words, we have proved the existence and uniqueness of the solution to the IVP (12) and its representation in the form of the convergent series (14) under the more relaxed conditions. So, we have derived all the necessary theoretical tools to explore the singular IVP associated with our BVP (10), (11) at the origin. We begin with the following differential problem:

(p − 1)|ρ ′ (r )|p−2 ρ ′′ (r ) + ρ(0) = ρ0 ,

N −1 r

|ρ ′ (r )|p−2 ρ ′ (r ) = fp (ρ),

(16)

lim r ρ ′ (r ) = 0

(17)

r →0+

where the right-hand function fp (ρ) is defined by (9). Assume further that, in a sufficiently small neighborhood of the point r = 0, the solution is represented in the form

ρ(r ) = ρ0 + Cr k (1 + o(1)),

as r → 0+ ,

(18)

where C and k are positive constants. Therefore, simple differentiation yields

ρ ′ (r ) = Ckr k−1 (1 + o(1)),

ρ ′′ (r ) = Ck(k − 1)r k−2 (1 + o(1)),

as r → 0+ .

(19)

It also follows from the ODE (16) that

 lim

r →0+

(p − 1)|ρ (r )| ′

p−2

ρ (r ) + ′′

N −1 r

|ρ (r )| ′

p−2



ρ (r ) = fp (ρ0 ). ′

(20)

Having substituted (19) in (20) we conclude that k = p/(p − 1) > 0, k − 1 = 1/(p − 1) > 0 and the constant

k−1

C = fp (ρ0 )/N /k. Next, taking (18) and the first formula in (19) into account we see that the solution ρ(r ) of problem (16) and its first derivative are continuous in the neighborhood of r = 0. However, according to the second formula in (19), its second derivative is not continuous since k < 2 for p > 2. This confirms the fact that the point r = 0 is singular. Now we want to boost the accuracy of the asymptotic representation (18) of the solution ρ(r ) near r = 0. For that, we perform the variable substitution ρ(r ) = ρ0 + Cr k (1 + y(r )) where the constants C and k have been determined above. The latter results in the following IVP with respect to the unknown function y(r ):





1 + y(r ) +

1 k

ry (r ) ′

p−2 

fp ρ0 + Cr k (1 + y(r ))



=

1 + y(r ) +

1 N



2(p − 1) +

N −1 k



ry (r ) + ′

p−1 Nk

r y (r ) 2 ′′





fp (ρ0 )

,

(21)

and the initial conditions are y(0) = 0 and limr →0+ ry′ (r ) = 0. Further, we are looking for a particular solution to the ODE (21) in the form of the following series: y(r ) =

∞  j=1

yj r jp/(p−1)

(22)

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

91

where 0 ≤ r ≤ r0 and r0 is a sufficiently small positive number. Formal substitution of (22) in (21) allows the coefficients yj , depending on ρ0 and ξ , to be determined. We recall that ξ is utilized in the right-hand side of Eq. (16) (see formula (9)). For instance, when j = 1 we obtain y1 (ρ0 , ξ ) =

CN (ρ0 (2ρ0 p + p − ρ0 ) − ξ (ρ0 p + 1))

(23)

2 (N (p − 1) + p) ρ0 (ρ0 + 1)(ρ0 − ξ )

where C is the constant fixed above. Finally, we prove the uniqueness of the solution (22). For that, we perform the variable substitution x1 (r ) := y(r ), x2 (r ) := ry′ (r ) and arrive at the system of two equations rx′ = Ax + φ(r , x) + g (r ) satisfying the initial condition limr →0+ x(r ) = 0. Here, we have utilized the following notation:





0 −Np/(p − 1)

A=

g (r ) =



g1 (r ) g2 (r )

1 , −N − p/(p − 1)

x(r ) =



x1 ( r ) , x2 ( r )



φ(r , x) =



f 1 ( r , x) , f 2 ( r , x)





where f1 (r , x) ≡ 0 ≡ g1 (r ), f2 (r , x) := h(r ,x1 ) (1 + x1 (r ) + (p − 1)/p)2−p − 1 − (2 − p)x1 (r ) − ((2 − p)(p − 1)/p) x2 (r ) − h(r , 0), h(r , x1 ) := fp ρ0 + Cr k (1 + x1 (r )) N /fp (ρ0 ), and g (r ) = h(r , 0) − 1. Since two eigenvalues of the matrix A are λ1 = −p/(p − 1) < 0 and λ2 = −N < 0, and the functions φ(r , x) and g (r ) satisfy all the conditions mentioned in Corollary 1, we conclude that (22) is the sole solution to the IVP (21). Having returned to the original variable we arrive at the following. Corollary 2. Suppose that N > 1 and p > 1. Then, for each ρ0 , the IVP (16), (17) has a unique solution ρ(r ) in a sufficiently small neighborhood of the point r = 0, and this solution can be presented in the following form:

ρ (r , ρ0 ) = ρ0 +

p−1



fp (ρ0 )

p

N

1/(p−1)

r p/(p−1) 1 + y1 r p/(p−1) + o r p/(p−1)







(24)

where the coefficient y1 is calculated by formula (23). We remark that the solution given by Corollary 2 is holomorphic near the origin when p = 2. However, this is not true for other values of the parameter p because derivatives of the function ρ(r ) of order higher than one become discontinuous at zero. The latter follows from formula (24). 3. Asymptotic approximation of the solution at infinity The next step is to derive an asymptotic expansion of the solution to the BVP (10), (11) as r → ∞. We apply the same approach as that published in [5,6] for the particular case p = 2. First of all we transform (10) to an asymptotically autonomous form; i.e., we represent it as a differential equation where all the coefficients tend to nonzero constants, as r → ∞. For that, we rewrite the ODE (10) in the following form:

(p − 1)ρ ′′ (r ) +

fp (ρ(r )) N −1 ′ ρ (r ) = ′ p−2 . r ρ (r )

(25)

Note that only monotonically increasing solutions ρ(r ) are considered in this paper. That is why the term |ρ ′ (r )|p−2 has been replaced with ρ ′ (r )p−2 in (25). Then, we make the variable substitution ρ(r ) = ξ + r a z (r ) where the constant parameter a := (1 − N )/(2(p − 1)). After that, Eq. (25) takes the form r a−2



r 2 (−1 + p)z ′′ (r ) −

= 2pλ2

(N − 1)(1 + N − 2p)z (r ) 4(p − 1)



(−r a z (r ))p−2 r a z (r ) (1 + ξ + r a z (r ))p−1 (ξ + r a z (r )) .  p−2 r a z ′ (r ) + ar a−1 z (r )

(26)

The boundary condition of the ODE (26) at infinity is lim z (r ) = lim z ′ (r ) = 0.

r →∞

r →∞

(27)

Having simplified the differential equation (26) with the boundary condition (27), we arrive at

(p − 1)z ′′ (r ) = 2pλ2

(−z (r ))p−2 z (r ) (ξ + a1 (r )) (ξ + 1 + a1 (r ))p−1 + a3 (r )z (r ) z ′ (r )p−2 (1 + a2 (r ))p−2

(28)

92

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

where a1 (r ) := r a z (r ), a2 (r ) := az (r )/ rz ′ (r ) , a3 (z ) := (N − 1)(1 + N − 2p)z (r )/ 4(p − 1)r 2 . Note that all the functions a1 (r ), a2 (r ), a3 (r ) satisfy the condition limr →∞ a1 (r ) = limr →∞ a2 (r ) = limr →∞ a3 (r ) = 0 and, hence, the resulting Eq. (28) is asymptotically autonomous. In order to obtain an asymptotic approximation to the solution of the differential problem (28), we consider



′′ (p − 1)z∞ (r ) = 2pλ2





z∞ (r )p−1 ξ p−1 (ξ + 1) ′ (r )p−2 z∞



(29)

with respect to an unknown function z∞ (r ). Further, we are looking for the solution to (29) in the form z∞ (r ) = c exp(τ r ) where c and τ are some constants. Having substituted the latter formula into the ODE (29) we derive the characteristic equation with two real zeros. They are:

τ1,2 = ±τ ,

(1 + ξ )ξ p−1 τ = 2pλ p−1 

2

1/p

.

(30)

Taking now into account the boundary condition (27) we conclude that the family of solutions to the ODE (29) with the mentioned boundary condition must have the form

    p−1 1/p 2 (1 + ξ )ξ z∞ (r ) = c exp − 2pλ r . p−1

(31)

It follows from (31) and the results published in [26] that the asymptotically autonomous Eq. (28) possesses a oneparameter family of solutions that satisfy the boundary condition (27) and can be expanded in the following asymptotic series: z (r ) =

∞ 

bk Ck (r )e−τ kr

(32)

k=1

where b is a positive real number and the functions Ck (r ) are computed by substitution of formula (32) into (28). The series (32) converges for any r ≥ r∞ where r∞ is sufficiently large, so that b e−τ r is small enough. Using the principal term of the series (32), we derive the sufficiently good asymptotic approximation to the needed solution of (28); i.e., z (r ) = b C1 (r )e−τ r (1 + o(1)), as r → ∞. Finally, having returned to the original variable ρ(r ), we conclude that Eq. (10) has a one-parameter family of solutions that satisfy the second boundary condition in (11) and are approximated asymptotically by the formula

ρ(r ) = ξ − b C1 (r )r a e−τ r (1 + o(1)),

as r → ∞,

(33)

where a is fixed above, τ is defined in (30) and C1 (r ) is the solution of the ODE

(p − 1)C1′′ (r ) − 2τ C1′ (r ) =

(N + 1)(N + 1 − 2p) C1 ( r ) 4(p − 1)r 2

(34)

satisfying the following conditions at infinity: limr →∞ C1 (r ) = 1 and limr →∞ C1′ (r ) = 0. As shown in [6], there exists a variable substitution that reduces Eq. (34) to a modified Bessel equation. Eventually, the asymptotic expansion of the function C1 (r ) is also known. We remark that formula (33) for p = 2 gives us the same asymptotic approximation as that published earlier in [5,6]. 4. Numerical approximation The complex numerical technique recommended for accurate solution of the BVP (10), (11) consists of a triple shooting method implemented within the bisection scheme and an efficient ODE solver with automatic global error control, as described below. As already mentioned in Section 1, we replace the discussed BVP having two singularities (at the origin and at infinity) with two auxiliary BVPs such that each of them has only one singularity. Then, we fix some constants r0 , R, and r∞ satisfying the condition r∞ > R > r0 > 0. Here, the value of r0 is close to zero and means the interval [0, r0 ] where the asymptotical representation (24) of the solution to the IVP (16), (17) works with sufficiently high accuracy (see Corollary 2). The parameter R denotes the point where the solution of the BVP (10), (11) with the right-hand side defined by (9) equals zero. We recall that the solution ρ(r ) grows monotonically in the interval −1 < ρ(r ) < ξ for r > 0 (see [13] for further explanation). So, the so-called bubble radius R is unique. The last number r∞ is large enough for an accurate asymptotic approximation of the solution given by (33) at infinity. Next, we divide the region [r0 , r∞ ] into two subintervals [r0 , R] and [R, r∞ ] and formulate two corresponding BVPs. Note that the exact value of R is unknown and has also to be computed in the course of integration of the BVP (10), (11). Further, we conduct the following scheme of computation.

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

93

First, we denote the unique monotone solution to (16) in the interval [r0 , R] by ρ− (r ). In other words, ρ− (r ) means the non-positive part of the solution ρ(r ) to the BVP (10), (11). Corollary 2 says that, with the accuracy up to higher-order terms, the boundary conditions are well approximated by

ρ− (r0 , ρ0 ) = ρ0 +

p−1



p

fp (ρ0 )

1/(p−1)

N

p/(p−1)

r0



p/(p−1)

1 + y1 r0



,

ρ− (R) = 0

(35)

where the coefficient y1 is determined in (23). We stress that the value of the parameter ρ0 in (35) is unknown at first and has to be found (together with the corresponding solution ρ− (r )) by the standard shooting method. We recall that the shooting method implies that we solve an IVP corresponding to the BVP (16), (35). For that, we compute the initial value of ′ the derivative ρ− (r0 ) by simple differentiation of the left-hand boundary condition with respect to the parameter r0 . Here, we also take into account that the parameter ρ0 in the boundary condition (35) must be in the range −1 < ρ0 < 0. The latter information is useful to determine properly an interval [ρ1 , ρ2 ] for the shooting method. Second, the function ρ+ (r ) denotes the unique monotone solution of (16) in the interval [R, r∞ ] (the non-negative part of the solution to the BVP (10), (11)). Clearly, with the accuracy up to higher-order terms, the solution ρ+ (r ) satisfies the following boundary conditions:

ρ+ (R) = 0,

a ρ+ (r∞ , b) = ξ − b r∞ C1 (r∞ )e−τ r∞ .

(36)

Notice that the second formula in (36) follows from the asymptotic behavior of the solution ρ(r ) at infinity, presented by (33) in Section 3, when all higher-order terms are neglected. Again, we apply the shooting method to compute the parameter b and the corresponding solution ρ+ (r ) satisfying (36). The value of the first derivative of the solution ρ+ (r ) at the point r∞ is found by differentiation of the second formula in the boundary conditions (36). We emphasize that the corresponding IVP is solved in the mentioned shooting method in the backward direction; i.e., starting from r∞ . When looking for an initial interval [b1 , b2 ] we should remember that the parameter b is positive. Third, the required function ρ(r ) in the region [r0 , r∞ ] is merely derived by

ρ(r ) =



ρ− (r ) ρ+ (r )

when r0 ≤ r ≤ R, when R ≤ r ≤ r∞ .

(37)

Unfortunately, the function ρ(r ) determined by formula (37) is not yet the solution to the ODE (16) in the interval [r0 , r∞ ] because, in general, its first derivative will be discontinuous. Therefore, the condition lim ρ ′ (r ) = lim ρ ′ (r )

r →R−

r →R+

(38)

must be provided additionally by varying the point R. For that, we calculate the difference

1(R) := lim ρ ′ (r ) − lim ρ ′ (r ) r →R−

r →R+

and apply the shooting method again to find the positive real parameter R that ensures condition (38); i.e., with 1(R) = 0. Eventually, we compute the bubble radius R and the solution ρ(r ) to the BVP (10), (11). However, we determine at first two values R1 and R2 such that 1(R1 ) · 1(R2 ) < 0. The described scheme of numerical integration can be presented in the form of the following pseudo-code. Algorithm 4.1. Computational Technique with Automatic Error Control Step 0. Initially, we set the accuracy parameter Tol, and find R1 , R2 , ρ1 , ρ2 , b1 , b2 ; ′ Step 1. Compute the initial vector (r0 , ρ− (r0 , ρ1 ), ρ− (r0 , ρ1 ))T ; ′ Step 2. Compute the numerical solution (r , ρ− (r , ρ1 ), ρ− (r , ρ1 ))T in [r0 , R1 ]; Step 3. Error := |ρ− (R1 , ρ1 )|; Step 4. If Error ≤ Tol, then ρ0 := ρ1 , go to Step 15; ′ (r0 , ρ2 ))T ; Step 5. Compute the initial vector (r0 , ρ− (r0 , ρ2 ), ρ− ′ Step 6. Compute the numerical solution (r , ρ− (r , ρ2 ), ρ− (r , ρ2 ))T in [r0 , R1 ]; Step 7. Error := |ρ− (R1 , ρ2 )|; Step 8. If Error ≤ Tol, then ρ0 := ρ2 , go to Step 15; Step 9. While Error > Tol do, Step 10. ρ0 := ρ1 + (ρ2 − ρ1 )/2; ′ Step 11. Compute the initial vector (r0 , ρ− (r0 , ρ0 ), ρ− (r0 , ρ0 ))T ; ′ Step 12. Compute the numerical solution (r , ρ− (r , ρ0 ), ρ− (r , ρ0 ))T in [r0 , R1 ]; Step 13. Error := |ρ− (R1 , ρ0 )|; Step 14. If ρ− (R1 , ρ0 ) · ρ− (R1 , ρ1 ) < 0, then ρ2 := ρ0 ;

94

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

else ρ1 := ρ0 ; end while; ′ Step 15. Compute the initial vector (r∞ , ρ+ (r∞ , b1 ), ρ+ (r∞ , b1 ))T ; ′ Step 16. Compute the numerical solution (r , ρ+ (r , b1 ), ρ+ (r , b1 ))T in [R1 , r∞ ]; Step 17. Error := |ρ+ (R1 , b1 )|; Step 18. If Error ≤ Tol, then b := b1 , go to Step 29; ′ Step 19. Compute the initial vector (r∞ , ρ+ (r∞ , b2 ), ρ+ (r∞ , b2 ))T ; ′ Step 20. Compute the numerical solution (r , ρ+ (r , b2 ), ρ+ (r , b2 ))T in [R1 , r∞ ]; Step 21. Error := |ρ+ (R1 , b2 )|; Step 22. If Error ≤ Tol, then b := b2 , go to Step 29; Step 23. While Error > Tol do, Step 24. b := b1 + (b2 − b1 )/2; ′ Step 25. Compute the initial vector (r∞ , ρ+ (r∞ , b), ρ+ (r∞ , b))T ; ′ Step 26. Compute the numerical solution (r , ρ+ (r , b), ρ+ (r , b))T in [R1 , r∞ ]; Step 27. Error := |ρ+ (R1 , b)|; Step 28. If ρ+ (R1 , b) · ρ+ (R1 , b1 ) < 0, then b2 := b; else b1 := b; end while; ′ Step 29. 1(R1 ) := ρ− (R1 , ρ0 ) − ρ+′ (R1 , b); Step 30. If |1(R1 )| ≤ Tol, then R := R1 , go to Step 94; ′ Step 31. Compute the initial vector (r0 , ρ− (r0 , ρ1 ), ρ− (r0 , ρ1 ))T ; ′ Step 32. Compute the numerical solution (r , ρ− (r , ρ1 ), ρ− (r , ρ1 ))T in [r0 , R2 ]; Step 33. Error := |ρ− (R2 , ρ1 )|; Step 34. If Error ≤ Tol, then ρ0 := ρ1 , go to Step 45; ′ Step 35. Compute the initial vector (r0 , ρ− (r0 , ρ2 ), ρ− (r0 , ρ2 ))T ; ′ Step 36. Compute the numerical solution (r , ρ− (r , ρ2 ), ρ− (r , ρ2 ))T in [r0 , R2 ]; Step 37. Error := |ρ− (R2 , ρ2 )|; Step 38. If Error ≤ Tol, then ρ0 := ρ2 , go to Step 45; Step 39. While Error > Tol do, Step 40. ρ0 := ρ1 + (ρ2 − ρ1 )/2; ′ Step 41. Compute the initial vector (r0 , ρ− (r0 , ρ0 ), ρ− (r0 , ρ0 ))T ; ′ Step 42. Compute the numerical solution (r , ρ− (r , ρ0 ), ρ− (r , ρ0 ))T in [r0 , R2 ]; Step 43. Error := |ρ− (R2 , ρ0 )|; Step 44. If ρ− (R2 , ρ0 ) · ρ− (R2 , ρ1 ) < 0, then ρ2 := ρ0 ; else ρ1 := ρ0 ; end while; ′ Step 45. Compute the initial vector (r∞ , ρ+ (r∞ , b1 ), ρ+ (r∞ , b1 ))T ; ′ Step 46. Compute the numerical solution (r , ρ+ (r , b1 ), ρ+ (r , b1 ))T in [R2 , r∞ ]; Step 47. Error := |ρ+ (R2 , b1 )|; Step 48. If Error ≤ Tol, then b := b1 , go to Step 59; ′ Step 49. Compute the initial vector (r∞ , ρ+ (r∞ , b2 ), ρ+ (r∞ , b2 ))T ; ′ (r , b2 ))T in [R2 , r∞ ]; Step 50. Compute the numerical solution (r , ρ+ (r , b2 ), ρ+ Step 51. Error := |ρ+ (R2 , b2 )|; Step 52. If Error ≤ Tol, then b := b2 , go to Step 59; Step 53. While Error > Tol do, Step 54. b := b1 + (b2 − b1 )/2; ′ Step 55. Compute the initial vector (r∞ , ρ+ (r∞ , b), ρ+ (r∞ , b))T ; ′ Step 56. Compute the numerical solution (r , ρ+ (r , b), ρ+ (r , b))T in [R2 , r∞ ]; Step 57. Error := |ρ+ (R2 , b)|; Step 58. If ρ+ (R2 , b) · ρ+ (R2 , b1 ) < 0, then b2 := b;

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

else b1 := b; end while; ′ Step 59. 1(R2 ) := ρ− (R2 , ρ0 ) − ρ+′ (R2 , b); Step 60. If |1(R2 )| ≤ Tol, then R := R2 , go to Step 94; Step 61. 1(R) := 1(R1 ); Step 62. While |1(R)| > Tol do, Step 63. R := R1 + (R2 − R1 )/2; ′ Step 64. Compute the initial vector (r0 , ρ− (r0 , ρ1 ), ρ− (r0 , ρ1 ))T ; ′ Step 65. Compute the numerical solution (r , ρ− (r , ρ1 ), ρ− (r , ρ1 ))T in [r0 , R]; Step 66. Error := |ρ− (R, ρ1 )|; Step 67. If Error ≤ Tol, then ρ0 := ρ1 , go to Step 78; ′ Step 68. Compute the initial vector (r0 , ρ− (r0 , ρ2 ), ρ− (r0 , ρ2 ))T ; ′ Step 69. Compute the numerical solution (r , ρ− (r , ρ2 ), ρ− (r , ρ2 ))T in [r0 , R]; Step 70. Error := |ρ− (R, ρ2 )|; Step 71. If Error ≤ Tol, then ρ0 := ρ2 , go to Step 78;

Step 72. While Error > Tol do, Step 73. ρ0 := ρ1 + (ρ2 − ρ1 )/2; ′ Step 74. Compute the initial vector (r0 , ρ− (r0 , ρ0 ), ρ− (r0 , ρ0 ))T ; ′ Step 75. Compute the numerical solution (r , ρ− (r , ρ0 ), ρ− (r , ρ0 ))T in [r0 , R]; Step 76. Error := |ρ− (R, ρ0 )|; Step 77. If ρ− (R, ρ0 ) · ρ− (R, ρ1 ) < 0, then ρ2 := ρ0 ; else ρ1 := ρ0 ; end while; ′ Step 78. Compute the initial vector (r∞ , ρ+ (r∞ , b1 ), ρ+ (r∞ , b1 ))T ; ′ Step 79. Compute the numerical solution (r , ρ+ (r , b1 ), ρ+ (r , b1 ))T in [R, r∞ ]; Step 80. Error := |ρ+ (R, b1 )|; Step 81. If Error ≤ Tol, then b := b1 , go to Step 92; ′ Step 82. Compute the initial vector (r∞ , ρ+ (r∞ , b2 ), ρ+ (r∞ , b2 ))T ; ′ Step 83. Compute the numerical solution (r , ρ+ (r , b2 ), ρ+ (r , b2 ))T in [R, r∞ ]; Step 84. Error := |ρ+ (R, b2 )|; Step 85. If Error ≤ Tol, then b := b2 , go to Step 92;

Step 86. While Error > Tol do, Step 87. b := b1 + (b2 − b1 )/2; ′ Step 88. Compute the initial vector (r∞ , ρ+ (r∞ , b), ρ+ (r∞ , b))T ; ′ Step 89. Compute the numerical solution (r , ρ+ (r , b), ρ+ (r , b))T in [R, r∞ ]; Step 90. Error := |ρ+ (R, b)|; Step 91. If ρ+ (R, b) · ρ+ (R, b1 ) < 0, then b2 := b; else b1 := b; end while; ′ Step 92. 1(R) := ρ− (R, ρ0 ) − ρ+′ (R, b); Step 93. If 1(R) · 1(R1 ) < 0, then R2 := R, 1(R2 ) := 1(R); else R1 := R, 1(R1 ) := 1(R); end while; ′ Step 94. Concatenate (r , ρ− (r , ρ), ρ− (r , ρ))T and (r , ρ+ (r , b), ρ+′ (r , b))T by formula (37); Step 95. Stop.

95

96

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

The presented numerical technique is designed for computing the solution ρ(r ) to the BVP (10), (11) with error not exceeding the user-supplied tolerance parameter Tol. However, the central role in Algorithm 4.1 is played by an ODE solver implied in Steps 2, 6, 12, 16, 20, 26, 32, 36, 46, 50, 56, 65, 69, 75, 79, 83, 89. Moreover, the numerical integration of the corresponding ODE is conducted in the half of these steps in the backward direction. So, the accuracy and efficiency of Algorithm 4.1 depends strongly on the ODE solver involved in the computation. In this paper, we recommend the Gauss-type nested implicit Runge–Kutta method of order 4 for implementation in the above algorithm because of its symmetry and many other useful properties (see [27] for more details). The adaptivity (i.e., an automatic stepsize selection) is organized here on the basis of effective local error estimation grounded in the embedded Runge–Kutta pair approach. In other words, our ODE solver implements the NIRK4(2) pair from [28]. The accuracy control of this solver includes the cheap global error estimation, which merely implies summation of the local error estimates that are available in the embedded Runge–Kutta pair under consideration almost for free, and the automatic local–global error control mechanism published in [28, Section 3]. The latter error control procedure is intended for producing numerical solutions with the global error not exceeding the tolerance parameter Tol (it is denoted by ϵg in the cited paper). This also follows from the theoretical results in [28] that the global error control theory works in the backward direction as well. The last fact to mention here is that NIRK4(2) is an implicit Runge–Kutta pair with excellent stability properties. Thus, we need additionally an iteration to treat arising nonlinear problems. Here, we utilize the simplified Newton iteration, as explained in [29, Section 6]. The cited paper proves that performing three iterates with trivial predictor in every step of the ODE solver is sufficient for the proper working of all the error estimators and error control algorithms mentioned here. Finally, we point out that Algorithm 4.1 computes the numerical solution only in the interval [r0 , r∞ ] ⊂ [0, ∞). So, its extension to the regions [0, r0 ) and (r∞ , ∞) is done on the basis of the asymptotic expansions derived in Sections 2 and 3. In the next section, we present and discuss numerical results obtained by our computational technique. We emphasize that Algorithm 4.1 has been coded in MATLAB and run on a personal computer with processor Intel Pentium IV, 3.0 GHz under operating system MICROSOFT WINDOWS XP. This is an ordinary home PC. Therefore, its power and the implemented arithmetic may be insufficient to treat accurately the BVP (10), (11) for some values of the parameters ξ and p. However, the mentioned problem is easily solved by using a better computational device. 5. Numerical results Our purpose in this section is to confirm that the derived numerical technique is well-suited to solving the considered singular BVP for various values of the parameters. We also observe and discuss properties of its solution depending on ξ and p. Special attention is paid to some data that have a known physical meaning. For instance, ρ0 is the gas density in the center of the bubble, and R is the bubble radius, as we have already mentioned in Section 4. We also calculate approximate values of the energy integral Jp (ρ) defined in (6), which can be written in the form Jp (ρ) =

∞

 0



ρ ′ (r )p p

 ρ(r )

 + Wp (ρ(r )) r N −1 dr

where Wp (ρ(r )) = ξ fp (u)du. The energy integral Jp (ρ) is computed by the composite trapezoidal rule on the mesh derived by Algorithm 4.1 for the solution ρ(r ). The inner integral Wp (ρ(r )) is evaluated by the composite Simpson’s rule on an equidistant mesh and with automatic error control for Tol/10. We recall that the function fp is fixed by (9) with λ = 1.0. So, the accuracy of the integral calculation is sufficiently high. Since the BVP (10), (11) has been solved many times for p = 2.0 (see, for example, [4–6]) we focus our attention on other values of p in the range from p = 1.5 up to p = 4.0. Concerning the parameter ξ , due to the condition 0 < ξ < 1, we consider ξ ∈ [0.1, 0.9]. The dimension N equals 3 in all our experiments. Table 1 displays numerical data collected from application of Algorithm 4.1 to the BVP (10), (11) for various parameters ξ and p. Here, we exhibit the computed values of ρ0 , b, R and Jp . We also show the length of the computation interval given by the parameter r∞ . Notice that the initial point r0 = 1.0e−3 in all the experiments of this paper. The last column in Table 1 means the accuracy of computation. More precisely, the data have been obtained for the numerical solution ρ(r ) with the magnitude of the error 1(R) in the first derivative (see Section 4 for the precise definition of this error) not exceeding the displayed number. Sometimes, the accuracy attained by our method is sufficiently low (see the data corresponding to p = 1.5 and ξ = 0.8 or to p = 2.0 and ξ = 0.9) or even absent (as for p = 1.5 and ξ = 0.9) because of the computer’s limitations mentioned in the concluding part of Section 4. It is also worthwhile to mention that the data computed for the particular case p = 2.0 are in agreement with those published in [4–6]. Some properties of the solutions are better seen in pictures. The plots of the density profiles ρ(r ) are displayed for ξ = 0.1, 0.2, . . . , 0.9 in Figs. 1 and 2. Each plot contains 6 graphics corresponding to p = 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 except for the case p = 1.5 in Fig. 2(c) (because of the computer’s limitations mentioned above). Fig. 3(a) and Fig. 4 exhibit the energy integrals and the bubble radiuses, respectively. Note that the integrals Jp (ρ) are plotted in the logarithmic format because they change quickly with ξ . For instance, for ξ ’s close to 1, this growth is nearly exponential. This follows from the almost linear behavior of the plotted integrals in Fig. 3(a).

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

97

Table 1 Numerical data for N = 3 and various values of p and ξ . p

ξ

ρ0

b

1.5

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

−0.6105636596680 −0.8604096984863 −0.9586869468689 −0.9913790631294 −0.9990432856750 −0.9999676963931 −0.9999999084199 −0.9999999999996

1.1424e+4 9.9341e+4 9.4891e+5 1.3672e+7 4.2310e+8 5.2715e+10 5.9129e+15 7.9214e+23

3.630270 3.105436 3.163314 3.475155 4.027547 4.929282 6.488611 9.662500

2.0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

−0.3046752929688 −0.5677530575562 −0.7707031250000 −0.9031249237061 −0.9711191558838 −0.9953001050949 −0.9997789021683 −0.9999995735622 −1.0+1.0e−15

3.9936 10.390 28.251 96.535 4.9390e+2 5.1169e+3 2.1375e+5 2.8020e+8 2.6484e+17

2.5

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

−0.2047297668457 −0.4093780517578 −0.6041381835938 −0.7737518310547 −0.8997742080689 −0.9712115764618 −0.9964619178772 −0.9999488453198 −0.9999999998790

3.0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

3.5

4.0

R

Jp

r∞

Accuracy

0.671493 1.261771 2.128140 3.536647 6.027463 10.91927 22.16883 56.29626

8 8 8 8 8 8 12 16

1.0e−6 1.2e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−4

3.321875 2.685732 2.582332 2.720984 3.070007 3.695892 4.816895 7.130994 14.12100

0.063317 0.213432 0.484458 0.969585 1.883850 3.774039 8.316614 22.63499 1.075669e+2

8 8 8 8 8 8 12 16 20

1.0e−5 4.7e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−3

0.2672 0.5298 1.0558 2.3855 6.9466 32.219 1.1360e+2 3.8275e+3 3.6884e+9

1.619730 1.430513 1.477381 1.677657 2.043221 2.619959 3.559356 5.397265 10.85258

0.006664 0.038554 0.118447 0.289453 0.648076 1.444708 3.466394 10.12723 51.10367

8 8 8 8 8 8 12 16 16

7.8e−7 4.5e−7 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6

−0.1552688598633 −0.3195709228516 −0.4898437500000 −0.6587505340576 −0.8114341735840 −0.9256558227539 −0.9848483676910 −0.9993912789345 −0.9999999651611

0.0575 0.1045 0.1860 0.3556 0.8054 2.5509 2.6559 14.689 4.3576e+5

2.112975 1.821796 1.763281 1.824774 2.006012 2.359473 3.025146 4.433972 8.765149

8.399469e−4 0.007943 0.032455 0.096503 0.249502 0.620513 1.623580 5.095253 27.33090

8 8 8 8 8 8 12 16 16

1.0e−7 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

−0.1253654479981 −0.2621994018555 −0.4103939056397 −0.5677215576172 −0.7264457702637 −0.8678085327148 −0.9627944183350 −0.9972033920288 −0.9999989060682

0.0202 0.0357 0.0605 0.1064 0.2111 0.5391 0.2697 0.5053 1.8825e+3

1.731190 1.534294 1.501674 1.557507 1.705295 1.991333 2.535901 3.699914 7.296338

1.224783e−4 0.001831 0.009797 0.035216 0.105003 0.291592 0.832711 2.808298 16.01524

8 8 8 8 8 8 12 16 16

2.3e−7 1.0e−7 1.0e−7 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

−0.1052383422852 −0.2223907470703 −0.3527099609375 −0.4967224121094 −0.6519195556641 −0.8066757202148 −0.9317675857544 −0.9921192970276 −0.9999888246345

0.0094 0.0164 0.0271 0.0455 0.0833 0.1859 0.0569 0.0517 48.671

1.452829 1.316882 1.302441 1.356099 1.482837 1.722894 2.181201 3.168555 6.233691

1.996156e−5 4.609515e−4 0.003193 0.013792 0.047318 0.146709 0.457504 1.658530 10.05654

8 8 8 8 8 8 12 16 16

1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6 1.0e−6

First of all we remark on an interesting property of all the numerical solutions to the BVP (10), (11) computed in this section; i.e., having fixed a certain p we observe that the bubble radius R increases both as ξ → 1 and as ξ → 0 (see Fig. 4). The latter implies that there must exist a minimum radius Rmin corresponding to some value of ξ . This critical value is referred to as ξcr . Such a behavior, which agrees with experimental evidence, has been observed for p = 2.0 in [30]. The present computations show that the same happens for the other values of p discussed in Section 5. The existence of the minimum radius Rmin for all the p’s is confirmed by Fig. 4. So, the value of ξcr has been computed with high accuracy for every p and collected together with other experimental data in Table 2. The behavior of the density profiles for the critical values of ξ is presented in Fig. 2(d). Comparing now the profiles corresponding to a certain value of p, but for different ξ ’s, we observe the same behavior as that noted in [5,6] for the case p = 2.0; i.e., for values of ξ greater than ξcr , the graphics feature an interior layer where the

98

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

(a) Density profiles for ξ = 0.1 and various p.

(b) Density profiles for ξ = 0.2 and various p.

(c) Density profiles for ξ = 0.3 and various p.

(d) Density profiles for ξ = 0.4 and various p.

(e) Density profiles for ξ = 0.5 and various p.

(f) Density profiles for ξ = 0.6 and various p.

Fig. 1. Density profiles for ξ = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 and various p.

solution ρ(r ) grows rapidly (see Fig. 1(c–f) and Fig. 2(a–c)). The latter corresponds to the interface between the phases and is referred to as the bubble interface (or the bubble wall). The thickness of this layer decreases monotonically as ξ → 1 (see Fig. 3(b)). For the remaining values of ξ , i.e., when ξ < ξcr , the solution ρ(r ) changes smoothly for r > 0 (see Fig. 1(a, b)).

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

(a) Density profiles for ξ = 0.7 and various p.

(b) Density profiles for ξ = 0.8 and various p.

(c) Density profiles for ξ = 0.9 and various p.

(d) Density profiles for ξα and various p.

99

Fig. 2. Density profiles for ξ = 0.7, 0.8, 0.9 and ξcr and for various p.

(a) Energy integrals for various ξ and p.

(b) Interface thickness for various ξ and p.

Fig. 3. The energy integrals (the left-hand plot) and the interface thickness (the right-hand plot) for various ξ and p.

At the same time, the derivative ρ ′ (r ) does not take large values in the same region. We also conclude that, for a fixed value of ξ , the bubble radius R grows in line with the decreasing p. Notice that the minimum radius Rmin declines with increasing p (see Table 2).

100

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

(a) Bubble radiuses for various ξ and p = 1.5.

(b) Bubble radiuses for various ξ and p = 2.0.

(c) Bubble radiuses for various ξ and p = 2.5.

(d) Bubble radiuses for various ξ and p = 3.0.

(e) Bubble radiuses for various ξ and p = 3.5.

(f) Bubble radiuses for various ξ and p = 4.0.

Fig. 4. The bubble radiuses according to different definitions for various ξ and p.

Besides the definition given above of the bubble radius R (the zero of the density profile), there exist other definitions used in the literature. For instance, in [30], the bubble radius for p = 2.0 and N = 3 is defined on the basis of a statistical

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

101

Table 2 Numerical values of ξcr and other data for N = 3. p

ξcr

ρ0

b

Rmin

Jp

r∞

Accuracy

1.5 2.0 2.5 3.0 3.5 4.0

0.2294810 0.2860938 0.2984375 0.2943750 0.2828500 0.2678000

−0.8999939002991 −0.7466207504273 −0.6012420654297 −0.4802055358886 −0.3842236638069 −0.3092498600483

1.8828e+5 24.332 1.0438 0.1798 0.0552 0.0230

3.088550 2.580018 2.113538 1.763097 1.500401 1.298910

1.479301 0.436897 0.116653 0.030311 0.007625 0.001836

8 8 8 8 8 8

2.9e−7 1.0e−7 1.0e−7 1.0e−7 1.0e−8 1.0e−8

Table 3 The interface thickness, the bubble radiuses and the surface tension according to different definitions and various values of ξ and p = 1.5, 2.0, 2.5. p

ξ

1.5

0.1 0.2

δ

σS

σIGS

3.630270 3.105436 3.088550 3.163314 3.475155 4.027547 4.929282 6.488611 9.662500

2.971784 2.724380 2.752794 2.909198 3.299513 3.905672 4.846145 6.434437 9.630673

2.397633 2.357848 2.425835 2.660124 3.131772 3.793251 4.771281 6.386326 9.602404

1.167813 0.710909 0.645369 0.534220 0.434618 0.369424 0.323617 0.289873 0.263993

0.308963 0.624362 0.699257 0.856220 1.050835 1.237412 1.427501 1.624573 1.829792

0.350865 0.680900 0.754159 0.902249 1.081781 1.256722 1.438986 1.630758 1.832553

3.321875 2.685732 2.580018 2.582332 2.720984 3.070007 3.695892 4.816895 7.130994 14.12100

2.962618 2.437234 2.381143 2.390235 2.572710 2.959495 3.617423 4.764681 7.099966 14.10745

2.065371 1.814443 1.892827 1.920844 2.228938 2.722594 3.461297 4.665711 7.042688 14.09898

1.156500 0.849679 0.726851 0.712104 0.628136 0.567122 0.519147 0.480583 0.449232 0.423306

0.035270 0.160112 0.314523 0.341208 0.534759 0.726040 0.922654 1.133718 1.363489 1.614176

0.044599 0.194488 0.365830 0.393907 0.585479 0.762432 0.945043 1.146131 1.369098 1.611749

2.642776 2.204561 2.113538 2.113556 2.192877 2.432404 2.889832 3.733765 5.499102 10.89785

2.505513 2.087669 2.006180 2.006349 2.096845 2.351733 2.827751 3.690545 5.472818 10.88594

1.619730 1.430513 1.475449 1.477381 1.677657 2.043221 2.619959 3.559356 5.397265 10.85258

1.000362 0.811850 0.722343 0.721267 0.664917 0.622936 0.587317 0.555851 0.528405 0.504641

0.005656 0.043533 0.129900 0.131666 0.265145 0.425816 0.604218 0.805154 1.035561 1.299719

0.007614 0.056519 0.160757 0.162802 0.308526 0.465711 0.631412 0.820836 1.042954 1.301723

R

ξcr

0.3 0.4 0.5 0.2 0.2 0.8 2.0

0.1 0.2

ξcr

0.3 0.4 0.5 0.2 0.2 0.8 0.9 2.5

0.1 0.2

ξcr

0.3 0.4 0.5 0.2 0.2 0.8 0.9

RIGR

RIGS

approach. Having extended that result to an arbitrary p and N we would arrive at the formula RIGR (ξ ) :=

 −1 

∞



ρ (r , ξ ) r ′

p

N −1

∞

r ρ ′ (r , ξ )p r N −1 dr .

dr

0

(39)

0

Another definition of the bubble radius accommodated to the case of the N-dimensional p-Laplacian comes from [31] and has the following form: RIGS (ξ ) :=



(N − 1)

1/N

∞



ρ (r , ξ ) r ′

p

N −1

dr



−1/N −Wp (ρ0 (ξ ))

(40)

0

where Wp (ρ0 ) is defined at the beginning of this section. In Fig. 4 and Tables 3 and 4, values of the bubble radiuses according to (39) and (40) are compared with our computed bubble radius R for various ξ and p. The parameter ξcr in these tables means the critical value of ξ for a certain p and its numerical approximation can be found in Table 2 for all the p’s considered here. Note that the data corresponding to p = 2.0 in Table 3 agree well with those computed and published earlier in [5,6]. For all the p’s, the existence of a minimal bubble radius is verified, and the value of ξcr does not depend significantly on the utilized definition of the bubble radius. In most cases, we observe that R(ξ ) > RIGR (ξ ) > RIGS (ξ ) for the same ξ (with some exceptions for p = 3.5 and p = 4.0). It is also worthwhile to remark that the values of the bubble radius according to all these definitions almost coincide as ξ → 1, i.e., when the thickness of the bubble interface becomes tiny. All of this is confirmed by Fig. 3(b), Fig. 4 and Tables 3 and 4.

102

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104 Table 4 The interface thickness, the bubble radiuses and the surface tension according to different definitions and various values of ξ and p = 3.0, 3.5, 4.0. p

ξ

R

RIGR

RIGS

δ

σS

σIGS

3.0

0.1 0.2

2.112975 1.821796 1.763097 1.763281 1.824774 2.006012 2.359473 3.025146 4.433972 8.765149

2.089558 1.783673 1.715038 1.714751 1.770608 1.952311 2.312934 2.990364 4.412047 8.755051

1.286295 1.158770 1.185978 1.190685 1.333428 1.610510 2.073581 2.838530 4.324888 8.716636

0.901737 0.774918 0.714420 0.711668 0.672140 0.643269 0.618244 0.594232 0.571531 0.550825

0.001084 0.013066 0.049664 0.052912 0.133206 0.253430 0.405028 0.587486 0.808833 1.076777

0.001523 0.017746 0.064649 0.068676 0.162825 0.288581 0.432942 0.604518 0.817217 1.079163

1.731190 1.534294 1.500401 1.501674 1.557507 1.705295 1.991333 2.535901 3.699914 7.296338

1.765841 1.542002 1.492352 1.490820 1.533112 1.672847 1.957850 2.508189 3.681539 7.287689

1.049708 0.964217 0.981635 0.992147 1.102556 1.320180 1.699910 2.342466 3.586582 7.245900

0.842129 0.751248 0.710941 0.704527 0.675256 0.654553 0.636981 0.619062 0.600541 0.582620

2.301472e−4 0.004208 0.017499 0.022191 0.068402 0.153068 0.276540 0.438023 0.646068 0.912457

3.334583e−4 0.005908 0.023739 0.029859 0.086907 0.180741 0.302723 0.455269 0.654943 0.915107

1.452829 1.316882 1.298910 1.302441 1.356099 1.482837 1.722894 2.181201 3.168555 6.233691

1.518114 1.351819 1.318725 1.315903 1.352519 1.466588 1.700104 2.159360 3.153053 6.226186

0.878712 0.820931 0.832300 0.847647 0.938111 1.115071 1.432223 1.983926 3.052579 6.182026

0.805764 0.737420 0.711110 0.701450 0.678887 0.663451 0.650972 0.637788 0.622770 0.607219

5.226709e−5 0.001425 0.005665 0.009636 0.035881 0.093696 0.191403 0.331930 0.524913 0.786584

7.755724e−5 0.002052 0.007951 0.013333 0.047016 0.114167 0.214565 0.348711 0.533957 0.789420

ξcr

0.3 0.4 0.5 0.2 0.2 0.8 0.9 3.5

0.1 0.2

ξcr

0.3 0.4 0.5 0.2 0.2 0.8 0.9 4.0

0.1 0.2

ξcr

0.3 0.4 0.5 0.2 0.2 0.8 0.9

Mathematically, the interface thickness is defined as follows:

δ(ξ ) :=

∞



(r − R(ξ ))2 ρ ′ (r , ξ )p r N −1 dr

1/p 

0

∞

ρ ′ (r , ξ )p r N −1 dr

−1/p

.

(41)

0

Formula (41) is an extension of the definition presented in [30] for p = 2.0 to the case of arbitrary p and N. Next, the surface tension is an important property of the bubbles, which can also be measured through the density profile ρ(r ) and its derivative. Having extended the formula given in [31] to arbitrary p and N, the surface tension of a bubble for the case of the p-Laplacian is now defined as follows: ∞

  σIGS (ξ ) := (N − 1)

ρ ′ (r , ξ )p r N −1 dr

1/N



(N −1)/N −Wp (ρ0 (ξ )) .

(42)

0

Another definition of the surface tension is used in cosmology and particle physics (see, for instance, [32]). Having applied that definition to our case we arrive at

σS (ξ ) :=

∞



ρ ′ (r , ξ )p dr .

(43)

0

Fig. 5 and Tables 3 and 4 exhibit values of the surface tension calculated by (42) and (43) for various p and ξ . As expected, the surface tension (according to these two definitions) increases in line with increasing ξ . Such a behavior does not depend on a specific value of p. On the other hand, we observe that the sign of the curvature of the graph changes (compare Fig. 5(a) with Fig. 5(c–f)). It also follows from Fig. 5 and Tables 3 and 4 that σS (ξ ) < σIGS (ξ ) for the same ξ . However, when ξ = 0.9 formulas (42) and (43) produce almost the same result. 6. Conclusion An extension of the density profile equation to nonhomogeneous fluids has been analyzed in this paper. More precisely, the classical Laplacian has been replaced with a degenerate one. The right-hand side of the mentioned equation is defined in such a way that the corresponding BVP possesses exactly one strictly monotone solution, as explained in [13]. Numerical simulations suggest that some properties of the solution known for the classical case p = 2.0 are also preserved for other values of p in the range 1.5 ≤ p ≤ 4.0.

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

(a) Surface tension for various ξ and p = 1.5.

(b) Surface tension for various ξ and p = 2.0.

(c) Surface tension for various ξ and p = 2.5.

(d) Surface tension for various ξ and p = 3.0.

(e) Surface tension for various ξ and p = 3.5.

(f) Surface tension for various ξ and p = 4.0.

Fig. 5. The surface tension according to different definitions for various ξ and p.

103

104

G.Yu. Kulikov et al. / Journal of Computational and Applied Mathematics 262 (2014) 87–104

In particular, we would like to point out the power of the numerical technique presented in Section 4 (see Algorithm 4.1). Our method is effective for numerical integration of the BVP (10), (11) for all the parameters ξ and p considered in this paper. Moreover, the accuracy of integration can be as high as we want provided that the method has been implemented on a computational device with a sufficiently long arithmetic. For instance, we are able to treat successfully the singular BVP with p-Laplacian for values of p and ξ that are not covered yet in the literature on this topic. Thus, our result has good potential for solving many important problems not only in fluid mechanics, but also in other areas of application where similar differential equations may arise. Finally, we emphasize that the existence and uniqueness result presented in [13] is applicable to the BVP (10), (11) only if p ≥ 2.0. On the other hand, the computations performed in Section 5 confirm that Algorithm 4.1 works for p = 1.5 as well. Thus, we plan to study theoretically the existence and uniqueness of the solution ρ(r ) to the BVP under consideration for p ∈ (1.0, 2.0) in the future. Acknowledgments The first two authors are grateful for the support from Portuguese National Funds through the Fundação para a Ciência e a Tecnologia (FCT) within the scope of the project PEst-OE/MAT/UI0822/2011. The authors are also grateful to the anonymous referees for their valuable remarks and comments on the paper. References [1] P. Seppecher, Moving contact lines in the Cahn–Hilliard theory, Internat. J. Engrg. Sci. 34 (1996) 977–992. [2] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy, J. Chem. Phys. 28 (1958) 258–267. [3] M.E. Gurtin, D. Polignone, J. Vinals, Two-phase binary fluids and immiscible fluids described by an order parameter, Math. Models Methods Appl. Sci. 6 (1996) 815–831. [4] G. Kitzhofer, O. Koch, P.M. Lima, E. Weinmüller, Efficient numerical solution of the density profile equation in hydrodynamics, J. Sci. Comput. 32 (2007) 411–424. [5] N.B. Konyukhova, P.M. Lima, M.L. Morgado, M.B. Soloviev, Bubbles and droplets in nonlinear physics models: analysis and numerical simulation of singular nonlinear boundary value problems, Comp. Maths. Math. Phys. 48 (2008) 2018–2058. [6] P.M. Lima, N.B. Konyukhova, N.V. Chemetov, A.I. Sukov, Analytical-numerical investigation of bubble-type solutions of nonlinear singular problems, J. Comput. Appl. Math. 189 (2006) 260–273. [7] I. Rachunkova, J. Tomecek, Strictly increasing solutions of a nonlinear singular differential equation arising in hydrodynamics, Nonlinear Anal. 72 (2010) 2114–2118. [8] I. Rachunkova, J. Tomecek, Bubble-type solutions of nonlinear singular problems, Math. Comput. Modelling 51 (2010) 658–669. [9] I. Rachunkova, L. Rachunek, J. Tomecek, Existence of oscillatory solutions of singular nonlinear differential equations, Abstr. Appl. Anal. 2011 (2011) Article Id. 408525. [10] B. Sciunzi, E. Valdinoci, Mean curvature properties for p-Laplace phase transitions, J. Eur. Math. Soc. 7 (2005) 319–359. [11] A. Petrosyan, E. Valdinoci, Density estimates for a degenerate/singular phase-transition model, SIAM J. Math. Anal. 36 (2005) 1057–1079. [12] V.O. Savin, E. Valdinoci, B. Sciunzi, Flat level set of p-Laplace phase transitions, Mem. Amer. Math. Soc. 182 (858) (2006). [13] G. Hastermann, P.M. Lima, M.L. Morgado, E.B. Weinmüller, Density profile equation with p-Laplacian: analysis and numerical simulation, Appl. Math. Comput. (to appear). [14] P. Drabek, R.F. Manasevich, P. Takac, Manifolds of critical points in a quasi-linear model for phase transitions, Contemp. Math. 540 (2011) 95–134. [15] P. Takac, Stationary radial solutions for a quasilinear Cahn–Hilliard model in N space dimensions, Electron. J. Differential Equations, Conf. 17 (2009) 227–254. [16] U. Ascher, J. Christiansen, R.D. Russel, A collocation solver for mixed order systems of boundary value problems, Math. Comp. 33 (1978) 659–679. [17] U. Ascher, U. Bader, A new implementation for a mixed order boundary value ODE solver, SIAM J. Sci. Stat. Comput. 8 (1987) 483–500. [18] L. Shampine, J. Kierzienka, M. Reichelt, Solving boundary value problems for ordinary differential equations in MATLAB with bv p4c, 2000. Available at ftp://ftp.mathworks.com/pub/doc/papers/bvp/. [19] J. Kierzienka, L. Shampine, A BVP solver that controls residual and error, J. Numer. Anal. Indust. Appl. Math. 3 (2008) 27–41. [20] L. Shampine, P. Muir, H. Xu, A user-friendly FORTRAN BVP solver, J. Numer. Anal. Indust. Appl. Math. 1 (2006) 201–217. [21] F. de Hoog, R. Weiss, Difference methods for BVPs with a singularity of the second kind, SIAM J. Numer. Anal. 13 (1976) 775–813. [22] W. Auzinger, G. Kneisl, O. Koch, E.B. Weinmüller, Manual and Code: SBVP, A MATLAB solver for singular boundary value problems, 2003. Available at http://www.math.ist.utl.pt/ewa. [23] W. Kitzhofer, G. Pulverer, O. Koch, C. Simon, E.B. Weinmüller, BVPSUITE, A new MATLAB code for singular implicit boundary value problems, 2009. Available at http://www.math.ist.utl.pt/ewa. [24] W. Kitzhofer, O. Koch, G. Pulverer, C. Simon, E.B. Weinmüller, The new MATLAB code BVPSUITE for the solution of singular implicit boundary value problems, J. Numer. Anal. Inal. Indust. Appl. Math. 5 (2010) 113–134. [25] N.B. Konyukhova, Singular Cauchy problems for systems of ordinary differential equations, USSR Comput. Maths. Math. Phys. 23 (1983) 72–82. [26] A.M. Lyapunov, Probleme General de la Stabilite du Mouvement, in: Annals of Mathematics Studies, vol. 17, Princeton University Press, Princeton, 1947. [27] G.Yu. Kulikov, S.K. Shindin, Adaptive nested implicit Runge–Kutta formulas of Gauss type, Appl. Numer. Math. 59 (2009) 707–722. [28] G.Yu. Kulikov, Cheap global error estimation in some Runge–Kutta pairs, IMA J. Numer. Anal. 33 (2013) 136–163. [29] G.Yu. Kulikov, A.I. Merkulov, S.K. Shindin, Asymptotic error estimate for general Newton-type methods and its application to differential equations, Russian J. Numer. Anal. Math. Modelling 22 (2007) 567–590. [30] F. Dell’Isola, H. Gouin, G. Rotoli, Nucleation and shell-like interfaces by second-gradient theory: numerical simulations, Eur. J. Mech. B/Fluids 15 (1996) 545–568. [31] F. Dell’Isola, H. Gouin, P. Seppecher, Radius and surface tension of microscopic bubbles second-gradient theory, C.R. Acad. Sci. Paris 320 (1995) 211–216. (serie IIb). [32] A.P. Linde, Particle Physics and Inflationary Cosmology, Harwood Academic, Chur, Switzerland, 1990.