A Uniﬁed Analytical Approach to Fixed and Moving Boundary Problems for the Heat Equation

: Fixed and moving boundary problems for the one-dimensional heat equation are considered. A uniﬁed approach to solving such problems is proposed by embedding a given initial-boundary value problem into an appropriate initial value problem on the real line with arbitrary but given functions, whose solution is known. These arbitrary functions are determined by imposing that the solution of the initial value problem satisﬁes the given boundary conditions. Exact analytical solutions of some moving boundary problems that have not been previously obtained are provided. Moreover, examples of ﬁxed boundary problems over semi-inﬁnite and bounded intervals are given, thus providing an alternative approach to the usual methods of solution.


Introduction
Moving boundary problems occur in numerous important areas of science and engineering [1][2][3]. Recent applications of such problems include modelling biological and tumour invasions [4][5][6][7], drug delivery [8] and melting of crystal dendrite [9]. The classical one-dimensional Stefan problem is the canonical moving boundary problem which was introduced by J. Stefan, a Slovenian physicist, in a series of four papers to model the melting of ice and to calculate the moving boundary where the phase transition from water to ice occurs; see Vuik [10] for some historical notes. Since Stefan's seminal work, there has been extensive research on moving boundary problems in a wide variety of areas. We refer the reader to the books by Gupta [11], Hill [12] and Crank [1], as well as to the comprehensive references therein.
There are many numerical and approximate analytical methods which can be employed to solve moving boundary problems. A comparison of different numerical methods for moving boundary problems was done by Furzeland [13]. One approximate analytical method for one-dimensional Stefan problems is the so-called heat balance integral method (HBIM). This method was first introduced by Goodman [14,15], and which in turn was adapted from the Kármán-Pohlhausen integral method for the study of boundary layers. Schlichting and Gersten's book [16] contains a translation of Kármán and Pohlhausen's work. Other approximate analytical methods for one-dimensional Stefan problems include the refined integral method (RIM) [12,17] and the combined integral method (CIM) [18].
Exact analytical solutions for some one-dimensional Stefan problems are reviewed in [1,12], for example. These are typically similarity solutions where the moving boundary varies as the square root of time. However, in general, such solutions have remained stubbornly difficult to obtain mainly because the moving boundary makes the problem nonlinear even if the underlying system is linear. This means that the standard methods for constructing exact analytical solutions via Green's functions or integral transforms such as those by Fourier or Laplace, amongst others, are usually not applicable to moving boundary problems. It should be noted that exact analytical solutions, while difficult to obtain, are also useful beyond the simple problem they solve since they are key to the validation of numerical or other approximate analytical methods [18][19][20].
On the other hand, the solution of an initial-boundary value problem (IBVP) for the one-dimensional heat equation with fixed boundaries can in many cases be done using (i) separation of variables when the domain is bounded or (ii) integral transforms when the domain is unbounded. However, with non-homogeneous boundary conditions (BCs), the above two methods are generally not applicable directly.
The goal of this article is to propose a unified approach to the solution of the onedimensional heat equation with either a fixed boundary or a moving boundary. The main idea is to embed the partial differential equation (PDE) and the initial condition of the original IBVP into a new initial value problem (IVP) whose spatial domain is all of R but with arbitrary unknown functions. For this new problem, the well-known formula for the IVP of a non-homogeneous heat equation in R is used. Then, the BCs of the original IBVP are imposed to determine the unknown functions. It should be noted that well-posedness of the moving boundary problem is an important issue but is outside the scope of this work. Instead, here we focus on deriving the formal solution of such a problem and give several interesting applications.
This article is organised as follows. In Section 2, we explain the proposed method as outlined above. Then, in Section 3, we illustrate the method with four examples of IBVPs with fixed boundaries on semi-infinite and bounded domains. Section 4 illustrates the method with two examples of Stefan problems. A brief discussion and conclusion are given in Section 5.

