Analytical Solutions to the Singular Problem for a System of Nonlinear Parabolic Equations of the Reaction-Diffusion Type

The paper deals with a system of two nonlinear second-order parabolic equations. Similar systems, also known as reaction-diffusion systems, describe different chemical processes. In particular, two unknown functions can represent concentrations of effectors (the activator and the inhibitor respectively), which participate in the reaction. Diffusion waves propagating over zero background with finite velocity form an essential class of solutions of these systems. The existence of such solutions is possible because the parabolic type of equations degenerates if unknown functions are equal to zero. We study the analytic solvability of a boundary value problem with the degeneration for the reaction-diffusion system. The diffusion wave front is known. We prove the theorem of existence of the analytic solution in the general case. We construct a solution in the form of power series and suggest recurrent formulas for coefficients. Since, generally speaking, the solution is not unique, we consider some cases not covered by the proved theorem and present the example similar to the classic example of S.V. Kovalevskaya.


Introduction
Parabolic partial differential equations [1] are used in mathematical modeling of various processes studied in the natural sciences. The most common tools for modeling are second-order equations, which are a differential form of physical laws, in particular, Darcy or Fourier ones. Linear models are the most popular, but with the comparative simplicity of research, which is their main advantage, they are not always accurate enough to describe real phenomena. In such cases, nonlinear and, in particular, quasilinear analogs are usually used. The monograph [2] presents the general theory of linear and quasilinear second-order partial differential equations (PDEs), as well as some systems of PDEs.
One of the classes of partial differential equations that allow us to model significant nonlinear effects are singular parabolic equations [3]. The most important example of such an equation is the porous medium equation [4]. In the presence of a source term, it has the form Here σ > 0 is a constant, T is an unknown function, f (T) is a sufficiently smooth specified function. Equation (1) is also called the nonlinear heat (filtration, diffusion) equation with a source term. The power type of the function K(T) = T σ (for the heat equation, it means a coefficient of thermal conductivity) plays a special role. This type of dependency is typical for a large number of applications. For instance, Equation (1) is used to describe the mechanisms of radiative heat conductivity [5,6], the processes of convection [7], diffusion and filtration of ideal polytropic gas in porous media [4,8], as well as the population dynamics process [9]. Sometimes Equation (1) with T = T(t, x) is considered as a specific case of a symmetric equation [10] or the generalized porous medium equation with a given exponent σ = σ(t, x) [11]. Mathematical models based on Equation (1) and its analogues, besides relatively high accuracy, are quite convenient for study, including the analytical one.
We can see that the parabolic type of Equation (1) degenerates at points where the unknown function T = T(t, x) vanishes. As a result, Equation (1) obtains some specific properties. In particular, it has solutions of the heat wave type (filtration wave, diffusion wave) propagating over a zero background with a finite velocity [5,6]. It is known that such properties are typical for hyperbolic equations and, generally speaking, atypical for parabolic ones. Geometrically, a heat wave is two integral surfaces (or hypersurface if we consider two or more spatial variables) of Equation (1) (trivial and non-negative solutions), joined along a certain line (or surface), which is called the wave front. Further, we will call such solutions as HDW-type solutions. Despite substantial physical and geometric interpretations, HDW-type solutions are relatively rare in the literature. We can highlight the fundamental monographs of Ya.B. Zeldovich [5], A.A. Samarskii and co-authors [6], as well as the works of A.F. Sidorov and his followers (see, for example, refs. [12,13]). Statements of boundary value problems on the initiation of filtration waves, as well as methods for constructing the solutions that have HDW-type, in the class of analytical functions, are proposed in [12,13]. The authors of this article also belong to the scientific school of A. F. Sidorov. Previously, we studied the analytical solvability of similar problems in various formulations: one-dimensional [14,15], symmetric [16,17] and non-one-dimensional [18].
A relevant and interesting area of research in partial differential equations is the construction of exact solutions with predetermined properties. There are quite a lot of methods for constructing such solutions that are applicable to Equation (1). We should first mention the group analysis method proposed and developed in the scientific school of L.V. Ovsyannikov [19,20]. We also point out the articles [10,21], in which exact solutions of symmetric Equation (1) are obtained by various variants of the generalized variable separation method. A detailed review of the approaches to construct exact solutions to nonlinear partial differential equations can be found, for example, in the handbook [22]. In the article [23], we found exact HDW-type solutions to Equation (1) with different fronts. The construction is reduced to the integration of ordinary differential equations (ODEs); their qualitative research is carried out.
As a rule, a successful mathematical model is based not on a single equation, but their systems. In particular, equations having type (1) in the case of a single spatial variable x form reaction-diffusion systems where σ, δ > 0 are constants. Such systems describe various processes in biochemistry, physical chemistry, chemical kinetics and thermodynamics [24,25]. For example, unknown functions T = T(t, x) and S = S(t, x) can represent concentrations of effectors (the activator and the inhibitor, respectively), which participate in the reaction [26]. Note that in the literature, you can find systems of a slightly different form, which are also based on second-order parabolic equations with power nonlinearity. In particular, reaction-diffusion systems where the functions f and p also depend on independent variables are considered in [27,28]; systems describing heat and mass transfer are studied in [29]. The objective of this research is the problem of constructing HDW-type solutions of the system in Equation (2) that has a known law of front motion (a problem with a given diffusion front). In the case of a system, it is more appropriate to talk about a diffusion wave, but we will keep the name "HDW-type solution" for simplicity. The power series method, widely used in the theory of differential, integral and operator equations, is chosen as the central research apparatus.
In modern mathematics, the classical power series method [30] usually requires adaptation to the studied problem. At the same time, unknown functions can have a wide variety of types. For example, in the work of N. A. Sidorov [31], the solution of the ODE is given as a power series with coefficients that depend on logarithms. We highlight the methods of characteristic [12] and special [12,13] series, which modifications we apply in this study.
In this paper, we generalize the results on the construction and research of HDW-type solutions for Equation (1), obtained earlier in the scientific school of A. F. Sidorov, to the case of the system in Equation (2). The problem statement is given and discussed in Section 2.
In Section 3, the solutions to the problem with a given diffusion front are constructed as power-law characteristic series with recursively determined coefficients. The proven existence theorem ensures their convergence. The proof is carried out by the majorant method using the Cauchy-Kovalevskaya theorem [30].
In Section 4, we obtain some exact HDW-type solutions of the system in Equation (2) for different fronts. Here, as in some cases (see, for example, refs. [23,32]), the construction of the solution is reduced to integrating ODEs. We point out the articles devoted to the construction of exact solutions of parabolic systems with power nonlinearity. In article [28], the solutions of the two-component reaction-diffusion system, which differs slightly from Equation (2), are constructed by the method of linear defining equations. Multicomponent systems are considered in [27], where solutions have the form of special matrix constructions. However, we could not find any publications deal with exact solutions of Equation (2) that have HDW-type.
Section 5 discusses several significant particular cases that are not subject to the proved statements. We give an example, which is a specialized analogue of the classical counterexample of S. V. Kovalevskaya. It shows that the sufficient conditions for analytical solvability of the considered problem (see Theorem 1 in Section 3) are close to the necessary ones.
Concluding the introduction, we note that the study of HDW-type solutions for systems of equations with power nonlinearity of the form in Equation (2), as far as we know, is performed for the first time.

