An integral method for the numerical solution of nonlinear singular boundary value problems Текст научной статьи по специальности «Математика»

MSC 65R20 DOI: 10.14529/mmp150401

AN INTEGRAL METHOD FOR THE NUMERICAL SOLUTION OF NONLINEAR SINGULAR BOUNDARY VALUE PROBLEMS

M. V. Bulatov, Matrosov Institute for System Dynamics and Control Theory of Siberian Branch of Russian Academy of Sciences; Irkutsk National Research Technical University, Irkutsk, Russian Federation, [email protected], P.M. Lima, Instituto Superior Tecnico, University of Lisbon, Lisbon, Portugal, [email protected],

Thanh Do Tien, Irkutsk National Research Technical University, Irkutsk, Russian Federation, [email protected]

We discuss the numerical treatment of a nonlinear singular second order boundary value problem in ordinary differential equations, posed on an unbounded domain, which represents the density profile equation for the description of the formation of microscopic bubbles in a non-homogeneous fluid. Due to the fact that the nonlinear differential equation has a singularity at the origin and the boundary value problem is posed on an unbounded domain, the proposed approaches are complex and require a considerable computational effort. This is the motivation for our present study, where we describe an alternative approach, based on the reduction of the original problem to an integro-differential equation. In this way, we obtain a Volterra integro-differential equation with a singular kernel. The numerical integration of such equations is not straightforward, due to the singularity. However, in this paper we show that this equation may be accurately solved by simple product integration methods, such as the implicit Euler method and a second order method, based on the trapezoidal rule. We illustrate the proposed methods with some numerical examples.

Keywords: density profile equation; singular boundary value problem; integro-differential equation; implicit Euler method.

1. Introduction

1.1. Density Profile Equation

The singular boundary value problem we discuss here originates from the Cahn-Hillard theory, which is used in hydrodynamics to study the behavior of non-homogeneous fluids. In this theory, an additional term involving the gradient of density (grad p) is added to the classical expression E0(p) for the volume free energy, depending оn the density p of the medium. Hence, the total volume free energy of a nonhomogeneous fluid can be written as

E (p,grad(p)) = Eo(p) + 2(grad(p))2, (1)

where E0(p) is a double-well potential, whose wells define the phases. The potential E0(p)

p

In [2], the density profile equation for the description of the formation of microscopical bubbles in a non-homogeneous fluid (in particular, vapor inside one liquid) is derived.

Let us briefly recall how this equation is obtained. The state of a non-homogeneous fluid (see [2] and [3]) is described by the following system of partial differential equations:

pt + div(pv) = 0, (2)

dv + V(p(p) - Y△p) = 0, (3)

where p, v denote the density and the velocity of the fluid, p represents its chemical

Y

potential and stationary, system (2), (3) is reduced to a single equation of the form

Y^p = p(p) — po, (4)

where p0 is a constant, depending on the state of the fluid. When searching for a solution of (4) with spherical symmetry which depends only on the variable r, we introduce as usual the system of spherical coordinates in lRn and equation (4) is then reduced to the following ordinary differential equation (ODE):

