ON THE CONSTRUCTION OF SOLUTIONS TO A PROBLEM WITH A FREE BOUNDARY FOR THE NON-LINEAR HEAT EQUATION Текст научной статьи по специальности «Математика»

DOI: 10.17516/1997-1397-2020-13-6-694-707 УДК 517.958:519.633

On the Construction of Solutions to a Problem with a Free Boundary for the Non-linear Heat Equation

Alexander L. Kazakov*

Matrosov Institute for System Dynamics and Control Theory SB RAS

Irkutsk, Russian Federation

Lev F. Spevak"

Institute of Engineering Science Ural Branch RAS Ekaterinburg, Russian Federation

Ming-Gong Lee*

Chung Hua University Hsinchu City, Taiwan

Received 08.06.2020, received in revised form 14.07.2020, accepted 10.08.2020 Abstract. The construction of solutions to the problem with a free boundary for the non-linear heat equation which have the heat wave type is considered in the paper. The feature of such solutions is that the degeneration occurs on the front of the heat wave which separates the domain of positive values of the unknown function and the cold (zero) background. A numerical algorithm based on the boundary element method is proposed. Since it is difficult to prove the convergence of the algorithm due to the non-linearity of the problem and the presence of degeneracy the comparison with exact solutions is used to verify numerical results. The construction of exact solutions is reduced to integrating the Cauchy problem for ODE. A qualitative analysis of the exact solutions is carried out. Several computational experiments were performed to verify the proposed method.

Keywords: non-linear heat equation, heat wave, boundary element method, approximate solution, exact solution, existence theorem.

Citation: A.L. Kazakov, L.F. Spevak, M.-G. Lee, On the Construction of Solutions to a Problem with a Free Boundary for the Nonlinear Heat Equation, J. Sib. Fed. Univ. Math. Phys., 2020, 13(6), 694-707. DOI: 10.17516/1997-1397-2020-13-6-694-707.

Introduction

We consider the non-linear parabolic heat equation [1] with a source (sink)

Tt = A^(T ) + $(T), (1)

which is also called the generalized porous medium equation [2]. If ^(0) = 0 and $(T) is power function Eq. (1) can be written as

ut = uAu + ■y(Vu)2 + auß. (2)

* [email protected] https://orcid.org/0000-0002-3047-1650 [email protected] https://orcid.org/0000-0003-2957-6962 [email protected] https://orcid.org/0000-0001-9405-2247 © Siberian Federal University. All rights reserved

Here, y > 0, a, ft = 0 are constants, a > 0 means the presence of a source and a < 0 corresponds to a sink. In what follows Cartesian coordinates are used.

Various forms of equation (2) are used to describe processes in continuum mechanics [2,3], plasma physics [1], etc. This mathematical object has a distinctive property related to the propagation of perturbations with finite velocity, which is not typical for parabolic equations.

For equation (2) the problem of initiating a heat wave is considered. The heat wave is a construction consisting of two hypersurfaces: u(t, x) > 0 and u(t, x) = 0 that continuously joined along some sufficiently smooth manifold r(t, x) = 0. The latter determines the front of the heat wave. Since the front is unknown in advance and it is determined simultaneously with the construction of the unknown function, we have a special problem with a free boundary [3], where u|r(tx)=0 = 0. The boundary conditions have the form

u|6(x)=0 = f (t, x), f (0, x) = 0, (3)

where b(x), f (t, x) are sufficiently smooth functions.

Previously, problem (2), (3) was already considered in the case a = 0, i.e., without a source (sink). Solutions were constructed both in the form of special series [4] and with the use of the boundary element method (BEM) [5]. In this paper, we propose an approximate method for constructing solutions to the problem of heat wave initiation based on the BEM.

There are various approaches to solve boundary value problems for parabolic equations using the BEM. The most natural one is to use the BEM based on time-dependent fundamental solutions [6]. This method is not suitable for solving problem (2), (3) because of the non-linearity of the right-hand side and because of the presence of a movable boundary (the heat wave front). Therefore, it is preferable to use a time-stepping BEM [7], where the boundary value problem for the elliptic equation is considered at each step, and the fundamental solution of the Laplace equation is used. In addition to the classical BEM, it is used in the method of fundamental solutions (MFS) [8], as well as in the interior field method (IFM) [9] and the null field method (NFM) [10] in annular domains for the case of circular symmetry. For problem (2), (3), we obtain the Poisson equation with a non-linear right-hand side at each time step. The dual reciprocity boundary element method (DRBEM) [11] is most suitable for solving this equation.