Formulation
To be able to consider the class of analytical functions for any valid values of the constants σ and δ, we will make the substitution u = T σ , v = S δ in Equation (2). Such transformations are standard [10] and are regularly used by the authors (see [14] and subsequent works). For the resulting system of equations we consider the boundary conditions It is assumed that a(t), F(u, v) and G(u, v) are known sufficiently smooth functions and In addition, let a (0) = 0 and F(0, 0) = G(0, 0) = 0. The problem in Equations (3) and (4) is atypical for second-order systems. First, it contains only one boundary condition for each unknown function, where, besides, the terms multiplying the higher derivatives vanish. Second, it is obvious that the trivial solution u ≡ 0, v ≡ 0 satisfies this problem. Thus, the classical existence and uniqueness theorems cannot be applied in this case. For the problem in Equations (3) and (4), the existence and uniqueness theorem of a nontrivial solution is valid. We formulate and prove it in the next section.

Existence Theorem
For convenience, we introduce a new independent variable z = x − a(t), which allows setting boundary conditions on the coordinate axis. The problem in Equations (3) and (4) takes the form Remark 1. It is easy to show that the substitution z = x − a(t) is non-degenerate (the Jacobian is equal to one). Therefore, Equations (3) and (4) and Equations (5) and (6) are equivalent.