Analytical Method
Consider the IBVP The parameter κ > 0 is the diffusion coefficient and a, b, c and d are constants such that |a| + |b| > 0, |c| + |d| > 0.
The functions x ± and g ± depend on t in general while the initial distribution u 0 is a function of x.
We can embed the PDE and initial condition of (1) into the IVP where f (x, t) = 0 for all x − (t) < x < x + (t) and t > 0. More specifically, let where H is the usual Heaviside function (i.e., H(x) = 1 if x ≥ 0 and H(x) = 0 if x < 0) and f ± are functions to be determined so as to satisfy the BCs in (1). Moreover, v 0 is an extension of u 0 in the sense that v 0 (x)| x − (0)≤x≤x + (0) = u 0 (x).
It is well known [21] that the formal solution of the IVP (2) can be expressed as where G is the heat kernel (or Green's function) Using (3), and assuming for simplicity that f ± depend on t but not on x, we obtain Direct calculations yield dy. Letting we have where erfc is the complementary error function with properties Note also that erfc(0) = 1 and erfc(∞) = 0. Thus, (5) becomes Hence, a function u that satisfies the PDE and initial condition of the IBVP (1) is where v 0 is a function such that v 0 (x)| x − (0)≤x≤x + (0) = u 0 (x). Note that f ± have yet to be determined. This is where the third and fourth equations in (1) will be needed. Next, we consider the BCs in (1). It is straightforward using (6) that 3 2 . Differentiating (7) with respect to x gives Substituting the expressions for u(x, t) and u x (x, t) into the BCs in (1), we have Equations (8) and (9) form a pair of coupled linear Volterra integral equations for the time-dependent functions f ± , which cannot be solved analytically in general. Hence, the formal solution of the IBVP (1) is given by (7)-(9).

Remark 1.
Although the extension of u 0 to v 0 is not unique, the formulas for f ± in (8) and (9) get adjusted, according to the choice of v 0 , in such a way that the BCs in (1) are always met.

Remark 2.
A particular case of (1) is a fixed boundary problem when x − (t) = x − and x + (t) = x + are both constants. Then, (8) and ( If we take the Laplace transforms of (10) and (11), and use the convolution property and the properties in (A1) (note that in the following, the formulas (A1), (A2) and (A3) can be found in Appendix A), we obtain respectively, wheref ± (s) andĥ ± (s) denote the Laplace transforms of f ± (t) and h ± (t), respectively. This is a linear algebraic system forf ± (s), which can be solved explicitly. However, the analytical Laplace inversion to get f ± (t) is not going to be straightforward in the general case. Nevertheless, in the next section, we shall illustrate how this can be done in some cases.

Examples: Fixed Boundary Problems
Here, we give illustrative examples of fixed boundary problems for semi-infinite and finite intervals. We remark that these examples are meant to give an alternative derivation of the solutions and, more importantly, to show how fixed and moving boundary problems can both be treated as special cases of our embedding method.

Example 2.
Next, we look at the IBVP This is an example of a fixed boundary problem with a constant flux BC. In this case, we have . As in Example 1, Equation (7) implies (15) and (16). Here, h − (t) = −1 and h + (t) = 0 from (12). Then, (13) becomes Hence, using (A1), we deduce that Equation (16) can be rewritten aŝ Taking the inverse Laplace transform of both sides finally gives as the solution of the IBVP (18). With a slight modification of the BC from u(∞, t) = 0 to u(∞, t) = −1, we can recover the exact solution given in [22] (p. 130) for the temperature during the pre-ablation stage for one-dimensional ablation; see also Example 6 below.
Example 3. Suppose that we have the IBVP This is an example of a fixed boundary problem with a Newton cooling-type BC. We identify . Similar to Examples 1 and 2, Equation (7) implies (15) and (16). We have that h − (t) = 1 and h + (t) = 0 from (12). Then, (13) simplifies to It is clear that f Hence, Having determinedf ± (s), Equation (16) can be rewritten aŝ However, Therefore, the solution of the IBVP (20) is where we have homogeneous Dirichlet BCs. Homogeneous Neumann BCs are treated similarly.
We identify Let v 0 be such that v 0 (x)| 0≤x≤L = u 0 (x) and v 0 (x) = 0 otherwise. Equation (7) gives Equation (12) yields with the help of (A1). The linear algebraic system (13) for this problem is Using (A1), the Laplace transform of (23) is thereforê The following lemma is a useful calculation; see Appendix B for the proof.
It was shown by Rodrigo and Worthy [23] that the inverse Laplace transform of (25) and we recover the well-known solution of the IBVP (22) using separation of variables and the theory of Fourier series. Equation (26) is useful for large times since the series can be truncated to a small number of terms because of the exponentially decaying factors, whereas (23) can be useful for small times and does not require any truncation since it is not in series form.

Examples: Moving Boundary Problems
In this section, we investigate some moving boundary problems using the same methodology as for fixed boundary problems. Note, however, that, in these problems, the boundary function s = s(t) (not to be confused with the Laplace transform parameter) is to be determined together with u = u(x, t).

Example 5. Consider the Stefan problem
where g ± = g ± (t) and u 0 = u 0 (x) are given. We wish to determine u = u(x, t) and s = s(t) that satisfy (27) and also the interface condition for some suitable univariate function F for instance. Following our previous notation, we have Then, (7) yields As u 0 is unspecified, we choose v 0 such that v 0 (x)| 0≤x<∞ = u 0 (x). The left BC (8) is while the right BC (9) is Moreover, from (29), we have The interface condition s (t) = F(u x (s(t), t)) is such that Hence, (28), (30) and (31) provide the determining equations for the time-varying functions f ± and s. Once f ± and s have been determined, u can be calculated from (29), thus obtaining the (formal) exact analytical solution of the Stefan problem (27) and (28).
Suppose that u 0 (x) = u 0 , a constant. Take v 0 (x) = u 0 for all −∞ < x < ∞, which satisfies v 0 (x)| 0≤x<∞ = u 0 (x). Furthermore, take f − (t) = αδ(t), where α is some constant to be determined and δ is the Dirac delta function, and f + (t) = 0. From (29), we obtain (28) and (31) yield respectively. Differentiating the first equation in (33) with respect to t gives A solution to this latter ordinary differential equation is s(t) = 2β √ t for some positive constant β. Substituting this into the equations in (33), we see that β satisfies the relations Thus, β must satisfy the transcendental equation Finally, the solution of (27), (28) is given by where erf is the error function erf(z) = 1 − erfc(z) for z ∈ R and β satisfies the relation (34). This recovers the well-known solution of this special Stefan problem [1,12], typically obtained by a similarity analysis.

Example 6.
Consider a Stefan problem involving a single-phase, semi-infinite, subcooled material. An application arises when determining whether ice melts or water freezes when hot water is thrown over cold ice [24]. This can be formulated as a Stefan problem involving a constant heat source term in the condition at the moving boundary [25]. A related industrial process is known as ablation, where a mass is removed from an object by vapourisation or other similar erosive processes [22,26]. The dimensionless model [20] is given by We want to find u = u(x, t) and s = s(t) that solve (35) and satisfy the interface condition where r is a given positive constant. To our knowledge, the exact analytical solution to this Stefan problem is not known. We identify The left BC (8) is and the (formal) exact analytical solution of the Stefan problem (35), (36) is given by (37).

Discussion and Conclusions
In this article, we proposed a unified approach to solving IBVPs with fixed or moving boundaries. For fixed boundaries, we showed through four examples how the same method is able to handle semi-infinite and bounded intervals. By comparison, IBVPs over semi-infinite and finite intervals are usually solved using separation of variables and integral transforms, respectively. We also derived exact analytical solutions to two examples of physically motivated moving boundary problems. The solutions involved time-varying functions that satisfy exact integral equations which in principle can be evaluated numerically, if not analytically. One possible numerical implementation is to adapt the procedure proposed by Guardasoni, Rodrigo and Sanfelici [27] for solving a Volterra integral equation of the first kind to systems of such integral equations. For instance, in current work by the authors of the present article, we consider the dynamics of outgrowth in a continuum model of neurite elongation and compare the exact analytical solution obtained from the embedding approach and the numerical solution obtained by solving the associated integral equations. As mentioned previously, the exact solution is "formal" in the sense that we assume that the moving boundary problem is well-posed. For example, sufficient hypotheses on the functions g ± , u 0 and F have to be given to ensure that (28), (30) and (31) have a unique (possibly weak) solution.
The embedding approach was introduced by Rodrigo [28] in the context of pricing American put and call options, which leads to a moving boundary problem. It was later applied to fixed boundary problems by Guardasoni, Rodrigo and Sanfelici [27], i.e., for pricing single and double barrier options. In both cases, the underlying equation is the Black-Scholes PDE. A follow-up to the current article is a unified approach to fixed and moving boundary problems for the Black-Scholes PDE. This would yield, in one go, pricing formulas for a variety of financial derivatives, not only American and barrier options. Although the Black-Scholes PDE can be transformed to the heat equation, the transformed initial condition and BCs are usually complicated. Hence, it is more convenient to embed into a final value problem of a non-homogeneous Black-Scholes PDE (instead of an IVP for a non-homogeneous heat equation) and use the Mellin transform techniques developed in [27][28][29][30][31][32].
Two-phase Stefan problems can also be handled with the embedding approach. For example, with a straightforward generalisation of the notation in Example 5, consider the moving boundary value problem where describes the evolution of u (2) = u (2) (x, t). Suppose that g for some suitable bivariate function F. For example, let u (1) (x, t) and u (2) (x, t) represent the temperatures for ice and water at position x and at time t, respectively. Then, x = s(t) is the location of the unknown boundary between the two phases. A typical choice for F is where 1 r is the latent heat of fusion times the density divided by the coefficient of heat conduction [33]. For an arbitrary function s = s(t), we can use (7)-(9) to write down expressions for u (1) (x, t) and u (2) (x, t) that depend on the unknown function s. These expressions for u (1) (x, t) and u (2) (x, t) are then substituted into (40) to obtain an integrodifferential equation involving only s(t). This would complete the specification of the (formal) exact analytical solution of the two-phase Stefan problem (38)-(40).
Our method can be adapted to other PDEs as well. For instance, suppose that we modify the PDE in (1) by including advection and reaction terms of the form where α, β ∈ R are constants and the same initial and boundary conditions in (1) are specified. Then, with the change of variable we obtain the heat (or diffusion) equation Similar initial and boundary conditions for w can also be deduced from the other equations in (1), with adjustments to the functions g ± and u 0 . The results of this article can then be used directly.
The embedding approach given here can also be applied to IBVPs for the wave equation where c is a nonzero constant and (two) initial conditions and BCs analogous to those in (1) are given. We can embed the IBVP into an IVP for a non-homogeneous wave equation defined on R but with arbitrary time-dependent functions. Then, d'Alembert's solution together with Duhamel's principle can be used to solve the IVP. As before, the arbitrary functions are then determined by imposing the BCs. A multilayer diffusion problem was studied in [23] by reformulating it as a sequence of one-layer diffusion problems with arbitrary time-dependent functions. Hence, the crucial step is to solve a one-layer diffusion problem, which is of the form (1). Thus, instead of the more involved Laplace transform approach used for (1) in [23], the solution of the one-layer problem given in (7)-(9) is much easier to derive compared to that given in [23] (Equation (3.22)). In Example 4, we showed the equivalence of the two formulas in the case of homogeneous Dirichlet BCs through (25).
It was mentioned in the Introduction that there are several approximate analytical techniques for one-dimensional Stefan problems (e.g., HBIM, RIM and CIM), all of which are based on the Kármán-Pohlhausen integral method for boundary layers. The basic idea of these approximate analytical methods is to assume in Example 5 that u is some convenient function of s (e.g., u(x, t) is a quadratic expression in x with time-dependent coefficients that depend on s(t)) and that it satisfies a weak form of the IVBP. It is conceivable that the above approximate analytical techniques could be combined with the embedding approach. Indeed, u in (29) already identically satisfies the PDE (and not just its weak form) and initial condition in (5) exactly for any f ± and s. Hence, instead of "guessing" the forms of u and s as in the HBIM, RIM and CIM, we may instead attempt to "guess" the approximate functional forms of s and the auxiliary functions f ± in (28), (30) and (31).
Finally, our embedding approach should also work for IBVPs for the higher dimensional heat (or diffusion) equation where Ω(t) is some domain in R n with a sufficiently smooth boundary with conditions specified there. The analogous IVP to (2) is where v 0 is an extension of u 0 (i.e., v 0 (x) = u 0 (x) for x ∈ Ω(0)) and f = f (x, t) is an arbitrary function such that f (x, t) = 0 for x ∈ Ω(t) and t > 0. The solution v of the above IVP is of course well known and is the higher dimensional analogue of (4). Then, f is determined in such a way that the boundary conditions for u are satisfied. Multilayer diffusion problems in higher dimensions can also be considered. The embedding approach is more straightforwardly adapted compared to the Laplace transform approach followed in [23], which was specific to a one-dimensional diffusion equation.
As can be seen in the discussion here, there are many potential extensions and applications of the embedding method proposed in this article. These are currently works in progress by the first author.
In summary, there are very limited known analytical solutions to Stefan problems and existing closed-form solutions strongly depend on prescribed initial and boundary conditions. As such, numerical simulations are mainly used for the study of moving boundary problems, both for linear and nonlinear equations [4][5][6][7][8][9][34][35][36]. However, a recent study adopted a deep neural network approach to solve Stefan problems [37]. The present article aims to provide a new method to solve one-dimensional Stefan problems for the linear heat equation analytically. This also includes the case where there are advection and reaction terms. Through our embedding method, we obtained a new exact solution to a Stefan problem involving a single-phase, semi-infinite, subcooled material. The exact solution (29) of the Stefan problem considered in Example 5 has also not been obtained previously with this generality. Exact solutions of two-phase Stefan problems can also be found with our approach. This theoretical progress for Stefan problems is a critical step to benchmark numerical results and to further our understanding of the mechanics of melting problems.

Appendix B. Proof of Lemma 1
Proof. We can rewrite (24) aŝ Breaking up the integrals giveŝ However, On the other hand, using the identity where a, b, c = 0, we see from (25) that v(x, s) = x 0 e (−x+y) Hence, we conclude thatû(x, s) =v(x, s) and this completes the proof of the lemma.