Application of the Homotopy Method for Fractional Inverse Stefan Problem

Abstract

The paper presents an application of the homotopy analysis method for solving the one-phase fractional inverse Stefan design problem. The problem was to determine the temperature distribution in the domain and functions describing the temperature and the heat flux on one of the considered area boundaries. It was demonstrated that if the series constructed for the method is convergent then its sum is a solution of the considered equation. The sufficient condition of this convergence was also presented as well as the error of the approximate solution estimation. The paper also includes the example presenting the application of the described method. The obtained results show the usefulness of the proposed method. The method is stable for the input data disturbances and converges quickly. The big advantage of this method is the fact that it does not require discretization of the area and the solution is a continuous function.

1. Introduction

The mathematical models with derivatives of the fractional order have recently found an application for modeling various kind of phenomena in physics, mechanics and economy (see for example [1,2,3,4,5,6]). In the literature, various definitions of the fractional derivatives can be found, the most commonly used are the definitions of Riemann-Liouville, Grünwald-Letnikov, Caputo or Riesz (for definitions see [6,7]). The Caputo derivative has found more applications in recent years. Its advantage is the fact that its Laplace transform contains only the initial values of the integer order derivatives. On the contrary, in case of the fractional derivative of Riemann-Liouville type, its Laplace transform usually contains the initial values of the fractional derivatives for which it is difficult to find a satisfactory physical interpretation. The mathematical theory for the mentioned derivatives is presented, for example, in the monographies [6,7].

In case of the heat conduction in the porous and composite materials it is justified to use the mathematical models with the derivatives of the fractional order [8,9,10]. Voller in the paper [11] demonstrates the usability of such models for modeling the thermal processes in the porous materials and presents two examples: the stationary problem of heat conduction and material melting illustrating the so-called anomalous heat conduction. This anomalous behavior may occur in case there are any impurities in the considered area which are the subareas with greater or less thermal conductivity than the rest of the area. We may expect this when the distribution of these impurities is chaotic (fractal). An example of such situation was also presented by Sierociuk et al. [12]. Brociek et al. [13] compare selected mathematical models for the inverse heat conduction problem on the basis of the measurements taken for the porous aluminum. In particular, the models with the fractional derivatives were compared with the models that use the classic derivative. The fractional derivatives used in the models were of Caputo and Riemann-Liouville type. By solving the proper inverse problems the considered models were compared, or rather their ability to fit into the measurement data. The obtained results show that the models with the fractional derivatives are a better match for the measurement data than the models with the classic derivatives. For the presented case, the best turned out to be the model with the fractional derivative of a Riemann-Liouville type.

The Stefan problem denotes the wide class of mathematical models describing the thermal processes characterized by the phase transitions. Group of such processes includes, for example, solidification of metals, creation of crystals, formation of igneous rock, freezing of food, freezing of water, melting of ice, molecular diffusion, solute transport and others. Solving the Stefan problem consists in determining the distribution of the temperature in considered region and, simultaneously, in determining the location of freezing front dividing the given region into subregions. In recent years, the generalization of this issue has been investigated in terms of application of the fractional derivatives. This type of model has found an application for describing the movement of the shoreline in a sedimentary ocean basin [14,15], the controlled release of a drug from slab matrices [16,17] and the heat conduction in the porous materials [11,18].

The analytical solution of the Stefan problem with the fractional derivative is known only for the simple one-dimensional case and moreover it is defined in semi-infinite region [16,19,20,21]. The papers describing the numerical methods of solving this type of problem are so far limited and apply to one-dimensional and one-phase cases [22,23]. The method for solving the two-phase direct problem is presented by Błasik in the paper [24].

The method using homotopy was previously used by Rejev et al. [15,25] to solve the direct classic Stefan problem and by Hetmaniok et al. [26,27] to solve the inverse problem. There are no papers, however, which would consider the inverse Stefan problem with the fractional derivative. In this paper we will present an application of the homotopy analysis method for solving the one-phase fractional inverse Stefan design problem. The problem consists in determining the temperature distribution in the domain and functions describing the temperature and the heat flux on one of the area’s boundaries.