We did not establish the convergence of the developed method. Therefore, to verify numerical results exact solutions in the form of travelling wave are used [2]. Construction of exact solutions is reduced to the solution of the Cauchy problem for a second-order ordinary differential equation with a singularity.

Finding exact solutions of non-linear partial differential equations is an important field of modern mathematics. There is a wide variety of methods to find exact solutions. Among these methods we emphasize the group analysis method that was proposed and developed by L. V. Ovsiannikov and his colleagues [12,13]. A review of methods for constructing exact solutions to equations of mathematical physics can be found, for example, in handbook [14]. Various generalizations and modifications of the method of separation of variables [15,16] are especially often used to construct exact solutions of non-linear parabolic equations having the form (1).

1. Formulation of the problem

In the case of one spatial variable, equation (2) can be written as

2 Vu p g

ut = uupp + Yup + -p- + a.up. (4)

Here, v = 0,1,2 corresponds to Cartesian, cylindrical and spherical coordinate systems, respectively; p = J = x2, where xi are the Cartesian coordinates. Condition (3) has the form

u\P=R = f (t), f (0) = 0, (5)

where R > 0 is some constant which must obviously be positive for v = 0.

The cases v = 0 and v = 1 are considered here. It follows from previous results [17] that for Cartesian and cylindrical coordinate systems problem (4), (5) has unique analytical heatwave type solution (in the form of a convergent Taylor series). The coefficients of the series are determined from the solution of systems of linear algebraic equations. However, the radius of convergence of the series is usually small, and it can be estimated only in some spacial cases. To solve this problem and obtain an approximate solution to the problem of heat wave initiating at a given time interval [0, t*], we usually use a step-by-step method based on the boundary element approach [5]. We also note that presence of additional term (source) requires significant modification of the previously developed approach.

2. Solution algorithm based on BEM

Problem (4), (5) is solved in the specified time interval t e [0, t*] by a step-by-step method based on the BEM. At each time step tk = kh, where h is the step size, we solve the spatial problem obtained from (4), (5) with t = tk. The solution domain where the unknown function u is positive is the interval p e [0, a(tk)), and p = a(t) is the equation of motion of the heat wave front, u(tk, a(tk)) = 0.

For certainty, it is assumed that the heat wave front moves from the origin. Since a(tk) is unknown in advance, the solution domain is also unknown at the moment t = tk. That is why, we interchange the desired function and the spatial variable p [5]. Equation (4) takes the following form

2

ptpu = upUu - ipu---+ auppu. (6)

p

We rewrite equation (6) in the form of the Poisson equation and obtain at t = tk the boundary value problem

puu = F(u,p,pt,pU), p\u=L = R, (7)

where F(u,p,pt,pu) = (ptpu + Ypu)/u + vpu/p + auQ-1 pu, p = p(tk,u) is the unknown function and L = f (tk). The unknown heat wave front for the original problem is defined by the condition

p\u=0 = a(tk).

At the front of the heat wave we have [18]

q(p)\ 0 = dp

q lu=0 dn

= -7TTT, (8)

u=0 a'(tk)

where q(p) is the flow of p(tk ,u), n is the external normal to the boundary of the solution domain, n(0) = -1, n(L) = 1. It follows from the results presented in [17] that a'(tk) = 0.

Thus, we arrive to the boundary value problem (7), (8) in the domain u e [0, L]. Using the boundary element method, we write the solution of this problem in the following form

f L

p(v) = qlpu*(v, 0) + q2ku* (v, L) - plkq*(v, 0) - Rq*(v,L) - F (u, p, pt, pu) u*(v,u)du. (9)

0

Here, v e (0, L), plk = p(ik, 0), p2k = p(ik,L) = R, g^ = q(p)(tk, 0), q^ = q(p)(tk,L), u*(v,u) is the fundamental solution of the one-dimensional stationary problem, q*(v,u) = du*(v,u)/dn [6]. Values of p1 k, qfk and q2Pk are not specified by the boundary conditions and should be determined.

Taking the limits v ^ 0 and v ^ L in equation (9), we obtain the system of two boundary integral equations

