On the Analytical and Numerical Study of a Two-Dimensional Nonlinear Heat Equation with a Source Term

: The paper deals with two-dimensional boundary-value problems for the degenerate nonlinear parabolic equation with a source term, which describes the process of heat conduction in the case of the power-law temperature dependence of the heat conductivity coe ﬃ cient. We consider a heat wave propagation problem with a speciﬁed zero front in the case of two spatial variables. The solution existence and uniqueness theorem is proved in the class of analytic functions. The solution is constructed as a power series with coe ﬃ cients to be calculated by a proposed constructive recurrent procedure. An algorithm based on the boundary element method using the dual reciprocity method is developed to solve the problem numerically. The e ﬃ ciency of the application of the dual reciprocity method for various systems of radial basis functions is analyzed. An approach to constructing invariant solutions of the problem in the case of central symmetry is proposed. The constructed solutions are used to verify the developed numerical algorithm. The test calculations have shown the high e ﬃ ciency of the algorithm.


Introduction
We consider a quasilinear parabolic equation in the case of power-law nonlinearity, which models the heat conduction process: Here, t is time, T is the required function (temperature), α, σ are positive constants, H(T) is a specified source function, and H(0) = 0, div, ∇ are spatial divergence and spatial gradient.
Proceeding from the physical meaning of the problem, we consider only the values u ≥ 0 and T ≥ 0. When u > 0 and T > 0, it is clear that Equations (1) and (2), are equivalent. If, simultaneously, T = 0

The Existence and Uniqueness Theorem
We write Equation (2) for the two-dimensional case in the polar coordinates ρ, ϕ as where η(0) = 0. The boundary condition is given as u| ρ=b(t,ϕ) = 0.
Symmetry 2020, 12, 921 3 of 15 We formulate and prove the theorem of the existence and uniqueness of a nontrivial solution to the problems (3) and (4) in the class of analytical functions with a simultaneous construction of the solution in the form of a converging power series. Theorem 1. Assume that the function b(t, ϕ) in the problems (3) and (4) possesses the following properties: 1.
b(t, ϕ) and 1/b(t, ϕ) are analytical with respect to ϕ when −π ≤ ϕ ≤ π and, in terms of t in some neighborhood of the initial time t = 0, η(u) is analytical with respect to u.
Proof. Theorem 1 is proved in two stages as follows: in Stage 1, the solution of the problems (3) and (4) is constructed as a power series with recurrently computed coefficients; in Stage 2, the convergence of the constructed series is proved. (3),

Stage 1. Make a substitution of the independent variables in Equation
With the substitution in (5), the derivatives become as follows: Transform Equation (6) to a more convenient form by grouping the terms: Here, B(τ, ψ, s) = 1 + b 2 ψ /(b + s) 2 . The boundary condition, (4), can then be written as The solution to the problems (7) and (8) is constructed in the form of a series with respect to the powers of s: where u k = ∂ k u ∂s k s=0 . It follows from the boundary conditions that u 0 (τ, ψ) = 0. In order to find u 1 , it is set that s = 0 in Equation (7). We obtain where B 0 = B| s=0 . Assume that u 1 0 (in the opposite case, we have a trivial solution). Then u 1 (τ, ψ) = −σβ(τ, ψ), where β(τ, ψ) = b τ (τ, ψ)/B 0 (τ, ψ). To find u 2 , we differentiate both sides, set s = 0 and represent u 2 in terms of the known quantities. We arrive at the following: and so on. By differentiating Equation (7) k times and setting s = 0, we have where Φ k depends on the coefficients u 0 , u 1 , . . . , u k (and their derivatives with respect to τ, ψ). Being too cumbersome, the expressions for Φ k are omitted here. It can be easily demonstrated that the multiplier for u k+1 in the latter equation is −(kσ + 1)b τ 0; i.e., all the coefficients of the series represented by Equation (9) are determined uniquely.