In the homotopy analysis method, the solution is sought as a series. After applying the method, it turns out that the elements of this series must satisfy some differential equation. The form of this equation is a consequence of the problem being solved. It turns out that if the series converges then its sum is a solution of the discussed equation. The sufficient condition of the convergence is also given as well as the error of the approximate solution estimation. An example illustrating the use of this method is also presented.

2. Formulation of the Problem

We are going to consider the problem defined in the area $\mathsf{\Omega }=\{(x,t);x\in [0,\xi \left(t\right)],t\in [0,t^{*})\}$ (see Figure 1). We split the domain’s $\mathsf{\Omega }$ boundary into three parts:

$$\begin{array}{cc}\hfill {\Gamma }_0& =\left\{(x,0);x\in [0,\xi (0\left)\right]\right\},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\Gamma }_1& =\left\{(0,t);t\in (0,t^{*})\right\},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\Gamma }_2& =\left\{(\xi \left(t\right),t);t\in (0,t^{*})\right\},\hfill \end{array}$$

(3)

where function $\xi$ describes the position of the moving interface.

In the domain $\mathsf{\Omega }$, we consider the following equation:

$$D^{\alpha }u(x,t)=a\frac{{\partial }^{2}u(x,t)}{\partial x^2},$$

(4)

where $\alpha \in (0,1)$ is the order of the derivative of Caputo type, $a=\overline{w}\overline{a}$ is the scaled (fractional) thermal diffusivity [${\mathrm{m}}^2/{\mathrm{s}}^{\alpha }$], that is the thermal diffusivity $\overline{a}$ [${\mathrm{m}}^2/\mathrm{s}$] multiplied by the scaling constant $\overline{w}$ with the numerical value of one and unit [${\mathrm{s}}^{\alpha -1}$], and u, t and x refer to the temperature, time and spatial location respectively. On the boundaries ${\Gamma }_0$ and ${\Gamma }_2$ we set the following conditions:

$$\begin{array}{cc}\hfill u(x,0)& =\phi \left(x\right),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill u\left(\xi \right(t),t)& =u^{*},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \kappa D^{\alpha }\xi \left(t\right)& =-k\frac{\partial u(x,t)}{\partial x}{|}_{x=\xi \left(t\right)},\hfill \end{array}$$

(7)

where k is the thermal conductivity [$\mathrm{W}/\left(\mathrm{m}\mathrm{K}\right)$], $\kappa =\overline{w}\overline{\kappa }$ is scaled (fractional) latent heat of fusion per unit volume [$\mathrm{J}{\mathrm{s}}^{\alpha -1}/{\mathrm{m}}^3$], that is the latent heat of fusion per unit volume $\overline{\kappa }$ [$\mathrm{J}/{\mathrm{m}}^3$] multiplied by the scaling constant $\overline{w}$ with the numerical value of one and unit [${\mathrm{s}}^{\alpha -1}$], $u^{*}$ is the phase change temperature [$\mathrm{K}$]. On the boundary ${\Gamma }_1$ we do not set any boundary condition as we are going to reconstruct it in the inverse problem.

The fractional derivative used in the above equations is of Caputo type. For $\alpha \in (0,1)$ this derivative is defined as follows [6]:

$$D^{\alpha }f\left(t\right)=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{f^{\prime }\left(t\right)}{{(t-s)}^{\alpha }}ds,$$

where $\Gamma (·)$ is the gamma function.

In case of the direct Stefan problem we know all of the functions and parameters describing the initial and boundaries conditions, as well as the material physical parameters and we are looking for the temperature distribution and the position of the moving interface. In case of the inverse problem, in turn, we do not know a part of the input information, e.g., the function describing the boundary conditions, but we do know, however, something about the solution. When considering the Stefan inverse problem, this known information is more often than not, the temperature in the chosen points of the domain or the position of the moving interface. The inverse problems for which the position of moving interface is known are the so-called inverse design problems.

