On an exact solution to the nonlinear heat equation with a source

The paper is devoted to the study of a singular nonlinear second-order parabolic equation, which is called the porous medium equation or the nonlinear heat equation. One of the important classes of its solutions is heat waves propagating over a zero background with a finite velocity. This property is not typical for parabolic equations and is a consequence of singularity. The main object of study is exact solutions of mentioned type. A new way of separating variables is used to represent them. We obtain conditions when it is possible to make a reduction to the Cauchy problem for an ordinary second-order differential equation with a singularity. It is shown that the Cauchy problem describes a heat wave whose front moves exponentially. We construct a solution to the Cauchy problem as a power series and determine the cases when the series breaks off, i.e. the solution has the form of a polynomial, and the corresponding heat wave can be written explicitly. If the Cauchy problem cannot be explicitly integrated, the solution is constructed numerically. An algorithm based on the boundary element method is proposed. We perform a computational experiment and conclude the properties of the found solutions. Besides, the accuracy of the calculation results is analyzed.


Introduction
We consider the nonlinear second order equation which is a well-known and popular mathematical object, since it has interesting mathematical properties and numerous applications. It is often named as "Generalized porous medium equation" [1], as well as "Nonlinear heat equation with a source" in Russian papers [2]. In physics and mathematical modeling, equation (1) is used to describe processes occurring in high-temperature plasma [3], gas filtration in porous media [4], population dynamics [5], etc. Mathematical properties of equation (1) are studied in a large number of publications. Among them, we point out the pioneer monograph of Russian scientists [6] and the fundamental work of a modern Spanish specialist [1]. One of the interesting mathematical features of equation (1), which is closely related to applications, is that in the case when Φ(0) = Ψ(0) = 0, it has one specific class of solutions. These solutions are non-trivial and non-negative, and they vanish on some manifold: U | a(t,x)=0 = 0, where x = (x 1 , x 2 , x 3 ) is a vector of spacial variables. Such solutions are called  [3] or filtration waves [7]. Apparently, such waves were first discovered in the middle of the last century in connection with studies of processes occurring after a nuclear explosion [3]. Later similar mathematical objects were applied to model gas filtration in a porous medium [4]. Much attention is paid to such solutions in the scientific school of A. F. Sidorov [7]. The method of special series [8], including multiple ones [9], is usually used as a research tool. Among the publications devoted to solutions of this type, we especially note the works of S. N. Antontsev and co-authors [10], in which equation (1) is considered in a general form.
The most common variant in the literature is when Φ(U ) and Ψ(U ) are power functions. Such a case is the main object of consideration in the above-mentioned fundamental books [1,2]. One of the key directions in the development of mathematical physics theory is the construction and study of exact solutions [11]. In connection with the exact solutions to equation (1) with power nonlinearities, the works of A. A. Samarsky and his followers [2,12], as well as our recent papers [13,14] can be mentioned.
The construction of exact solutions allows us to study in more detail the nature of the processes that describe the considered PDEs. Another essential application of exact solutions is the verification of approximate methods for solving boundary value problems. For more than ten years, we have been developing algorithms for solving singular nonlinear parabolic equations in various formulations based on the boundary element method [15,16]. Each class of problems requires verification of algorithms, and, consequently, the corresponding exact solutions. This paper is devoted to the study of one interesting exact solution to equation (1).
If β is a natural number and the function a(t) is sufficiently smooth, problem (2), (3) satisfies the theorem proved in [13]. The theorem ensures the existence and uniqueness of a non-trivial analytical solution that changes the sign for ρ = a(t). The positive part of this solution with trivial u ≡ 0, whose existence for β > 0 is obvious, forms the already mentioned heat wave propagating over the zero background with a finite velocity.
Since the proof of the theorem is constructive and gives a recurrent procedure for constructing a solution in the form of a characteristic series [17], it can be used to verify the results of numerical calculations. We have experience in verifying calculations performed with the boundary element method (BEM) by using segments of series [15]. However, for non-integer values of β > 0, the theorem is not applicable, and the radius of convergence of the characteristic series is usually unknown and small. In this context, the problem of constructing exact solutions with the type of heat wave becomes especially relevant.
Previously, we studied this problem both in the special case α = 0 [18,19] and for the general form of equation (2) [13,14]. We obtained exact generalized self-similar solutions u = u 1 (t)u 2 (ρ/a(t)) and solutions in the form of a generalized traveling wave u = u 1 (t)u 2 (ρ−a(t)) as well. They are constructed as power series, for which estimates of convergence radii are found in some particular cases [18]. The properties of some solutions are analyzed by the methods of 3 the qualitative theory of ordinary differential equations (ODE) [14,19]. Besides, we propose an approach to the construction of generalized self-similar solutions and solutions in the form of a generalized traveling wave by the BEM [20].
However, there are still some open questions of interest. We can consider other ways of separating variables, for example, when the function u 1 = u 1 (ρ), i.e. it depends not on t, as in published articles, but on a spatial variable. In this article, we carry out a study of such solutions, obtain the conditions for their existence, and study their properties.