Remark 2.
Analytical function at a point means that it coincides in some neighborhood of this point with its Taylor series expansion.
Proof of Theorem 1. The proof of the theorem is divided into two stages.
1. Constructing a solution in the form of the Taylor series. 2. Proving of convergence of a series by the majorant method.
We construct the solution to the problem in Equations (5) and (6) in the form of Taylor series with respect to powers of the variable z = x − a(t), Since the line x = a(t) is obviously a characteristic for the system in Equation (5), the series in Equation (8) are characteristic [30]. The series in Equation (8) will be constructed when all their coefficients are found. Let us search for them by induction by order of differentiation n.
To obtain u 1 and v 1 we set z = 0 in Equation (5). From Equation (6), we have that u| z=0 (8)), we get a system of two square equations with respect to u 1 and v 1 which has 4 solutions: By the condition of the theorem u z (t, 0), v z (t, 0) ≡ 0, hence, u 1 = −a σ, v 1 = −a δ and the other pairs from Equation (9) are extraneous roots.
Differentiating Equation (5) by z for z = 0, we find the coefficients Thus, the induction base has been created. Now let the coefficients of the series in Equation (8) up to the number n inclusively be known. We differentiate Equation (5) n times by z and set z = 0, then solve the resulting algebraic equations with respect to unknown quantities u n+1 and v n+1 . As a result, we get the expressions For convenience here we use the notation To simplify further formulas and to avoid repeating, we introduce the notation The right-hand sides in Equations (10) and (11) depend on the values that are known by the induction assumption. Thus, the series in Equation (8) are constructed. For u 1 , v 1 ≡ 0, its coefficients are determined uniquely by Equations (10) and (11). Accordingly, the solution is unique.
Let us now move on to the second stage of proof. The proof uses the majorant method, which is close to the classical method of the Cauchy-Kovalevskaya theorem proving [30]. Here, the construction of the majorant is more difficult since the system we have is not the Cauchy-Kovalevskaya type. Therefore, it cannot be solved with respect to the higher derivative. To apply it, we first obviate the singularity by a partial expansion of the unknown function in a power series. Then we remove the singularity by introducing majorant equations and differentiate them by z to resolve this with respect to the higher (third) derivatives. As a result, we obtain the Cauchy problem of the third-order system having Cauchy-Kovalevskaya type.
Before constructing the majorant problem, we make the substitution in the system in Equation (5) The substitution in Equation (12) is a partial Taylor series expansion of the functions u and v. In this case, Equation (6) is always satisfied. After the substitution, we collect like terms and divide the first and second equations in Equation (5) by za σ and za δ, respectively. The problem takes the form Functions f i , i = 0, 1, 2, 3, are found by formulas Functions g i , i = 0, 1, 2, 3, can be found by the similar formulas if we change σ by δ, U by V, V by U and F by G. Here F i,j , G i,j are coefficients of the Taylor series that exist and converge by the assumptions of the theorem.
You can see that the functions f i , g i , i = 0, 1, 2, 3 are analytical in their variables in a neighborhood of the origin. We sequentially differentiate Equation (13) and set z = 0. It leads us to the solution in the form of the Taylor series whose coefficients are determined by the formulas U n (t) = u n+2 (t) (n + 1)(n + 2) , V n (t) = v n+2 (t) (n + 1)(n + 2) .
It is easy to see that the right sides of Equation (15) can be found using Equations (10) and (11). It follows from Equation (15) that the series in Equations (14) and (8) have the same radius of convergence. Thus, the local convergence of the series in Equation (14) ensures that the series in Equation (8) converges in the same neighborhood.
We now construct a majorant for the series in Equation (14). Note that U 0 , V 0 , U 1 , V 1 , f i , g i , i = 0, 1, 2, 3 are analytical functions. Therefore, we can find majorants for them. Let Let us show that the problem is the majorant problem for Equation (13). According to the definition of the majorant problem [30], we need to prove that U n , V n W n , where W n are the Taylor series coefficients The coefficients U n and V n have already been found. As for W n , we determine by a similar procedure, consistently differentiating Equation (16) and setting z = 0. It follows from the known properties of majorant that By similar reasoning, we get V 2 W 2 and V n W n , n ≥ 3. For brevity, we use the notation α n = 1 2(1 + 1/σ) + (4 + 1/σ)n + n(n − 1) .
The coefficients of the series in Equation (17) majorize the coefficients of the series in Equation (14). Finally, we reduce Equation (16) to a problem of the Kovalevskaya type [30]. The necessity for such reduction may also arise in the study of abstract operator equations [33]. We differentiate Equation (16) by z and resolve it with respect to W zzz . It brings the problem in Equation (16) to the form Equation (18) is a Kovalevskaya type problem with analytical data. According to the Cauchy-Kovalevskaya theorem, it has a unique analytical zero-majoring solution. This solution is a majorant for the series in Equation (14) by construction. The local convergence of Equation (14), as already noted, entails the convergence of the series in Equation (8), which is a solution to the problem in Equations (5) and (6).
The analytical solution constructed in the proof is nontrivial. The existence of a trivial solution for the problem in Equations (5) and (6) is obvious. We prove that it is unique, at least in the class of analytical functions. Corollary 1. Let a(t) be an analytical function at t = 0, F(u, v) and G(u, v) be an analytical functions at (0, 0). Let the following identity also hold, Then the problem in Equations (5) and (6) has the unique analytical solution u(t, z) ≡ 0, v(t, z) ≡ 0 for t = z = 0.
Proof of Corollary 1. To prove the corollary, it is enough to construct the solution in the form of the Taylor series and show that all their coefficients will be identically equal to zero. From the conditions in Equations (6) and (19) it follows that u 0 , v 0 ≡ 0 and u 1 , v 1 ≡ 0. We have that u 1 = u 14 , v 1 = v 14 (see Equation (9)). Using the procedure for constructing the solution, as in the proof of Theorem 1, we obtain that U n = 0, n = 2, ... This means that the problem in Equations (5) and (6) has only trivial solution.