The presented inverse Stefan problem consists in finding a function u which describes the temperature distribution in domain $\mathsf{\Omega }$ and reconstruct functions $\theta$ and q describing the temperature and the heat flux on the boundary ${\Gamma }_1$, which will satisfy Equations (4)–(7). We assume that all of the other functions and parameters are known.

3. Solution of the Problem

For solving the discussed problem, the homotopy analysis method will be used. This method was elaborated by the Chinese mathematician Shijun Liao and is designed for solving various kind of linear and nonlinear operator equations. For the first time the method was described in 1992 in its author PhD thesis [28] and has found an application in various fields since then. It allows us to solve the operator equation:

$$N\left(u\right(x,t\left)\right)=0,(x,t)\in \mathsf{\Omega },$$

(8)

where N is operator (in particular it can be the nonlinear operator), whereas u is unknown function. First we define the homotopy operator $\mathcal{H}$ as [29]:

$$\mathcal{H}(\mathsf{\Phi },p)\equiv (1-p)L\left(\mathsf{\Phi }(x,t;p)-u_0(x,t)\right)-phN\left(\mathsf{\Phi }(x,t;p)\right),$$

(9)

where $p\in [0,1]$ is an embedding parameter, $\mathsf{\Phi }$ is an auxiliary function, $h\ne 0$ the convergence control parameter, $u_0$ an initial approximation of the solution of problem (8) and L an auxiliary linear operator.

Considering the equation $\mathcal{H}(\mathsf{\Phi },p)=0$, we get the so-called zero-order deformation equation:

$$(1-p)L\left(\mathsf{\Phi }(x,t;p)-u_0(x,t)\right)=phN\left(\mathsf{\Phi }(x,t;p)\right).$$

(10)

For $p=0$ we have $L(\mathsf{\Phi }(x,t;0)-u_0(x,t))=0$, from which we obtain $\mathsf{\Phi }(x,t;0)=L^{-1}\left(L\left(u_0(x,t)\right)\right)$. In turn, because for $p=1$ we get $N\left(\mathsf{\Phi }\right(x,t;1\left)\right)=0$, then $\mathsf{\Phi }(x,t;1)=u(x,t)$, where u is the sought solution of the Equation (8). This way the change of the parameter p from zero to one is corresponding to the change from the trivial problem to the original problem.

By expanding the $\mathsf{\Phi }$ into the Maclaurin series with regards to the parameter p and transforming it we get:

$$\mathsf{\Phi }(x,t;p)=u_0(x,t)+\sum _{m=1}^{\infty }{u}_m(x,t)p^m,$$

(11)

where

$$u_m(x,t)=\frac{1}{m!}\frac{{\partial }^m\Phi (x,t;p)}{\partial p^m}{|}_{p=0}={\mathcal{D}}_m\left(\mathsf{\Phi }(x,t;p)\right),m=1,2,3,\dots ,$$

(12)

where ${\mathcal{D}}_m$ is mth-order homotopy-derivative operator [29]:

$${\mathcal{D}}_m\left(f(x,t;p)\right)=\frac{1}{m!}\frac{{\partial }^{m}f(x,t;p)}{\partial p^m}{|}_{p=0}.$$

If the series that occurs in the Formula (11) is convergent in the proper area, then for the $p=1$ we obtain the sought solution:

$$u(x,t)=\sum _{m=0}^{\infty }{u}_m(x,t).$$

(13)

As the approximate solution we can choose the partial sum of the above series:

$${\widehat{u}}_n(x,t)=\sum _{m=0}^n{u}_m(x,t).$$

(14)

To find the exact or approximate solution of the discussed equation it is necessary to determine the function $u_m$. For this we need to differentiate the left and right side of the Formula (10) m-times, with regards to the parameter p, then we divide the obtained result by $m!$ and substitute $p=0$ obtaining the formula that is called the mth-order deformation equation ($m>0$):