Reduction to ODE
Let us look for an exact generalized self-similar solution (2), (3) in the form Note that we have not previously considered this method of finding solutions of the desired type. As a result of substitution, we get the equation Next, we change the independent variable, s = − ln z, in equation (4), that arrives us to the equation After reducing similar ones and multiplying both sides by ρ 2 /ϕ 2 (ρ) we get Proof. In order to prove the proposition we equate all coefficients in equation (6) to constants. As a result, we have an overdetermined system that is solvable if and only if conditions a)-c) are satisfied.
Next, we consider a(t) = exp(Bt), since a(t) = B leads to the well-known stationary solution [11]. Under Proposition 1, equation (6) takes the form Let A = B > 0. Then equation (7) is simplified as The first of them follows from (3), since the condition u| ρ=exp(Bt) = 0 leads to v(0) = 0. The second one ensures the consistency of the equation and the existence of a non-trivial solution.
An important advantage of problem (8), (9) is that the equation is autonomous and does not explicitly depend on s. However, in the case of plain symmetry, it is more convenient to refuse substitution s = − ln z. This makes it possible to consider the solution for ρ < 0. It makes sense for ν = 0, but not for ν = 1, 2.

Construction a wave with an exponential front as a series
We consider the case of plane symmetry. To construct and study the solution in the form of a power series, let us return to equation (4).
Let the conditions of Proposition 1 be satisfied. Then the ODE has the form vv For (10) We are interested in solving problem (10), (11) on the interval z ∈ (0, 1]. You can see that (10) has a singularity at the point z = 0. To expand it, we introduce a new function p and a new independent variable ζ: p = z 2 v, ζ = 1 − z.
After substitution, reduction of similar terms and multiplying by (ζ − 1) 4 we get the following equation Recall that A = B > 0, α > 0. Conditions (11) take the form It is easy to see that problem (12), (13) does not have any other singularities besides vanishing at ζ = 0 of the term that multiplies p . So it can be considered in −∞ < ζ < +∞. We construct a solution in the form of the power series The coefficients of series (14) are determined by induction on n. From conditions (13) we know that p 0 = 0, p 1 = σ. To find p 2 , we differentiate both parts of (12) with respect to ζ and put ζ = 0. By resolving the obtained equality with respect to p 2 , we have Similarly, we get that And so on. Assume coefficients up to n are found. To find p n+1 , we differentiate both parts of (12) n times with respect to ζ and put ζ = 0. We get the equality and resolve it with respect to p n+1 . Then we have the formula Since σ > 0, then p n+1 can be uniquely found from (17). Local convergence of the series (14) follows from previously proved theorems (see [13]). An interesting question is the value of the radius of convergence R. If we can show that R ≥ 1, this means that there is an analytical solution to problem (2) (15) it follows that p 2 = 0. It can be shown using (17) that this yields the equalities p 3 = p 4 = . . . = p n = . . . = 0. In other words, the series (14) ends on the second term and the solution to problem (12), (13) has the form p(ζ) = p 0 + p 1 ζ = σζ.