Remark 3.
We can assume from continuity that the derivative a (t) preserves the sign in the considered neighborhood of zero. Therefore, on one side of the line x = a(t) solution in Equation (8) is positive, and on the other side it is negative. The positive part of the solution in Equation (8) and the trivial solution u(t, z) ≡ 0, v(t, z) ≡ 0 form a continuous, piecewise smooth (on the line x = a(t) derivatives are discontinuous) HDW-type solution of the problem in Equations (5) and (6).

Exact Solutions
As we have shown previously [23], the construction of HDW-type solutions for Equation (1) can be reduced to Cauchy problems for second-order ODEs in some cases. Usually, such solutions to nonlinear partial differential equations are called exact solutions [22]. If the ODE can be integrated in quadratures [23], then we obtain exact solutions in the traditional sense. In this section, we will find conditions that allow us to perform such a reduction for system (3).
Let, in the problem in Equations (5) and (6), which, recall, is equivalent to Equations (3) and (4), terms that do not contain derivatives have the form F(u, v) = αv γ , G(u, v) = βu γ , α, β, γ > 0 are constants. The problem becomes When γ ∈ N, a(t) is an analytical function at t = 0 and a (0) = 0, then the problem in Equations (20) and (21) comes within Theorem 1 and has a unique nonzero analytical solution.
The following theorem shows that in some particular cases, this solution is exact, and its construction is reduced to the Cauchy problem for a second-order ordinary differential equation.
Then the problem in Equations (20) and (21) allows being reduced to the Cauchy problem for a system of second-order ODEs (1) for a(t) = c 1 e c 3 t , if γ = 1, (2) for a(t) = c 3 ln (c 1 t + c 2 ), if γ = 2, (3) for a(t) = (c 1 t + c 2 ) (γ−2)/(2γ−2) , if γ ≥ 3. (20) and (21) has a unique analytical solution in the form of the series in Equation (8) with respect to the powers of z and with coefficients depending on t. Let us show the possibility of reduction to the Cauchy problem for the system of ODEs by analyzing the structure of the series in Equation (8).