Stage 2.
Let us now prove the convergence of the series in Equation (9).
All u k (τ, ψ) are analytical in some neighborhood of τ = 0 and at all −π ≤ ψ ≤ π. This follows from the construction procedure and the theorem conditions; namely, Φ k denotes analytic functions by construction since they have been obtained from the differentiation of analytical functions; (kσ + 1)b τ 0 is an analytic function according to condition (2) of the theorem. It also follows from the theorem statement that η(u) = uη * (u), where η * (u) is a function analytical in the neighborhood of u = 0.
The proof follows the classical majorant method [2], which is generally applied to proving the Cauchy-Kovalevskaya theorem and its analogs. In this case, the construction of the majorant is complicated by the fact that the equation we have is not the Cauchy-Kovalevskaya type, since it cannot be solved with respect to the highest derivative. In this connection, we first eliminate the singularity by means of partial expansion of the required function in a Taylor series. Then we construct a majorant equation, unsolved with respect to the highest derivative, but already lacking singularity. Finally, we resolve this equation with respect to the highest derivative by differentiating both sides by the variable s.
The substitution represented by Equation (13) brings the problems (7) and (8) to the form The cumbersome expressions for the known functions g, G i and i = 1, 2, 3, are omitted here. Note, however, that they possess the following properties (provided that the conditions of Theorem 1 are met): 1.

2.
G i and i = 1, 2, 3 are quadratic polynomials with respect to the required function v and its derivatives, whose coefficients are functions of the independent variables τ, ψ, s, analytical in some neighborhood of τ = 0, s = 0 and at all −π ≤ ψ ≤ π. Equation (14) has a unique analytical solution in some neighborhood of the manifold s = 0. There is no need to specify the Cauchy conditions for it when s = 0, since, due to degeneration, Equation (14) is solvable for the only pair of functions v| s=0 , v | s=0 .
The coefficients in Equation (15) are the found by successive differentiation of Equation (14) followed by substitution of s = 0. This results in the following recurrent relationships: . (16) It is obvious that, when k ≥ 2, the right-hand sides in Equation (16) depend on the coefficients v 0 , v 1 , . . . , v k−1 , and this enables the series represented by Equation (15) to be uniquely constructed. Therewith, it follows from the condition of theorem 1 and the properties of analytical functions that all v k (τ, ψ) are analytical when −π ≤ ψ ≤ π and in some neighborhood of the initial time τ = 0.
To prove the existence of the unique analytical solution, majorizing zero, for the problems (18) and (19), one can differentiate the both sides of Equation (18) with respect to s and solve the resulting equation for V sss . The following relationship ensues: where the function Ψ is analytical, majorizing zero, and, in terms of the variable ψ, the analyticity domain contains the whole interval −π ≤ ψ ≤ π.
Equation (24) is of the third order; therefore, to write the Cauchy conditions, one can supplement conditions (19) with a condition for the second derivative V ss | s=0 , which can be easily obtained by setting s = 0 in the both sides of Equation (18). As a result, Equation (24) is resolved with respect to the highest derivative. Its right-hand side and condition (25) are analytical functions; consequently, the problems (24) and (25) has a unique analytical solution. Therewith, it follows from the properties of the right-hand side of Equation (24) that this solution a) majorizes zero; b) in terms of the variable ψ, the analyticity domain contains the whole interval −π ≤ ψ ≤ π. Hence, according to the construction of the problems (24) and (25), it follows that the statement of Theorem 1 is true.