Y (V + ~r--p^ = p(p) - Po, r E (0, to). (5)

Since we consider the case of a spherical bubble, ODE (5) is closed with the boundary conditions

p'(0) = 0 (6) (following from spherical symmetry) and

lim p(r) = pi > 0, (7)

r—y^G

where pl is the density of the liquid surrounding the bubble. In the simplest models for

p

the difference p — p0 has 3 real roots. Taking into account that p(pl) = p0, the right-hand side of (5) may be written in the form

p(p) — po = 4a(p — pi)(p — p2)(p — pi), 0 < pi <p2 < pi, a > 0. (8)

Finally, in order to diminish the number of parameters in the equation we introduce the new variable

p — P2 x =-,

P2 — Pi

define the positive constant A = (p2 — pi) , and denote £ = ^Z^ > 0. Then, without loss of generality, instead of (5)-(7) we can investigate the boundary value problem

n_1

x''(r) +--x'(r) = 4A2(x(r) + 1)x(r)(x(r) — £), (9)

x'(0) = 0, x(to) = (10)

A

chosen as A = 1 without restriction of generality, n is the dimension of the problem, which

in the physically meaningful case is equal to 3, and which is varied in the range [0,1] such as to reflect different physical situations.

Note that problem (9), (10) always has the constant solution x(r) = £, which physically corresponds to the case of a homogeneous fluid (without bubbles).

We are interested in computing a monotonously increasing solution for 0 < r < to, the so called "bubble-type solution". When such a solution exists it has exactly one zero R in that interval, where R is interpreted as the bubble radius. Furthermore, it can be shown that —1 < x(0) < 0 and —1 < x(r) < £, r > 0. The derivative of the solution attains a maximum at some value f < R, and tends to 0 at infinity. Finally, it turns out that the solution features an interior layer, which becomes sharper for £ ^ 1. All these properties have been discussed in [4] (see also [5] and [6]).

1.2. Existence and Uniqueness of Solution

It is worth to remark that the existence of a strictly increasing solution to the problem (9), (10) is far from being a simple question. In [4], it was shown (using a variational approach developed in [7]), that such a solution can exist only if £ satisfies 0 < £ < 1. Furthermore, based on the results of [8], it is possible to show that this restriction on £ is also a sufficient condition for the existence of such a solution. These results agree with the experimental evidence and the numerical simulations reported, for example, in [2].

It is worth to remark that the density profile equation can be extended to a more general context, where the free energy of the mixture of fluids is given by

c

E (p, grad(p)) = Eo(p) + - |grad(p)|p, (11)

p

where p > 1 and c> 0. In the case p = 2, expression (11) reduces to (1).

In the general case, the differential operator on the left-hand side of (9) has the form p

r

l-n(rn-l\x\r)\p-2x\r))' = fp(x), r> 0, (12)

where fp is a function which has the same roots and the same sign as the right-hand side

p

look for a strictly increasing solution of (12), which satisfies the boundary conditions (10). The existence and uniqueness of solution of this problem was discussed in [9], where the numerical solution of this boundary value problem by collocation methods was described.

A different numerical approach to the solution of problem (9), (10), (12) was introduced in [10], where the singular boundary value problem is reduced to a sequence of auxiliary intitial value problems, which are solved by means of computational methods with global error control.

2. Integral Formulation

Equation (9) may be written in the form

r1-n(rn-1x' (r))' = f (x(r)), (13)

where f (x) = 4A2(x — £)(x + 1)x, 0 < £ < 1, A is a positive number (tipically A = 1). According to the integral method we rewrite equation (13) in the form

Note that (14) is a Volterra integro-differential equation of the first kind with a singular kernel. The numerical integration of such equations is not straightforward, due to the singularity.

In the theory of singular integral and integro-differential equations, product integration methods have been often used to solve problems of this type. These methods are recommended when the considered integrand function is the product of two parts, one of which is singular. This approach was first introduced by Weiss in [11,12]. A detailed description of the methods is given in the monographs [13-15]. Their history can be found in the survey paper [16].

Taking into account the form of the integral on the right-hand side of equation (14), we have decided to apply product integration methods to its numerical solution.

x

x' (0) = 0 and

The first of these conditions is satisfied by any solution of (14). In order to satisfy the second boundary condition we need to know that equations (14) and (13) have only 3 kinds of solutions:

1. If x(0) < x*, then the solution x(r) blows up at a finite r;

2. If x(0) > x*, then the solution x(r) is oscillatory and limr—TO x(r) = 0;

3. If x(0) = x*, then the solution x(r) is monotonic and limr—TO x(r) = £.

Obviously, what we need is to find a solution of the third type. Moreover, we know that the value x* is determined uniquely for each £ and satisfies x* E [— 1, 0]; we have x* ^ — 1 as £ ^ 1, and x* ^ 0, as £ ^ 0 (see [4]).

3. Numerical Methods

In order to solve the boundary value problem (14), (15) we have implemented a first order and a second order method.

3.1. First Order Method

In the first case, we use the implicit Euler method to approximate equation (14). First we introduce a uniform mesh on the interval [0,T]: ri = ih,i = 1,..., N, with stepsize h, such that Nh = T. Then we approximate the solution x % a vector xh = (xo,x\,... ,xN), such that xi ~ x(ri).

The components of this vector must satisfy the equation:

r > 0.

(14)

lim x(r) = £.

(15)

r—y^G

(16)

This results from approximating x'(r) by (x(r) — x(r — h))/h and using the right rectangles rule to approximate integral on the left-hand side of (14).

At each step, for a given i we solve a nonlinear equation for xi+i. This is done by the fixed point method, using the intial approximation x°+1 = xi.

x*

According to this method, we start with a certain interval [a, b] C [—1, 0], such that a) if x(0) = a , then the solution x(r) (approximated by the implicit Euler method) is of the first