Plk - q[kL = Qi, Plk + q{kL = Q2, (10)

L L

where Qi = R — J F (u,p,pt,pu) u*(0,u)du, Q2 = R + / F (u,p,pt,pu) u*(L,u)du. 0 0

Using quadratic approximation of the function p(t, 0) = a(t) on the interval [tk-i,tk], one

can write condition (8) in the following approximate form

pik — pi(k-l)

Yh(q[Pk + qft - i)) 2q(p)q(p)

2qik qi(k-i)

(11)

Expressing pik from (11) and substituting it into the first equation (10), we obtain an equation q( )

q k

in the unknown q( )

L

(q( k)2 + ^Q i — p i (k-1) —

Yh

2q(p) 2q i(k-i)

q(k + Yjh = 0.

(12)

Since the derivative of the unknown function in (4), (5) is negative along the front up |p=a(t) < 0 (for chosen direction of the movement of the heat wave front), the inverse function obeys the inequality pu|u=o = 1/(up|u=o) < 0. Therefore, q[pk> = n(0)pu|u=o = —pulu=o > 0, and a suitable solution of equation (12) has the form

q ( k)

1

2L

n 1 Yh

p i (k-i) — Q i +

2qi(k- i) + \

Q i — p i (k-i) —

Yh

2q(p) 2q i (k-i)

+ 2YhL

(13)

Substituting (13) into (10), we can find pik and

p = q(p) T + Q q(P) = Q2 — pik pi k = qi kL + Q u q2k =—l—'

(14)

Since function F(u, p, pt, pu) in (7) depends on the unknown function and its derivatives, we use the following iterative procedure to solve problem (7), (8). Let us take p(0) = R and Qi = 0, Q2 =0 as the initial values. Then the i-th iteration of the solution has the form

p(i) (v) = q[PkKi)u* (v, 0) + q(2Pkmu* (v, L) —

—pfkq*(v, 0) — Rq*(v,L) — I F(u,p(i-i),prL',p(:-i))u*(v,u)du.

~(i-i) p(i-i) Jii-i)),

(15)

The boundary values of the unknown function and flow can be found according to (13), (14)

as

q(p)(i) = q1k =

2L

(i-1) Yh

pi(k-i)— Qi +:rw

2qi(k-i) + \

(i-1) Qi — pi(k-i)

Yh

2q(p) 2q1(k-1)

+ 2YhL

2

L

0

2

1

pfk = + Qii-1), = (Q2i-1) - p^/L.

Here, q-, g(/-1) are calculated from the previous iteration:

fL fL

q(- 1 ) = R — F(i-()u*(0,u)du, Q2-1 ) = R +/ F(i-1 )u*(L,u)du, (16) Jo Jo

where F(i-1) = F (^u, p(i-1), p[i-1), pU-1)^J. The iteration process is terminated at the n-th

iteration if

(Jn) P(n-i))/P(n) (pik - Pik )/Pih

< e, where e is a given constant. Then the approximate solution

of (7), (8) at t = tk is p(tk, u) = p(n)(u). Since this solution is continuous, the solution of problem (4), (5) at t = tk, u(tk,p) can be determined without a loss of accuracy.

The developed algorithm allows us to construct a solution that is continuous with respect to a spatial variable for problem with a free boundary (4), (5) at each given time step.

L

To calculate integrals f F (u, p, pt, pu) u*(v,u)du in (15), (16), we use the dual reciprocity

o

method [11] based on the expansion of F (u, p, pt, pu) in terms of radial basis functions (RBF)

n

F (u,p,pt,pu) = a.k(u). (17)

k = 1

Functions depend on the distance between the current point and collocation points u1,u2 ,...,un that belong to the interval [0, L]: $k (x) = $ (rk), where rk = \u — uk |. For each function $k there is a function wk such that $k = Awk. After substituting expansion (17) into the integrand and twice integrating by parts, we obtain the following equality

,L

/ F (u,p,pt,pu) u*(v,u)du =

n o (18)

= 53 ak[—Wk(v)+ Pk(0)u*(v, 0)+ pk(L)u*(v,L) — wk(0)q*(v, 0) — Wk(L)q*(w, L)],

k=1

where pk(x) = dwk(u)/dn. The coefficients ak for each iteration are determined from the system of equations obtained from (17) for the current iteration p(i)(u) at the collocation points

n

F (uj, p(i) (uj), pti) (uj), pU (uj)) =J2 ak Mu ),j = 1, 2,...,n.

k=1

The use of the simplest functions = rk as RBFs [5,18] results in stable convergence of iterative processes and good accuracy of solutions. However, it is rather complicated to use these functions to solve problems with the source term because convergence of iterative processes is unstable and depends on the parameters of the problem. Obviously, the additional non-linear term requires a more precise expansion.

It is difficult to analyse the influence of RBFs on convergence analytically. Therefore, we perform a numerical analysis of the influence of used RBFs on convergence and accuracy of the solution. We consider linear function = 1+ rk, polyharmonic spline = r| and multi-quadric function = + er'2. Stable convergence is observed when the last two functions are used. At the same time, multi-quadric function ensures better accuracy of the solution so we use it in calculations. The results are shown in Section 4.

3. Construction and study of exact solutions

It is problematic to prove convergence of the algorithm based on the BEM. One possible way is to construct and study exact solutions and then use them to verify numerical results.

Plain geometry. We consider the non-linear heat equation with a source for v = 0

ut = uuxx + yuX + aug. (19)

We assume that in this case the heat wave front is defined by following equation

u(t,x)lx=a(t) =0. (20)

If functions a(t) and ug are analytical (ft € N) then it follows from the previously proved theorems (see, for example, [16]) that problem (19), (20) has unique analytical solution in form of a power series with respect to variable z = x — a(t) with recurrently determined coefficients. We consider the case when the condition of the analyticity of the source is not satisfied. Let ft > 0, ft € R in equation (19). For non-integer ft the source function can not be expanded into a Taylor series with respect to powers of u. We construct the solution in the form of a travelling wave u = v(z), z = x — pt — n, p > 0, n > 0. Due to the invariance of equation (19) one can take n = 0.

Substituting v(z) in (19), we obtain the following ordinary differential equation

vv" + y(v')2 + pv' + avg =0. (21)

Solving this equation with the condition v(0) = 0, we obtain the solution of problem (19), (20) which is a heat wave with the front x = a(t) = pt + n (if it exists). Properties of solutions of this type for ft € N were previously studied [16]. For non-integer values of ft, as far as we know, this problem was not previously considered and it is now studied for the first time.

In this case one cannot impose arbitrary condition for the derivative at z = 0. If we set z = 0 in both parts of Eq. (21) then we obtain the quadratic equation

Y (v ' (0))2 + pv' (0) = 0,

which has two roots v'(0) = —p/Y and v'(0) = 0. For other values of v'(0) equation (21) is incompatible. Thus, we obtain from the continuity condition that either

v(0) = 0, v (0) = —p/Y, or (22)

v(0) = 0, v (0) = 0. (23)

Problems (21), (22) and (21), (23) have specific properties that they inherit from original problem (19), (20). In particular, at the starting point z = 0 function v(z) (that is the term at the higher derivative) turns to zero. It means that problems cannot be written in normal form, and the classical theorems on existence of solutions of the Cauchy problem are not applicable. It is obvious that (21), (23) has trivial solution v = 0. Special study is required to prove the existence and uniqueness of the solution of (21), (22) and the existence of a non-trivial solution of (21), (23).

In what follows we consider only the case a > 0. It means that there is an influx of energy (matter) into the system. This case is widely encountered in applications.

Theorem 1. The Cauchy problem (21), (22) for a > 0, ft ^ 1 has unique solution v(z) £

C[z0,0] n Cfz0t0) for some z0 <

Proof. Using substitution v' = p, we lower the order of equation (21). Then it takes the form

vp-j- + YP2 + hp + a.v{i = 0. (24)

dv

Using the change of variable w = in equation (24), we obtain

„ dp aw

ftw-f- + YP + H +-=°. (25)

dw p

Since for ft > 0 the equality |^=0 = 0 holds then for equation (24) conditions (22) correspond to

-(0) = -H/Y. (26)

Let us perform a qualitative analysis of equation (21). To do this, we consider the equivalent dynamic system

dw ft dp Y 2 a

— = -wp, — = —p - p--w. (27)

d£ h d£ h H

Here the parametrization is performed in such a way that dz = Hwd£. System (27) is very close to system (4.3) from [19]. The first equations differ by a positive constant multiplier on the right-hand side. The second equations are the same.

System (27) has two singular points M1(0, -h/y) and M2(0,0). It follows from qualitative analysis [19] that M1 has the topological type "saddle", and M2 is the "saddle-node" with one nodal and two saddle sectors. Considering condition (26), we are primarily interested in phase trajectories that enter point M1 and/or leave it. In this case, there is no need to consider the left phase half-plane w < 0 (except ft = k/m, where k, m are natural odd numbers) because w = .

There is a separatrix S that tends to the point M1 as £ ^ +0. Phase trajectories which are located to the right of S bypass the nodal sector bounded by S and the coordinate axis Op. The phase trajectories inside the nodal sector tend to M2 as £ ^ (Fig. 1).

The separatrix S corresponds to the solution v = v*(z) of problem (21), (22). Taking into account the conditions of the theorem, the solution has the following properties: 1) v is defined and continuous on some interval [z0,0]; 2) v is twice continuously differentiable on the interval (zo, 0); 3)@v*(0) = v*(zo) = 0, v*(z) > 0; 4) <(z) changes sign once, v^(zmax) = 0, v^(zm&x) = vmax, lim v'(z) = 5) v''(z) < 0. Let us note that condition ft > 1 guarantees that