The BEM Solution Algorithm
In the Cartesian coordinates x 1 and x 2 , Equation (3) and the boundary condition (4) have the form Here, the equation b(t, x) = 0 at each moment of time t determines the zero front of a heat wave S (t) , i.e., a closed smooth line bounding the domain V (t) containing the origin; . The problem consists in the determination of the function u = u(t, x) satisfying Equations (26) and (27) where Ω (t) is the spatial domain bounded by S (0) and S (t) .
When condition (27) is met, the following relationship is valid (see Appendix A): where q = ∂u ∂n is heat flux and n is the vector of the external normal to the domain boundary at time t. To solve the problems (26) and (27) on the basis of the boundary element method, we propose the following time-step algorithm. At each (k-th) time step, for t = t k = kh, we solve a boundary value problem for an elliptic-type equation. Namely, in the domain Ω (t k ) , the Poisson equation is considered: with the following boundary conditions corresponding to Equations (27) and (28): The iterative solution of the problems (29)-(31) by the BEM leads to the relationship where u (n) (x) is the n-th iteration of the solution at the step t k , ξ is an internal point of the domain Ω (t k ) , is the fundamental solution of the Laplace equation, and q * (ξ, x) = ∂u * (ξ,x) ∂n [29]. To compute the integral over Ω (t k ) , which is involved in the right-hand side of Equation (32), we apply the dual reciprocity method. The expression 1 where f i denotes radial basis functions (RBF) whose values depend on the distance between the current point and the specified collocation points y 1 , y 2 , . . . , y K lying in the domain Ω (t k ) and on the boundary S k , i.e., f i = f (r i ), where r i = x − y i . In view of Equation (33), the integral over Ω (t k ) can be transformed into an integral over the boundary S k : Here,û i (ξ) is such a function that ∂n ; the coefficients α (n) i are determined on each iteration from the system of linear equations At the first iteration, α In view of Equation (34), the following equation is true for the boundary point x 0 ∈ S k : Since both boundary conditions, (30) and (31), are specified on the boundary S (t k ) , we need to divide this boundary and the boundary S (0) into equal numbers of boundary elements. Dividing each of these boundaries into N rectilinear boundary elements with constant approximation and writing Equation (36) for each node, we arrive at a system of linear equations for the determination of the nodal values of temperature and flux on the boundary S (0) : Here, e m , m = 1, . . . , 2N denotes the boundary elements; x m is the node corresponding to the element e m ; u m and q m are the nodal values of temperature and flux.
Thus, with the use of Equations (32)-(37), the problems (26) and (27) is solved in time steps. At each step, t = t k , the following iteration procedure is performed: We assume u (0) ≡ 0 as the initial iteration; 2.
Solving system (37) The determined nodal values give us the (n+1)-th iteration of the solution

5.
We go to the next iteration.
The iteration process stops with the (n+1)-th iteration, close enough to the n-th one, according to the condition and the (n+1)-th iteration is assumed as the solution at step t = t k , The time derivative u (n) t (y j ) in system (35) at stage 4 is calculated with the use of finite differences as follows: Here, h is the size of the time step, h * j = t k − t * j , and t * j is the instant of time at which the point y j belongs to the zero front, i.e., b t * j , y j = 0. At the first time step, the derivative is calculated by the lower formula in the right-hand side of Equation (41) for all the collocation points y j .

Exact Solutions
We seek the exact solution of Equation (26) in the form u = W(t)U(z), where z = x 2 1 + x 2 2 /a(t). We used similar representations earlier in [17] for similar problems without a source term. Then, the equality U(1) = 0 corresponds to the condition expressed by Equation (27), when Substituting the expression W(t)U(z) for u into Equation (26), reducing similar terms, and multiplying both parts by a/4W 2 , we arrive at the following relationship: The examination of Equation (43) has enabled us to find that the following equalities are sufficient conditions for it to become an ordinary differential equation (ODE) with respect to U(z): Here, γ, λ, µ, ω, and R are constants, γ > 2 and R > 0. Then, Equation (43) becomes Finally, having grouped the terms and added the initial conditions, we obtain the Cauchy problem for the ODE: Here, the condition for the derivative results from the requirements of compatibility of problem (46) and solution nontriviality.
Problem (46) is degenerate when z = 0 since the multiplier at U is zeroed, i.e., the classical existence and uniqueness theorems are inapplicable to it. Nevertheless, according to Theorem 1, it has a unique nontrivial locally analytical solution when γ = 1, 2, . . .. If γ > 1 and is a non-integer, problem (46) does not fall within the scope of the theorem, and the issue of its solvability remains a challenge. We here suppose that the classical (i.e., twice continuously differentiable) solution does exist in this case, although it is not analytical. When γ < 1, it is obvious that problem (46) has no classical solution.
The computational experiment shows that the domain of convergence of an analytical solution to problem (46) in the form of the series which, in fact, is a particular case of series (9), significantly exceeds that in the general case. The issue whether condition (44) is necessary remains beyond consideration so far. It follows from the initial conditions that U 0 = 0 and U 1 = σ/4. The other coefficients of the series (47) are determined by successively differentiating both sides of Equation (46) and subsequently substituting 1 for z.