x(0) = b x(r)

have x* E [a, b]. Then, as usual, we construct a sequences of intervals [ak, bk]k = 1, 2,... such that [ak,bk] C [ak_1,bk_1] , bk — ak = (bk-1 — ak-1)/^d x* E [ak,bk]■ The iteration process stops when bk — ak < e, for a given e.

3.2. Second Order Method

In this case we approximate equation (14) by the second order scheme

3xi+2 — 4xi+1 — xi = 2(2 £ (rj)n-1f (xj) + rr+_21f(xi+2^ , i = 0, ...,N — 2. (17)

x1

h2

x1 = xo + — f (x°), (18)

x

For each value of i (starting with i = 0) we determine xi+2 by solving the nonlinear equation (17) by the fixed point method.

4. Numerical Results

In this section we present the results of some numerical experiments we have carried out to test the performance of the proposed methods. All the programs were implemented in MATLAB.

0<£<1

£

values of x(0) obtained by the implicit Euler and the second order method, with h = 0, 001. For comparison, we give also the values presented in [10]. Note that the numerical scheme used in this work envolves an ODE solver with variable step size, with control of the global error at each step, so that it provides accuracy of about 8 digits. We see that the results obtained by the second order method are in good agreement with the ones presented in [10] (they have in general 3 common digits); in the case of the implicit Euler method, about 2

£1

difficult to approximate the solution. The convergence of the implicit Euler method was tested in various examples. Some results are displayed in Table 2, where we consider the approximation of x(1), in the case £ = 0, 5, with different stepsizes. With the purpose of checking the convergence order, we compute the following coefficient:

K = l ( \xh - xhn\ \

2 V\Xh/2 - Xh/4 \J

Table 1

x(0)

£

t x(0) (1st order) x(0) ( 2nd order) Result in [10]

0,1 -0,2999 -0,3046 -0,3047

0,2 -0,5597 -0,5679 -0,5672

0,3 -0,7636 -0,7708 -0,7707

0,4 -0,8990 -0,9031 -0,9031

0,5 -0,9696 -0,9712 -0,9711

0,6 -0,9950 -0,9953 -0,9953

0,7 -0,9998 -0,9998 -0,9998

0,8 -0,99992178 -0,99992178 -0,9999995

These results indicate that the implicit Euler method has the first order of convergence, when applied to this problem, as it could be expected. The same conclusion follows from the results displayed in Table 3, where the implicit Euler method is applied to the computation of x(0), for the same value of

Table 2

Convergence order of the approximations of x(1) by the implicit Euler method in the case £ = 0, 5

h x(l) K

1/10 -0,6252 0,9615

1/20 -0,6215 1,078

1/40 -0,6196 1,169

1/80 -0,6187 -

1/160 -0,6182 -

Finally, the same value was approximated using the method described in Sec. 3.2. Once again, the numerical results displayed in Table 4 confirm that this method has the second order of convergence, as expected.

Table 3

x(0)

by the implicit Euler method in the case £ = 0, 5

h x(0) K

1/100 -0,96958965 1,00002

1/200 -0,97036934 1,00000

1/400 -0,97075918

1/800 -0,97095410

Table 4

Convergence order of the approximations of x(0) by the method of Sec. 3.2 in the case £ = 0, 5

h х(0) К

1/100 -0,97113484 1,88

1/200 -0,97112361 1,66

1/400 -0,97112057

1/800 -0,97111961

5. Conclusions and Future Work

In this work we have analysed new algorithms for the numerical solution of a nonlinear singular boundary value problem, arising in the mathematical modelling of mixtures of fluids. The integral formulation provides an alternative approach to analyze the problem and obtain numerical approximations. The advantage of this approach is that the resulting integro-differential equation can be efficiently solved by simple methods, with low computational complexity. So far, we have applied only the first and the second order methods, and therefore the accuracy of the results is not so high as in the case of the more sophisticated methods used in other works. In the future we intend to apply the integral approach with discretization methods of higher order. Another possible direction of future research is the extension of this approach to the case of the density profile equation with degenerate Lapplacian (see (11)).

Research was supported by RFBR grant No 15-01-03228. References

1. Gurtin M.E., Polignone D., Vinals J. Two-Phase Binary Fluids and Immiscible Fluids Described by an Order Parameter. Mathematical Models and Methods in Applied Sciences, 1996, vol. 6, pp. 815-831. DOI: 10.1142/S0218202596000341

2. Dell'Isola F., Gouin H., Rotoli G. Nucleation of Spherical Shell-Like Interfaces by Second Gradient Theory: Numerical Simulations. European Journal of Mechanics - B/Fluids, 1996, vol. 15, pp. 545-568.

3. Gavrilyuk S.L., Shugrin S.M. Media with Equations of State that Depend on Derivatives. Journal of Applied Mechanics and Technical Physics, 1996, vol. 37, pp. 177-189. DOI: 10.1007/BF02382423

4. Lima P.M., Chemetov N.V., Konyukhova N.B., Sukov A.I. Analytical-Numerical Investigation of Bubble-Type Solutions of Nonlinear Singular Problems. Journal of Computational and Applied Mathematics, 2006, vol. 189, pp. 260-273. DOI: 10.1016/j.cam.2005.05.004

5. Kitzhofer G., Koch O., Lima P.M., Weinmuller E. Efficient Numerical Solution of the Density Profile Equation in Hydrodynamics. Journal of Scientific Computing, 2007, vol. 32, pp. 411-424. DOI: 10.1007/sl0915-007-9141-0

6. Konyukhova N.B., Lima P.M., Morgado M.L., Soloviev M.B. Bubbles and Droplets in Nonlinear Physics Models: Analysis and Numerical Simulation of Singular Nonlinear Boundary Value Problems. Computational Mathematics and Mathematical Physics, 2008,

vol. 48, no. 11, pp. 2018-2058. DOI: 10.1134/S0965542508110109 [Пузырьки и капельки в моделях нелинейной физики: анализ и численное решение сингулярной нелинейной краевой задачи / Н.Б. Конюхова, П.М. Лима, М.Л. Моргало. М.Б. Соловьев // Журнал вычислительной математики и математической физики. - 2008. - Т. 48. - №11. -С. 2019-2023.J

7. Derrick G. Comments on Nonlinear Wave Equations as Models for Elementary Particles. Journal of Mathematical Physics, 1965, vol. 5, pp. 1252-1254. DOI: 10.1063/1.1704233

8. Gazzola F., Serrin J., Tang M. Existence of Ground States and Free Boundary Problems for Quasilinear Elliptic Operators. Advances in Differential Equations, 2000, vol. 5, pp. 1-30.

9. Hastermann G., Lima P.M., Morgado M.L., Weinmuller E.B. Density Profile Equation with p-Laplacian: Analysis and Numerical Simulation. Applied Mathematics and Computation, 2013, vol. 225, pp. 550-561. DOI: 10.1016/j.amc.2013.09.066

10. Kulikov G.Yu., Lima P.M., Morgado M.L. Analysis and Numerical Approximation of Singular Boundary Value Problems with the p-Laplacian in Fluid Mechanis. Journal of Computational and Applied Mathematics, 2014, vol. 262, pp. 87-104. DOI: 10.1016/j.cam.2013.09.071

11. Weiss R., Anderssen R.S. A Product Integration Method for a Class of Singular First Kind Volterra Equations. Numerische Mathematik, 1972, vol. 18, pp. 442-456. DOI: 10.1007/BF01406681

12. Weiss R. Product Integration for the Generalized Abel Equations. Mathematics of Computation, 1972, vol. 26, pp. 177-186. DOI: 10.1090/S0025-5718-1972-0299001-7

13. Brunner H. Collocation Methods for Volterra Integral and Related Functional Equations. Cambridge, University Press, 2004. DOI: 10.1017/CB09780511543234

14. Brunner H., Van Der Houven P.J. The Numerical Solution of Volterra Equations. Amsterdam, North-Holland, CWTI Monographs 3, 1986.

15. Linz P. Analytical and Numerical Methods for Volterra Equations. SIAM, Philadelphia, 1985. DOI: 10.1137/1.9781611970852

16. Brunner H. 1896-1996: One Hundred Years of Volterra Integral Equations of the First Kind. Applied Numerical Mathematics, 1997, vol. 24, pp. 83-93. DOI: 10.1016/S0168-9274(97)00013-5

Received May 5, 2015

УДК 519.62 DOI: 10.14529/mmpl50401

ИНТЕГРАЛЬНЫЙ МЕТОД ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ НЕЛИНЕЙНЫХ СИНГУЛЯРНЫХ КРАЕВЫХ ЗАДАЧ

М.В. Булатов, П.М. Лима, До Тиен Тхань

В статье предложены численные методы решения нелинейной краевой задачи для обыкновенного дифференциального уравнения второго порядка, заданного на полуоси и неразрешенного относительно главной части. Такие задачи описывают плотность микроскопических пузырьков в неоднородной жидкости. В связи с тем, что исходное нелинейное дифференциальное уравнение неразрешено относительно главной части, и

краевая задача рассматривается на полуоси, то ранее разработанные подходы являются сложными и требуют значительных вычислительных затрат. Именно этот факт послужил мотивацией для данной статьи, где мы описываем альтернативный подход, в котором предложено записать исходную задачу в виде интегро-дифференциального уравнения типа Вольтерра с особенностью в ядре. Итак, исходную задачу мы записали в виде интегро-дифференциального уравнения типа Вольтерра с сингулярным ядром и, в виду специфики исходной задачи, условием на правом конце. Численное интегрирование таких уравнений также достаточно сложная задача. В данной работе мы предлагаем специальные методы решения таких уравнений первого и второго порядков. Приведены численные расчеты модельных примеров по предлагаемым алгоритмам. Данные расчеты показали перспективность дальнейшего развития такого о одхода.

Ключевые слова: уравнение плотности; сингулярная краевая задача; интегро-дифференциальное уравнение; неявный метод Эйлера.

Михаил Валерьянович Булатов, доктор физико-математических наук, профессор, Институт динамики систем и теории управления имени В.М. Матросова СО РАН; Иркутский национальный исследовательский технический университет (г. Иркутск, Российская Федерация), [email protected].

Педро Мигель Лима, Университет Лиссабона (г. Лиссабон, Португалия), [email protected].

До Тиен Тхань, аспирант, Иркутский национальный исследовательский технический университет (г. Иркутск, Российская Федерация), [email protected].

Поступила в редакцию 5 мая 2015 г.