$$L\left(u_m(x,t)-{\chi }_m{u}_{m-1}(x,t)\right)=h{\mathcal{D}}_{m-1}\left[N\left(\sum _{i=0}^{\infty }{u}_i(x,t)p^i\right)\right],$$

(15)

where

$${\chi }_m=\left\{\begin{array}{cc}0\hfill & m\le 1,\hfill \\ 1\hfill & m>1.\hfill \end{array}\right.$$

(16)

By the proper choice of the h parameter, which was introduced in the Formula (9), we can impact the region of convergence of the series (13) and its convergence rate [29,30,31]. This is why it is called the convergence control parameter [32,33,34]. One of the methods of choosing the value of this parameter is the so-called “optimization method” [29]. In this method the residue (deviation) function is defined for the considered equation and determined approximate solution ${\widehat{u}}_n$:

$$E_n\left(h\right)=\underset{\mathsf{\Omega }}{\int }{\left(N\left[{\widehat{u}}_n(x,t)\right]\right)}^{2}dxdt.$$

(17)

The optimum value of the convergence control parameter is obtained by finding the argument for which the $E_n$ function reaches its minimum.

The effective region of the convergence control parameter is also defined as follows:

$$\mathbf{R}=\left\{h:\underset{n\to \infty }{lim}{E}_n\left(h\right)=0\right\}.$$

(18)

If we choose the value of the convergence control parameter different than the optimum but still from the effective region $\mathbf{R}$, the obtained series will also be convergent but the rate of the convergence will be lower. The version of the method including the described above selection of the optimum value of the convergence control parameter is called the basic optimal homotopy analysis method [29].

In order to speed up the computation, Liao [29] suggests to substitute the integral in the Formula (17) by its approximate value using the appropriate numerical integration method. In the examples demonstrated by him, the obtained values of the convergence control parameter did not differ much from the values obtained by the application of the Formula (17).

In the discussed problem the N operator has the form:

$$N\left(u\right)=a\frac{{\partial }^{2}u}{\partial x^2}-D^{\alpha }u.$$

(19)

However, as the linear operator L we can choose the operator of form:

$$L\left(u\right)=\frac{{\partial }^{2}u}{\partial x^2}.$$

(20)

Assuming that the series ${\displaystyle \sum _{i=0}^{\infty }}{u}_i(x,t)p^i$ is convergent then for the derivative of a Caputo type we get:

$$\begin{array}{c}{D}^{\alpha }\left(\sum _{i=0}^{\infty }{u}_i(x,t)p^i\right)=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{\partial }{\partial t}\left(\sum _{i=0}^{\infty }{u}_i(x,t)p^i\right)\frac{1}{{(t-s)}^{\alpha }}ds=\hfill \\ =\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\left(\sum _{i=0}^{\infty }\frac{\partial u_i(x,t)}{\partial t}{p}^i\right)\frac{1}{{(t-s)}^{\alpha }}ds=\\ \hfill =\sum _{i=0}^{\infty }{p}^i\left(\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{\partial u_i(x,t)}{\partial t}\frac{1}{{(t-s)}^{\alpha }}ds\right)=\sum _{i=0}^{\infty }{D}^{\alpha }{u}_i(x,t)p^i.\end{array}$$

(21)

By using it we obtain:

$$\begin{array}{c}{\mathcal{D}}_{m-1}\left(D^{\alpha }\left(\sum _{i=0}^{\infty }{u}_i(x,t)p^i\right)\right)={\mathcal{D}}_{m-1}\left(\sum _{i=0}^{\infty }{D}^{\alpha }{u}_i(x,t)p^i\right)=\hfill \\ =\frac{1}{(m-1)!}\left(\frac{{\partial }^{m-1}}{\partial p^{m-1}}\left(\sum _{i=0}^{\infty }{D}^{\alpha }{u}_i(x,t)p^i\right)\right){|}_{p=0}=\\ \hfill =\frac{1}{(m-1)!}\left((m-1)!D^{\alpha }{u}_{m-1}(x,t)+\sum _{i=m}^{\infty }w\left(i\right)D^{\alpha }{u}_i(x,t)p^{i-m+1}\right)){|}_{p=0}=D^{\alpha }{u}_{m-1}(x,t),\end{array}$$