z^zo+0

properties 1 and 2 hold. The schematic representation of function v(z) is shown in Fig.2.

It seems impossible to find the exact values z* ,zmax,vmax. Further we discuss how to find interval estimates for them.

Remark 1. For ft = k/m, where k, m are natural odd numbers, k > m, the left phase half-plane w = < 0 can also be considered. The solution of problem (21), (22) can be continued to the right from the point z = 0 as well as problem (21), (23) at z < 0 has a non-trivial solution. Both solutions are negative so they are meaningless from the physical point of view.

Interval estimates are constructed for the key parameters of the solution of (21), (23) in the particular case ft =1 to simplify the study. We follow the procedure suggested in [19]. Let us use linear substitution of the unknown function and independent variable

5 = Az,v = Bv, A = aY/H(Y + 1)], B = aY2/[H2(Y + 1)]. (28)

P 0 -Wy - -0 "max

Fig. 1. Phase portrait Fig. 2. The traveling wave configuration

Then problem (21), (23) takes the form

vv'' + y(v')2 + yv' + (y + 1)v = 0, v(0) = 0, v'(0) = — 1. (29)

Here the tilde is omitted.

We analyse the properties of solution of (29) and find estimates for z*,zmax, vmax. In the proof of Theorem 1 we show that for z £ [zmax, 0] function v decreases, —1 < v' < 0 and v'(zmax) = 0. Then from (29) we have that —1 — y < v'' < —1, and v''(0) = —1, v''(zmax) = —y — 1.

Integrating the upper bound for v'' on the interval [z, 0] and taking into account the Cauchy conditions (29), we obtain that

v' < —(y + 1)z — 1, v > —0.5(y + 1)z2 — z, zmax < z < 0. (30)

It follows from the first inequality (30) that zmax < — 1/(y +1). The right-hand side of the second inequality at zi = — 1/(y +1) ^ zmax takes the maximum value that satisfies the inequality v(z1) > 1/[2(y + 1)]. Since v(z1) < vmax we obtain that vmax > 1/[2(y + 1)]. Integrating the lower bound for v on the interval [z, 0] we obtain that

v' ^ —z — 1, v ^ —0.5z2 — z, 0 < z ^ zmax. (31)

It follows from (31) that zmax > —1 and vmax < 0.5.

When z £ [z0,zmax] function v increases, i.e., v' > 0. Taking into account this inequality, it follows from (29) that v'' < —y — 1. Integrating it twice, we obtain vmax—v > (y+1)(zmax — z)2/2, z0 ^ z ^ zmax. Since vmax ^ 0.5, for z = z0 we have 0 < zmax — z0 ^ 1 /^/y + 1 whence it follows that —1/^/y + 1 — 1 ^ z0 < — 1/(y + 1). So, we have all the required estimates. Let us apply the transformation inverse to (28). Then we obtain the following inequalities

mt+1) . . M /—it, ^ M ^ p (1 + y) iqo,