5.
Construction of a wave with an exponential front by the BEM As already noted, the radius of convergence of series (14) is unknown, so the relevant task remains to find solutions to problem (12), (13) on the segment ζ ∈ [0, 1]. Let's represent (12), (13) as Here, P (ζ, p, p ) = p − ζp − (p ) 2 /σ /p+2−α/B, q = ∂p/∂n is a flow, i.e. a normal derivative at the boundary points. Using an approach to solving one-dimensional problems of potential theory [21], we can obtain a solution to problem 1 by the boundary element method. The iterative procedure is as follows Here p (k) 2 and q (k) 2 are k-th iterations of the boundary values of the unknown function and flow at ζ = 1; p (k) (ξ); ξ ∈ (0, 1) is k-th iterations of the solution; u * (ξ, ζ) and q * (ξ, ζ) are fundamental solution and its normal derivative [21]; Decomposing P (ζ, p (k) , [p (k) ] ) by the system of radial basis functions (RBF) [22,23] we transform the integral J (k) in accordance with the dual reciprocity method [24,25], Here, f i (ζ) = f i (|ζ − ζ i |) are RBF whose values depend on the distance between the current point and the specified collocation points ζ 1 , ζ 2 , . . . , ζ m lying on the segment [0, 1]. For each function f i (ζ) there is a function u i that f i = ∂ 2 u i /∂ζ 2 ; q i = ∂ u i /∂n; coefficients α (k) i , i = 1, ..., m, are determined from the system of equations obtained from equality (22) written at the collocation points for (k − 1)-th iteration.
We verified the iterative algorithm described above by comparing it with the exact linear and quadratic solutions given in Proposition 2. It is easy to show that in the first case, for α/B = 1, the exact solution is obtained already at the first iteration. of the BEM solutions for α/B = σ + 2 relative to the quadratic exact solutions for σ = 1 for different numbers of collocation points. Multiquadric functions [26,27] are used as the RBF.
The results show good accuracy of the BEM solutions, as well as convergence relative to the number of collocation points.
Let us now consider case of α/B ∈ {1; σ + 2}. Figure 1 shows graphs of the BEM solutions for σ = 1 and various values of α/B. The convexity of the graphs corresponds to the conclusions made in the previous section, and over the entire considered area. Analysis of numerical solutions showed the following interesting facts.
1. For 1 < α/B < σ + 2, the solution is monotone, and for α/B ≥ σ + 2, it has a maximum at ζ max ∈ (0, 1], which shifts to the left when the parameter α/B increases. 2. For α/B = 2σ + 2, the solution is zero for ζ = 1. If α/B > 2σ + 2, then the positive solution to problem (12), (13) is determined on the segment ζ ∈ [0, ζ 0 ], ζ 0 < 1, and p(ζ 0 ) = 0. In this case, we can construct a solution to problem (2), (3) in a standard way only on a finite , the solution is positive on the segment ρ ∈ [a(0), a(t)). When t ≥ t * it is equal to zero for ρ = a(0), i.e. the heat wave has the form of a solitary wave. Note that solutions to problem (2), (3) probably have the same property for 1 < B/α < 2σ + 2, if we consider them on the all numeric plane. In any case, the solution to equation (19), as can be easily seen, is a solitary wave defined at −∞ < t < +∞, − exp(Bt) ≤ ρ ≤ exp(Bt), that takes zero values at the boundaries ρ = ± exp(Bt). The continuous form of these solutions allows us to find their derivatives, so the accuracy  (12). The largest absolute value of the expression on the left side of equation (12) on the interval (0, 1), i.e. the residual of the equation, is used as an error estimate. Table 2 shows the errors of the BEM solutions for σ = 1 for different numbers of collocation points, as well as similarly calculated errors of solutions in the form of a series for different degrees of its segments. The values in parentheses in the last two lines are errors on the intervals (0, ζ max ). Thus, solutions in the form of a series have satisfactory accuracy only in this part. The results of the calculations show that the proposed algorithm based on the BEM allows us to find solutions to problem (12), (13) with sufficient accuracy. Their error decreases with increasing number of collocation points and weakly depends on the numerical parameters of the problem, in contrast to solutions in the form of series segments.

Conclusions
In this article, we have studied one exact solution, which has the heat-wave type, of the nonlinear heat equation with a source. Its construction has been reduced to the Cauchy problem for a second-order ODE with the singularity. The solution of the Cauchy problem has been obtained as a power series with recursively determined coefficients. We have considered in detail the cases when the series breaks off, and the solution has the form of a polynomial. The result has been interpreted from the point of view of the original problem. We have proposed an algorithm for numerical solution of the Cauchy problem by the boundary element method. The computational experiment has demonstrated good accuracy of numerical solutions and allowed us to perform a numerical analysis of their properties. The most interesting result is that for some values of the input parameters, the solution increases monotonically and indefinitely, while for others it has an extremum and takes positive values on a finite interval. Further research can be related to the analysis of equation (8) by methods of the qualitative theory of differential equations and the determination of the properties of its solutions.