(22)

for $m=1,2,\dots$ and where $w\left(i\right)\in \mathbb{N}$ for $i=m,m+1,\dots$. Likewise we get:

$${\mathcal{D}}_{m-1}\left(a\frac{{\partial }^2}{\partial x^2}\left(\sum _{i=0}^{\infty }{u}_i(x,t)p^i\right)\right)=a\frac{{\partial }^2{u}_{m-1}(x,t)}{\partial x^2}.$$

(23)

Using this in mth-order deformation Equation (15) for $m=1$, (19) and (23), we obtain:

$$L\left(u_1(x,t)\right)=h\left(a\frac{{\partial }^2{u}_0(x,t)}{\partial x^2}-D^{\alpha }{u}_0(x,t)\right),$$

(24)

and for $m>1$:

$$L\left(u_m(x,t)\right)=L\left(u_{m-1}(x,t)\right)+h\left(a\frac{{\partial }^2{u}_{m-1}(x,t)}{\partial x^2}-D^{\alpha }{u}_{m-1}(x,t)\right).$$

(25)

Taking into account the L operator we obtain the set of partial differential equations that allow us to determine the functions $u_m$. For $m=1$ we have:

$$\frac{{\partial }^2{u}_1(x,t)}{\partial x^2}=h\left(a\frac{{\partial }^2{u}_0(x,t)}{\partial x^2}-D^{\alpha }{u}_0(x,t)\right),$$

(26)

and for $m>1$:

$$\frac{{\partial }^2{u}_m(x,t)}{\partial x^2}=\left(1+ah\right)\frac{{\partial }^2{u}_{m-1}(x,t)}{\partial x^2}-h{D}^{\alpha }{u}_{m-1}(x,t).$$

(27)

To make the solution unambiguous, we need to supplement the above partial differential equations with additional conditions. For this we will use the conditions () and (). First equation we supplement with the conditions of form:

$$\begin{array}{cc}\hfill u_0& (\xi \left(t\right),t)+u_1(\xi \left(t\right),t)=u^{*},\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill -k& \frac{\partial (u_0+u_1)(\xi \left(t\right),t)}{\partial x}=\kappa D^{\alpha }\xi \left(t\right).\hfill \end{array}$$

(29)

For the rest of the equations ($m>1$), on the other hand, we set the conditions of form:

$$\begin{array}{cc}\hfill u_m(\xi \left(t\right),t)& =0,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill -k\frac{\partial u_m(\xi \left(t\right),t)}{\partial x}& =0.\hfill \end{array}$$

(31)

The above conditions provide us with the fact that any approximate solution constructed as the partial sum of the series (13) will meet the specified boundary conditions. As the initial approximation we can take the function that determines the initial condition:

$$u_0(x,t)=\phi \left(x\right).$$

(32)

Therefore, the problem was reduced to solving a series of partial differential Equations (26) and (27), with the conditions (28) and (29) and (30) and (31) respectively. The obtained equations are easier to solve than the initial problem. By knowing the function u or its approximation ${\widehat{u}}_n$, we can easily determine the missing boundary conditions:

$$\begin{array}{cc}\hfill \theta \left(t\right)& =u(0,t),\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill q\left(t\right)& =-k\frac{\partial u(x,t)}{\partial x}{|}_{x=0},\hfill \end{array}$$

(34)

or their approximations:

$$\begin{array}{cc}\hfill {\widehat{\theta }}_n\left(t\right)& ={\widehat{u}}_n(0,t),\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\widehat{q}}_n\left(t\right)& =-k\frac{\partial {\widehat{u}}_n(x,t)}{\partial x}{|}_{x=0}.\hfill \end{array}$$

(36)

It can be shown that the sum of the constructed series is a solution of the discussed equation.

Theorem 1.