- < zmax <--;--w Y + 1+ Y +1) < z0 <zmax; ~-2 < vmax < -~-2-. (32)

aY aY aY 2aY2 2aY2

We omit cumbersome transformations that are required to construct estimates of the form (32) for the general case. Now a solution of equation (19) can be found. The solution has the form of a heat wave that propagates with the constant velocity u(t, x) = v*(x — pt — n).

An interesting feature of this solution is that it is a soliton (solitary wave). We point out that if a = 0 (no source) then solution has explicit form u = —p(x — pt — n)/Y. It is easy to show that it is not a soliton.

Derivation of travelling wave solution. Since one can not obtain analytical solution of problem (21), (22) we solve it numerically. It follows from the proof of Theorem 1 that solution should be constructed on the interval [z0,0], and z0 < 0 is unknown. We only know that v(z0) = 0 and v'(z) ^ to when z ^ z0. This is the problem with a free boundary, and its formulation is non-standard for solving by the boundary element method. Note that we prefer to use the BEM because, unlike difference methods, it allows one to construct a continuous solution.

For convenience, we use substitution V(z) = v(-z) and resolve equation (21) with respect to the higher derivative. Taking into account (22), we have the Cauchy problem

V'' = V

-V' - y (V')2 - aV, V(0) = 0, q(0) = --. (33)

J y