Test Examples
As Example 1, we consider a sourceless problem (η ≡ 0) with a symmetric zero front of the form (42) when In this case, Equation (46) is integrated in quadratures, and the exact solution (48) has the form where θ = 4 + 4/σ, χ = 4/θ, µ = θ/C, and C is an arbitrary constant. In all the examples with a zero front symmetric relative to the origin, the two-dimensional BEM solutions are also symmetric. Therefore, they are compared with the exact solutions along the axis Ox 1 . Figure 1 shows a comparison between the BEM solution and the exact solution in the form of Equation (50) with the following values of the parameters: σ = 3, R = 1, C = 5. In all the examples with a zero front symmetric relative to the origin, the two-dimensional BEM solutions are also symmetric. Therefore, they are compared with the exact solutions along the axis 1 Ox . Figure 1 shows a comparison between the BEM solution and the exact solution in the form of Equation (50) with the following values of the parameters:    In Example 2, the BEM solution is compared with the exact solution, (48), for the following parameter values: λ = 1, γ = 3, R = 1, µ = 1, ω = 0.5. The comparison is shown in Figure 2. In Example 2, the BEM solution is compared with the exact solution, (48), for the following parameter values:   (27). Therefore, the relative error along   0 S is used to estimate the accuracy of the BEM solution for Example 2.  The calculation results for various values of the parameters agree well with the exact solution, (48). In all the calculation variants, the greatest deviation of the numerical solution from the exact one at different times occurs on the boundary S (0) , and this agrees well with the form of the boundary condition (27). Therefore, the relative error along S (0) is used to estimate the accuracy of the BEM solution for Example 2. Table 1 shows the relative errors of the BEM solutions on the boundary S (0) for various numbers of the boundary elements for h = 0.1 and h = 0.05. The linear functions f i = 1 + r i are used in these calculations as RBFs. The calculations were made for ε = 10 −7 . The calculation results demonstrate a high efficiency in solving problems (26) and (27) by the boundary element method, as well as convergence with an increasing number of boundary elements.

Analysis of Using Radial Basis Functions
A crucial issue of the effectiveness of the presented algorithm is to select a system of RBFs and their parameters as the most suitable for the class of problems under study. To analyze the use of various kinds of RBF, the constructed algorithm is implemented in the form of a program. The analysis is based on comparing the calculation results for the Example 2, considered in the previous section, with the exact solution in the form of Equation (48).
The form of the RBFs, their shape parameters and the time step are taken as analysis factors. Five RBF systems are considered, namely, the polyharmonic splines f i = r i and f i = r 3 i , the linear functions f i = 1 + r i , the scaled linear functions f i = 1 + δr i , and the multiquadric functions f i = 1 + (δr i ) 2 . The analysis is aimed at estimating the effect of the form of the RBF and the value of the shape parameter δ on the solution accuracy.
The solutions obtained at different RBFs and parameter values are fairly close. The solution accuracy estimation results are shown in Table 2. The table shows [31], and δ 2 = 0.8 √ N/D, where D is the diameter of the smallest circle containing all the collocation points [32], are taken as the values of the shape parameter in the scaled linear functions and in the multiquadrics.   The calculations have shown that the solution accuracy increases as the time step decreases, and this testifies to the convergence of the algorithm developed. The comparison among the results of using different RBF systems enables us to make the following inferences. The BEM solutions obtained with the use of all the selected RBFs have small relative errors, so that the algorithm is stable with respect to the choice of RBFs. The best results have been obtained with the use of the multiquadric functions.
Thus, the developed algorithm, based on BEM with the application of the dual reciprocity method, allows one to solve, with a high accuracy, problems of the form expressed by Equations (26) and (27) in a finite time interval.

Conclusions
The paper has discussed the problem of constructing a heat wave for a nonlinear heat conduction equation with a source (sink) in the case of two spatial variables. In order to solve this problem, the following two different approaches are applied: the power series method, enabling one to prove the existence and uniqueness of the solution in the neighborhood of the initial time moment and to construct the solution in the form of a converging power series, and the boundary element method, allowing the problem to be solved on a specified finite time interval. Calculations showing a stable convergence of the algorithm have been made. Exact solutions have been used to verify the numerical method, their construction being reduced to the integration of the Cauchy problem for the second-order ODE with a singularity. The comparison of the BEM calculations with the here-obtained solutions has testified to a practically complete coincidence, which allows the conclusion that the BEM algorithm works correctly.
Further research may involve consideration of other heat wave initiation problems in a two-dimensional statement, as well as an increase in the problem dimensionality.