Let the functions $u_m$, $m\ge 1$ be the solutions of the Equations (26) and (27) with the conditions (28) and () and (30) and () respectively. Then if the series ${\sum }_{m=0}^{\infty }{u}_m(x,t)$ is convergent in the domain Ω, then its sum designates the solution of the Equation (4).

The next theorem gives the sufficient condition for the convergence of the series constructed in the method.

Theorem 2.

Let the functions $u_m$, $m\ge 1$ be the solutions of the Equations (26) and (27) with the conditions (28) and (29) and (30) and (31) respectively. Then if the parameter h is selected in such a way that there exist such constants ${\gamma }_h\in (0,1)$ and $m_0\in \mathbb{N}$ that for each $m⩾m_0$ the inequality

$$\parallel u_{m+1}\parallel ⩽{\gamma }_h\parallel u_m\parallel$$

(37)

is satisfied in the domain Ω, then the series ${\sum }_{m=0}^{\infty }{u}_m(x,t)$ is uniformly convergent in Ω.

The last theorem gives the approximate solution error estimation, which we obtain using the partial sum of the series.

Theorem 3.

If assumptions of Theorem 2 are satisfied and additionally $n\in \mathbb{N}$ and $n⩾m_0$, then the estimation of error of the approximate solution is described by the following formula:

$$\parallel u-{\widehat{u}}_n\parallel ⩽\frac{{\gamma }_h^{n+1-m_0}}{1-{\gamma }_h}\parallel u_{m_0}\parallel ,$$

(38)

where ${\widehat{u}}_n$ is determined by the Equation (14).

The proofs for all of the above theorems are similar to those for the classic Stefan inverse problem and may be found in the paper [27].

4. Example

We will now present the method outlined in the previous chapter using an example. We will consider the problem for which under discussion we take the following data: $\overline{a}=1$, $\alpha =\frac{1}{2}$, $t^{*}=1$, $u^{*}=273$, $k=2$, $\overline{\kappa }=2$, $\phi \left(x\right)=274-x$, $\xi \left(t\right)=t+1$. Let us mention that all of the calculations presented here were carried out with the aid of Mathematica software [35].

Assuming that $u_0(x,t)=\phi \left(x\right)=274-x$ in the first step, we get:

$$u_1(x,t)=\frac{\sqrt{t}(1+t-x)}{\sqrt{\pi }}+x-1,$$

and

$$u_2(x,t)=-\frac{1}{24}h\left(5t-2x+2\right){\left(t-x+1\right)}^2.$$

The optimal value of the convergence control parameter equals to $-1.08347$. Figure 2 displays the graph of logarithm of the squared residual for $n=5$. This way for this value of h parameter we obtain the following approximate solutions:

$${\widehat{u}}_1(x,t)=\frac{\sqrt{t}(1+t-x)}{\sqrt{\pi }}+273$$

and

$${\widehat{u}}_2(x,t)=\frac{\sqrt{t}(1+t-x)}{\sqrt{\pi }}+0.04514458\left(5t-2x+2\right){\left(t-x+1\right)}^2+273.$$

Figure 3 displays the reconstructed boundary conditions, which are the temperature distribution and the heat flux on the boundary ${\Gamma }_1$. The obtained approximate temperature distribution for the whole $\mathsf{\Omega }$ area whereas is presented in the Figure 4. The next Figure 5, however, presents the convergence of the sequence of the partial sums which are the consecutive approximate solutions. In the picture the graphs of functions $R_i\left(t\right)=|{\widehat{\theta }}_i\left(t\right)-{\widehat{\theta }}_{i-1}\left(t\right)|$ and $Q_i\left(t\right)=|{\widehat{q}}_i\left(t\right)-{\widehat{q}}_{i-1}\left(t\right)|$ are presented using a logarithmic scale for $i=4,6,8,10$. The functions ${\widehat{\theta }}_i$ and ${\widehat{q}}_i$ determine the approximate temperature distribution ${\widehat{\theta }}_i$ and the heat flux ${\widehat{q}}_i$ on the boundary ${\Gamma }_1$, and are defined by Formulas (35) and (36) respectively. The logarithmic scale was used because all of those functions take very small values. At points where a peak down appears on the graph the corresponding function for which the absolute value is taken changes its sign. This way the absolute value equals to zero and its logarithm tends to minus infinity. Graphs of the $R_6$ and $R_{10}$ functions without a logarithmic scale are presented in Figure 6. Figure 5 shows as well that for the obtained approximate solutions the inequality (37) holds. We may assume that the next approximations will behave the same way. This way according to Theorem 2 the series forming the approximate solution is convergent, so its sum is the sought solution (see Theorem 1).