It is difficult to find its solution with satisfactory accuracy on the interval z £ [0, z*], where z* = —z0, because of specific properties (mentioned above) of the unknown function. A correct solution can be found in two stages. In the first stage, problem (33) is solved on the interval z £ [0,L], where L < z* is selected in such a way that V'(L) < 0, i.e., L > —zmax. Then the BEM iterative procedure is used to obtain continuous solution of this problem (see [20]).

One cannot construct a solution of problem (33) on the interval z £ [L, z*] in the original formulation because the derivative is unbounded. Then, in the second stage, we interchange the independent variable and the unknown function, just as we did in Section 2. As a result, we obtain the inverse Cauchy problem for the unknown function z(V).

1

z'' = y [y - -z' + aVP(z')2] , z(L*) = L, qz(L*) = —, (34)

where V £ [0, L*], qz = dz/dn, L* = V(L); V(L) and q(L) are obtained in the first stage.

Now problem (34) can be solved with the use of the iterative BEM. As a result, in particular, the value z(0) = z* is found. The continuousness of the found function allows us to determine V(z) for z £ [0, z*] and the solution of problem (21), (22) without loss of accuracy. Note that estimates (32) can be used to select parameter L.

Cylindrical geometry. For v =1, there are no exact travelling wave solutions for equation (4). It is known [16] that for 3 =1 equation (4) has a self-similar solution u(t,p) = = a(t)a'(t)v(r), where r = p/(Reet), and R > 0,0 are constant. Function v(r) satisfies the Cauchy problem