Proof of Theorem 2. Theorem 1 ensures that the problem in Equations
Spatial variable substitution s = (a(t) − x)/a(t) is successfully applied in the construction of exact HDW-type solutions of nonlinear parabolic equations [17,23]. Such a replacement, first, allows you to transform the condition on the diffusion front into a condition at s = 0, and second, it corresponds to the structure of characteristic series [12,17,30]. This makes it possible to further perform the separation of independent variables.
We introduce a new independent spatial variable s = (a(t) − x)/a(t) = −z/a(t). The second independent variable does not change; we keep the notation t for it. The problem in Equations (20) and (21) takes the form Based on the type of the system in Equation (22) and taking into account that F(v), G(u) are power functions, we make (by analogy with [17]) one more substitution u = a θ (t)p(t, s), v = a θ (t)q(t, s).
To eliminate variable t from the terms multiplying p, p s and q, q s , we should choose such a(t) that is a solution to the following ODE where c = 0 is a constant. Note that the case c = 0 leads to a trivial solution.
Functions a(t) = c 1 e c 2 t (for γ = 1) and a(t) = (c 1 t + c 2 ) (γ−2)/(2γ−2) (for γ ≥ 3) are solutions to Equation (27). Therefore, for all natural γ = 2, the system in Equation (26) can be rewritten as It remains to verify that p t = 0 and q t = 0, in other words, p = p(s) and q = q(s). In fact, it follows from the proof of Theorem 1 that there is a unique nontrivial analytical solution to the problem in Equations (28) and (25) that can be written as series Their coefficients are determined by the following formulas: Let us show that all the coefficients are constant. It is easy to see that p 0 , p 1 , p 2 , q 0 , q 1 , q 2 do not depend on t. Setting n = 2 in Equations (30) and (31), we obtain that p 3 , q 3 do not also depend on t. The base of induction has been created. We now assume that all coefficients p 3 , ..., p n and q 3 , ..., q n have the same property. Then the terms p n a 2−2γ 2−γ and q n a 2−2γ 2−γ vanish. Therefore, p n+1 and q n+1 are also independent of t. In accordance with the principle of mathematical induction, all coefficients of the series in Equation (29) are constant, and the unknown functions have the form p = p(s), q = q(s).
Adding initial conditions p(0) = p 0 , q(0) = q 0 , p (0) = p 1 , q (0) = q 1 that provide nontriviality of the solution, we get the Cauchy problem Its solution helps us to construct a solution to the original the problem in Equations (20) and (21) in the form where a(t) = c 1 e c 2 t (for γ = 1) and a(t) = ( Let us now consider the special case γ = 2. In this case F(u, v) = αv 2 , G(u, v) = βu 2 , so We make the substitution u = a (t)p, v = a (t)q in Equations (20) and (21) and get a condition for a(t), which ensures p = p(z), q = q(z). The problem in Equations (20) and (21) takes the form 1 p(t, s)| s=0 = q(t, s)| s=0 = 0.
To eliminate variable t from the terms multiplying p and q we should choose such a(t) that is a solution to the following ODE a = c(a ) 2 , c = 0 − const, The case c = 0 will be considered separately. It is easy to make sure that the function a(t) = c 3 ln(c 1 t + c 2 ) is the solution to Equation (36) and c = −1/c 3 . Therefore, the system in Equation (34) can be written as 1 a p t + cp = p z + pp zz + 1 σ p 2 z + αq 2 , 1 a q t + cq = q z + qq zz + 1 The solution to the problem in Equations (37) and (35) is constructed as a series with the coefficients F n and G are obtained from (32) for γ = 2. By analogy with the general case, we can show that all coefficients of the series in Equation (38) are constant, and the unknown functions have the form p = p(z), q = q(z). Initial conditions p(0) = p 0 , q(0) = q 0 , p (0) = p 1 , q (0) = q 1 that provide nontriviality of the solution bring us to the Cauchy problem By its solution, we find a solution to the original the problem in Equations (20) and (21) for where a(t) = c 3 ln(c 1 t + c 2 ).
Thus, all the cases specified in the theorem condition are considered. The construction of the solution to the problem in Equations (20) and (21) is reduced to the integration of Cauchy problems for ODEs having the form either Equation (33) in cases 1 and 3 or Equation (39) in case 2.
Next, we consider the case of a simple wave that has a constant propagation velocity. The constraints on the functions F(u, v), G(u, v), in this case, can be significantly weakened compared to Theorem 2.

Proof of Proposition 1. Consider the problem
which is a particular case of Equations (5) and (6), where a(t) is a linear function and z = x − a 1 t + a 0 . All the conditions of Theorem 1 are satisfied; hence the problem in Equations (40) and (41) has a unique solution in the form of the series in Equation (8). We rewrite Equations (10) and (11) obtained for coefficients in the general case, taking into account the conditions of Proposition 1: Following the procedure used in the proof of Theorem 2, we can show that all the coefficients of the series in this case are constant. Therefore, the unknown functions have the form u = u(z), v = v(z).
Taking into account the conditions for derivatives at z = 0, which provide nontriviality of the solution, the problem in Equations (40) and (41) reduces to the following Cauchy problem: It is easy to see that the case γ = 2, c = 0, which was omitted in the proof of Theorem 2, leads to such a problem.
The study of the main case that leads to HDW-type solutions is finished. In the next section, for completeness, we consider cases when none of the conditions in Equations (7) and (19), are satisfied, i.e., u 1 = u 12 , v 1 = v 12 and u 1 = u 13 , v 1 = v 13 .

Special Case
Assume now u 1 = u 12 , v 1 = v 12 (see Equation (9)). Case u 1 = u 13 , v 1 = v 13 is considered similarly up to the notation and therefore is omitted.

Proposition 2.
Let a(t) be an analytical function at t = 0, F(u, v) and G(u, v) be an analytical functions at (0, 0). Let the following condition hold Then the problem in Equations (5) and (6) has a unique solution for t = z = 0 in the form of formal power series.
Assuming that G u (0, 0) = 0, we can construct a non-trivial solution to the problem. The other coefficients of the series in Equation (43) are as follows: It is easy to see that all coefficients of the series in Equation (43) are determined uniquely.

Remark 4.
If G(u, v) = G(v), that is, the second equation of the system in Equation (5) does not depend on the unknown function u, then G n = 0, n = 1, 2, . . . and, therefore, v ≡ 0.
The proven proposition leaves open the question of convergence of the series in Equation (43). The research has shown that it is impossible to prove a general existence and uniqueness theorem similar to Theorem 1 in this case. The convergence or divergence of the series is determined by the parameters of the problem. Next, we give examples confirming this.
We proceed to the limit as n → ∞ and get Substituting the last inequality in Equation (48) we get that R = 0. In other words, the series in Equation (47) converges only if z = 0, and the problem does not have an analytical solution.

Example 2.
Here we consider the problem in Equations (44) and (45) in the case u 1 = u 13 ≡ 0, v 1 = v 13 = δc. Then, according to Remark 4, the problem u t + cu z = uu zz + 1 σ u 2 z , u(t, z)| z=0 = 0, under the condition u 1 ≡ 0 has only a trivial solution u ≡ 0. Then, Equation (44) takes the form It, as shown in Example 1, has a unique non-trivial solution v = δcz. The original problem in Equations (44) and (45) is uniquely solvable in the class of analytical functions.
Thus, we can see that in the particular case discussed in this section, the solution in the form of a formal series can always be constructed. However, the convergence directly depends on the specific parameters of the problem.

Remark 5.
In the considered case, Example 1 is an analog of the classical counterexample of S. V. Kovalevskaya for the linear heat equation [30], which was presented 150 years ago. Just as the counterexample of S. V. Kovalevskaya emphasizes the significance of all the conditions of the Cauchy-Kovalevskaya theorem, Example 1 shows the importance of all the conditions of Theorem 1 proved in this paper.

Conclusions
Summarizing the study, we highlight that the obtained results generalize the research on the construction of HDW-type solutions for nonlinear parabolic equations with singularity, performed earlier in the scientific school of A. F. Sidorov, for the case of a parabolic system with power nonlinearity. We have considered the problem with a given diffusion front. The solution has been constructed as two characteristic series with recursively determined coefficients. We have proved the convergence of series by the majorant method using the Cauchy-Kovalevskaya theorem. This means that the existence and uniqueness theorem of a nontrivial solution in the analytical functions class has been proved.
The obtained non-zero solution with a trivial one forms a piecewise analytical heat wave. We have considered special cases where the construction of HDW-type solutions has been reduced to the Cauchy problem for a second-order ordinary differential equation.
It has been shown that the solution can also be constructed in the form of formal Taylor series if the conditions of the proved theorem are not satisfied. However, Example 1, which is an analog of the well-known counterexample of S. V. Kovalevskaya, shows that it is impossible to prove a more general statement about the convergence of the specified series, as compared to Theorem 1.
Solution representations in the form of series can be directly applied for numerical calculations [12]. However, the series convergence domain is the neighborhood of the point t = 0, and this neighborhood is usually small. In this regard, we plan using an original approach based on the boundary element method (BEM) to solve the problem of the diffusion wave movement at a given time interval. Previously, we successfully applied the BEM in combination with the dual reciprocity method for constructing HDW-type solutions for the porous medium equation [15,16]. The main point here is the elimination of singularity in the multipliers at the highest derivatives; this can be done by the series constructed in this paper.
Further research in this field may be associated with the following directions. The first one is the proof of existence and uniqueness theorems for the problem of heat wave initiation by a boundary regime for the system in Equation (3). The second one is the development of a numerical-analytical method based on Theorem 1 and the boundary element method [15,16]. The third one is increasing the number of spatial variables in the problem [18].
Systems of nonlinear parabolic equations with singularity will further be applied to study natural processes occurring in Lake Baikal. This will allow us to improve the accuracy and adequacy of modeling and forecasting. Additionally, the modeling and simulation will help us to study this system in more detail that, in turn, assists in the protection of the environmental system of the lake included in the list of UNESCO world heritage sites. More specifically, the developed analytical methods in combination with computational technologies, which are planned to be developed, will be applied to modeling polluting substances near settlements, changes in the thickness of the ice and the evolution of the Baikal biota populations.

Conflicts of Interest:
The authors declare no conflict of interest.