In the considered inverse problem, the information compensating for the lack of part of the input data is the known position of the moving interface. Therefore, in the example the stability of the procedure was also examined by terms of the errors of this position determination. The location of the moving interface was given with the random 1%, 2%, 5% and 10% errors. The relative errors of the boundary conditions reconstruction on the boundary ${\Gamma }_1$ and the temperature distribution in the area $\mathsf{\Omega }$ determined by the norm $L_2$ are displayed in the

Table 1. The relative errors of reconstruction the boundary conditions and temperature distribution for the disturbed input data.

Noise	${\widehat{\mathit{\theta }}}_{10}[\%]$	${\widehat{\mathit{q}}}_{10}[\%]$	${\widehat{\mathit{u}}}_{10}[\%]$
$\phantom{1}1\%$	$0.76$	$0.82$	$0.74$
$\phantom{1}2\%$	$2.62$	$2.75$	$2.61$
$\phantom{1}5\%$	$9.73$	$10.19$	$9.66$
$10\%$	$14.23$	$14.90$	$14.11$

. The relative errors were calculated for the shifted temperature ($u(x,t)-273$). In this case they achieve the greatest values then. For the output temperature the relative errors do not exceed the $0.08$% (for the greatest disturbance). In case of the heat flux, the shift does not affect the relative error. In the case when the input data gets disturbed by the $1\%$ and $2\%$ errors, the errors of reconstruction are on a similar level. Only for the bigger disturbances of the input data the errors of reconstruction increase by about $5\%$ relative to the value of the input data disturbance. These errors are not yet very drastic, especially if looking at the course of the reconstructed curves representing the boundary conditions (see Figure 7). As can be seen, the shape of the proper curve has been maintained for each case. The graphs of the absolute errors of the boundary conditions reconstruction for various disturbances of the input data are displayed in the Figure 8. The Figure 9 just like before presents the convergence of the approximate solutions in case of the input data disturbed by the $10\%$ error.

The Figure 10 presents the reconstructed boundary conditions for the different values of the order of the Caputo derivative $\alpha$. For $\alpha \in (0,\frac{1}{2}]$ this method converges quickly. The obtained solutions get bigger values with increasing value of the $\alpha$ parameter. On the other hand, for $\alpha >\frac{1}{2}$ the improper integrals that occur during the calculations do not converge. Thus, the method cannot be used in this case.

5. Conclusions

We use the homotopy analysis method to solve the one-phase fractional inverse Stefan problem. The mathematical model considered in the paper consists of a differential equation with a Caputo fractional derivative and the conditions on the boundaries of the considered domain, with the temperature distribution and heat flux unknown on one of the boundaries. In order to find these unknown functions, in accordance with the homotopy analysis method, an appropriate series is constructed which, if it converges, is the solution of the considered equation. We present theorems concerning the series convergence as well as the estimation of the error of the approximate solution. We also present a numerical example illustrating the application of the described method, its convergence and stability to the disturbances of the input data. As can be seen in the figures and the table presented in Section 4, the method converges quickly and the errors of the approximate solution are relatively small, especially for the input data disturbed by the 1% or 2% error. An important advantage of the method is the fact that the solution is obtained in the form of a continuous function and the considered domain does not need to be discretized. In the future, we plan to apply the method to two-phase and multidimensional problems, as well as to models with other types of fractional order derivative.