vv'' + y(v')2 + (r + +(R - 2) v = 0, v(1) = 0, v'(1) = - -. (35)

This solution corresponds to a heat wave with an exponential law of motion of the heat wave front a(t) = Reet. Since equation (35) is not autonomous, it is difficult to perform a qualitative analysis in this case, and we don't consider this task here.

In conclusion, we note that solution of problem (35) by the iterative boundary element method [20] does not require additional modification because the unknown function in the solution domain r £ (0,1] monotonically decreases. In this case, to solve problem (4), (5) in the time interval t £ [0, t*], it is sufficient to solve problem (35) in the segment r £ [a(0)/a (t*), 1].

4. Computational experiment

Verification of the BEM algorithm. Here we test the BEM algorithm developed in Section 2 by comparing the calculation results with the exact solutions presented in Section 3 for various values of parameters.

Example 1. We consider the problem in the case of plane geometry (v = 0) with 7 = 1/3, a = 1,p = 1, n = 0. Parameter 3 is varied: 3 = 0.8; 1; 1.2; 1.4; 1.6; 1.8; 2.

The accuracy of solving problem (4), (5) using the BEM algorithm is tested as follows. The exact travelling wave solution ue(t, x) is found by solving problem (21), (22) using the boundary element method described in Section 3. When the exact solution is found condition (5) that corresponds to this solution is

u\x=o = ue(t, 0). (36)

Next, numerical solution of problem (4), (36) is calculated using BEM. Then it is compared with the exact solution ue(t, x). Numerical and exact solutions for 3 = 1.4 are shown in Fig. 3.

u 4

3

2

1

0

0.0 0.5 1.0 1.5 2.0 *

Fig. 3. Travelling heat wave for 3 = 1.4

To estimate the accuracy of the numerical solution we compare the law of motion of the heat wave front found by the BEM and the law of motion x = pt + n of the exact solution. The relative error in determining the heat wave front is shown in Tab. 1.

Table 1. The relative error in determining the heat wave front

t 3 = 0.8 3 = 1 3 = 1.2 3 = 1.4 3 = 1.6 3 = 1.8 3 = 2

0.5 4.96E-03 1.77E-03 6.29E-04 2.10E-04 7.22E-05 2.19E-05 2.86E-05

1 4.63E-03 1.61E-03 5.85E-04 2.10E-04 8.88E-05 4.97E-05 4.37E-05

The results show that the accuracy of the solution decreases as parameter 3 decreases. The results are not acceptable for 3 < 0.8. This seems to be related to the fact that the term u3-1 (see (7)) has a singularity at u = 0 if 3 < 1. Note that Theorem 1 is also valid only for 3 ^ 1. Nevertheless, one should note that the acceptable numerical results are obtained not only under the conditions of Theorem 1 but also for 0 C 3 < 1.

Example 2. We consider the problem in the case of cylindrical geometry (v = 1) with 7 = 1/3, a =1,3 = 1, R = 1,d =1.

The exact solution ue(t, p) is found by solving problem (35) using the BEM [20]. The boundary condition for problem (4), (5) is

u\p=R = ue(t,R). (37)

Comparison of the BEM solution of problem (4), (37) and the exact solution is shown in Fig. 4.

Thus, the results demonstrate the effectiveness of the developed algorithm for solving the problem of heat wave initiating for a non-linear heat equation with a source.

Fig. 4. Heat wave with exponential front

Traveling wave. The second part of the computational experiment is devoted to the numerical analysis of estimates obtained for the parameters of the solution of problem (21), (22). Tab. 2 shows values of z*,zmax and vmax, as well as boundaries of their interval estimates z_ ^ zmax ^ z+, z_ ^ z0 ^ zmax, v_ ^ vmax ^ v+ (see (32)) for a = p =1 and various Y. One can see that the estimates are relatively accurate for zmax and vmax. The values of zmax are closer to the left border of the interval, and the values of vmax are closer to the right border. For the wave length z* = -z0 the rough estimate was obtained, and further improvement is needed.

The results of calculations show that obtained analytical estimates can be probably improved. So, for all calculations performed (both presented in Tab. 2 and not included

in it) the inequalities \zmax — zM\/\zM\ < 0.05, v_ + 0.7 Av < vmax < v_ + 0.8 Av are valid. Here, zM = (z_ + z+)/2, Av = v+ — v_. It follows from (28) and (29) that if we have values of z0(Y,a, p)\a=^=1, zmax(Y, a, p)\a=^=i and vmax(Y, a, p)\a=M=i, we can find z0(Y,a,p), zmax(Y, a, p) and vmax(Y, a, p) for all a and p as z0(Y,a,p) = = pa_1zo (y, 1, 1), zmax (Y, a, p) = pa_1zmax(Y, 1, 1), vmax(Y, a, p) = p^ a_1 v max (Y, 1, 1).

Table 2. Estimates of travelling wave parameters

Y a M Z- zmax z+ z zo v- vmax v+

1 1 1 -2 -1.547918 -1 -3.414214 -2.328672 0.5 0.858849 1

0.5 1 1 -3 -2.575608 -2 -5.449490 -4.265408 2 2.744421 3

1/3 1 1 -4 -3.590301 -3 -7.464102 -6.229110 4.5 5.643519 6

0.2 1 1 -6 -5.611311 -5 -11.477226 -10.194475 12.5 14.4864942 15

Conclusions

In this study we consider the problem of a heat wave initiating for a non-linear heat equation with a source and construct the solution on a specified finite time interval. We develop a step-by-step algorithm based on the iterative BEM using the dual reciprocity method. We choose

systems of radial basis functions that ensure convergence of the iterative process at each step of the solution. To verify the developed algorithm we use exact travelling wave solutions. Their construction is reduced to solving the Cauchy problem for the ODE with a singularity. For this Cauchy problem we prove the existence and uniqueness theorem of the classical solution that does not have to be analytical. A qualitative study of the solution properties was performed, and some estimates for the amplitude and wave length were obtained. To solve the Cauchy problem numerically, we develop an iterative algorithm based on the BEM. It allows one to determine correctly the boundary of the solution domain where the derivative of the unknown function tends to infinity. The performed calculations show the effectiveness and accuracy of the developed computational algorithm.

Further research can be directed towards expanding the proposed approach to other types of problems and to multidimensional equations. It is also interesting to apply new methods such as the method of fundamental solutions, the interior field method and the null field method to solve the considered elliptic equations and to compare these methods with the BEM.

The study was funded by RFBR (research project no. 20-07-00407) and by RFBR and MOST (research project no. 20-51-S52003.

References

[1] A.A.Samarskii, V.A.Galaktionov, S.P.Kurdyumov, A.P.Mikhailov, Blow-up in quasilinear parabolic equations, Berlin-New York, Walter de Gruyte, 1995.

[2] J.L.Vazquez, The Porous Medium Equation: Mathematical Theory, Clarendon Press, Oxford, 2007.

[3] V.K.Andreev, Yu.A.Gaponenko, O.N.Goncharova, V.V.Pukhnachev, Mathematical models of convection, Berlin, De Gruyter Studies in Mathematical Physics, 2012.

[4] M.Yu.Filimonov, L.G.Korzunin, A.F.Sidorov, Approximate methods for solving nonlinear initial boundary-value problems based on special construction of series, Rus. J. Numer. Anal. Math. Modelling, 8(1993), no. 2, 101-125.

[5] A.L.Kazakov, O.A.Nefedova, L.F.Spevak, Solution of the problem of initiating the heat wave for a nonlinear heat conduction equation using the boundary element method, Comp. Math. Math. Phys., 59(2019), no. 6, 1015-1029. DOI: 10.1134/S0965542519060083

[6] P.K.Banerjee, R.Butterfield, Boundary Element Methods in Engineering Science, London, McGraw-Hill, 1981.

[7] M.Tanaka, T.Matsumoto, S.Takakuwa, Dual reciprocity BEM for time-stepping approach to the transient heat conduction problem in nonlinear materials, Comput. Methods Appl. Mech. Eng., 195(2006), 4953-4961.

[8] Z.C.Li, M.G.Lee, H.T.Huang, J.Y.Chiang, Neumann problems of 2D Laplace's equation by method of fundamental solutions, Appl. Num. Math., 119(2017), 126-145.

DOI: 10.1016/j.apnum.2017.04.004

[9] Z.C.Li, J.Y.Chiang, H.T.Huang, M.G.Lee, The interior field method for Laplace's equation in circular domains with circular holes, Eng. Anal. Boundary Elem., 67(2016), 173-185. DOI: 10.1016/j.enganabound.2016.03.006

[10] M.G.Lee, Z.C.Li, L.Zhang, H.T.Huang, J.Y.Chiang, Algorithm singularity of the null-field method for Dirichlet problems of Laplace's equation in annular and circular domains Eng. Anal. Boundary Elem., 41(2014), 160-172. DOI: 10.1016/j.enganabound.2014.01.013

[11] L.C.Wrobel, C.A.Brebbia, D.Nardini, The dual reciprocity boundary element formulation for transient heat conduction, In: Finite elements in water resources VI, Berlin, SpringerVerlag, 1986, 801-811.

[12] L.V.Ovsiannikov, Group analysis of differential equations, Academic Press, New York-London, 1982.

[13] V.K.Andreev, O.V.Kaptsov, V.V.Pukhnachov, V.V.Rodionov, Applications of group-theoretical methods in hydrodynamics, Dordrecht, Kluwer Academic Publishers, 1998.

[14] A.D.Polyanin, V.F.Zaitsev, Handbook of nonlinear partial differential equations, Boca Raton-London-New York, Chapman and Hall/CRC, 2011.

[15] V.V.Pukhnachev, Exact multidimensional solutions of the nonlinear diffusion equation, J. Appl. Mech. Тechnical Phys., 36(1995), no. 2, 169-176.

[16] A.L.Kazakov, On exact solutions to a heat wave propagation boundary-value problem for a nonlinear heat equation, Sib. Electronic Math. Reports, 16(2019), 1057-1068 (in Russian). DOI: 10.33048/semi.2019.16.073

[17] A.L.Kazakov, P.A.Kuznetsov, A.A.Lempert, On a Heat Wave for the Nonlinear Heat Equation: An Existence Theorem and Exact Solution, Continuum Mechanics, Applied Mathematics and Scientific Computing: Godunov's Legacy, 2020, Springer, 223-228.

DOI: 10.1007/978-3-030-38870-6_29

[18] A.L.Kazakov, L.F.Spevak, Numerical and analytical studies of a nonlinear parabolic equation with boundary conditions of a special form, Appl. Math. Model., 37(2013), 6918-6928. DOI: 10.1016/j.apm.2013.02.026

[19] A.L.Kazakov, Sv.S.Orlov, S.S.Orlov, Construction and study of exact solutions to a nonlinear heat equation, Sib. Math. J., 59(2018), no. 3, 427-441.

DOI: 10.1134/S0037446618030060 Study of Exact Solutions to A Nonlinear Heat Equation. Sib Math J 59, 427-441 (2018). https://doi.org/10.1134/S0037446618030060

[20] A.L.Kazakov, P.A.Kuznetsov, L.F.Spevak, Analytical and numerical construction of heat wave type solutions to the nonlinear heat equation with a source, Journal of Mathematical Sciences, 239(2019), no. 1, 111-122. DOI: 10.1007/s10958-019-04294-x

О построении решений задачи со свободной границей для нелинейного уравнения теплопроводности

Александр Л. Казаков

Институт динамики систем и теории управления имени В.М. Матросова СО РАН

Иркутск, Российская Федерация

Лев Ф. Спевак

Институт машиноведения УрО РАН Екатеринбург, Российская Федерация

Минг-Гонг Ли

Университет Чунг Хуа Город Синьчжу, Тайвань

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

Ключевые слова: нелинейное уравнение теплопроводности, тепловая волна, метод граничных элементов, приближенное решение, точное решение, теорема существования.