A Semi-explicit Algorithm for Parameters Estimation in a Time-Fractional Dual-Phase-Lag Heat Conduction Model

1. Introduction

Fourier’s law is a fundamental phenomenological relation in the heat transfer theory that describes the heat conduction in solids and fluids. However, it is not applicable for modelling the heat transfer in ultrafast processes (e.g. during laser heating and cooling [1]), on micro/nanoscales (e.g. heating of carbon nanotubes [2]), in some complex media (e.g. heat conduction in biological tissues [3] and materials with non-homogeneous inner structure [4]) . For this reason, various linear and nonlinear generalizations of Fourier’s law have been proposed by many researchers through the past two centuries (see books [5,6,7] and reviews [8,9,10,11]).

The thermal conductivity is a unique property of a material in Fourier’s law whereas non-Fourier models usually include several material thermal properties. For example, the frequently used dual-phase-lag (DPL) heat conduction model, proposed by Tzou [12], incorporates the thermal conductivity, the thermal relaxation time and the thermal retardation. As a result, the possibility of using a non-Fourier model in real-world applications is based on an ability to obtain necessary thermal properties of a material. From the mathematical point of view, the problem of parameters estimation can be considered as an inverse coefficient problem which is usually ill-posed problem. Thus, the development of theoretical techniques for estimation of thermal properties in non-Fourier models is a challenging problem of mathematical modelling in the heat transfer theory.

The method of time integral characteristics (TIC) is an efficient technique for estimation of constant parameters in linear models. It was proposed by Shatalov [13] at the end of last century for solving inverse coefficient problems of heat conduction theory. The method is based on integral transformation of the initial-boundary value problem for the considered linear model on the time variable and solving the corresponding inverse coefficient problem in the transformed space. The absence of necessity to perform the inverse integral transform is a main advantage of this method. Later, this method was extended to time-fractional diffusion models [14,15].

In this study, we focus on solving the problem of parameters estimation in a time-fractional generalization of DPL heat conduction model by the method of TIC. Such model allow us to take into account power-law memory thermal effects in a material by incorporating time derivatives of fractional orders (see, e.g., [16,17]) into the thermal constitutive relation. Thus, in this model the orders of fractional differentiation are additional parameters that also have to be estimated.

Time-fractional dual-phase-lag (TFDPL) model was proposed by Xu and Jiang [18] to interpret the experiment results for processed meat. The Caputo time-fractional derivatives of two different orders are used in this model. The authors obtained analytical solution of the corresponding bioheat transfer equation and solved the inverse problem for estimation of model parameters by applying the nonlinear least-square method. The same model was used in [19] for treating thermoelastic response of skin subjected to sudden temperature shock. In [20], a fundamental solution of the TFDPL heat conduction problem was obtained. Also, a TFDPL model with a single order of fractional differentiation was considered by several scholars. In [21], such a model was used for describing the heat conduction in a multi-layered spherical medium with azimuthal symmetry. In [22], a similar model with temperature jump boundary condition was utilized for numerical simulation of heat transfer in transistors. Numerical schemes for solving several TFDPL heat conduction problems was proposed in [23,24].

However, it is necessary to note that in all papers mentioned above the considered TFDPL models have been obtained from the classical DPL model by formal replacing of the integer order time derivatives by their fractional analogues. In this paper we overcome this weakness by proposing the derivation of TFDPL model from a general linear constitutive relation for the heat transfer by conduction.

The paper is organized as follows. Section 2 contains a brief description of the time integral characteristic method. Section 3 is devoted to derivation of the TFDPL heat conduction model and corresponding non-Fourier heat conduction equation. A proposed semi-explicit algorithm for TFDPL model parameters estimation is described in Section 4. An illustrative example of using this algorithm is presented in Section 5.

2. Preliminaries

This section gives a brief description of the TIC method. This method was proposed for solving inverse problems of parameters estimation in linear evolution equations.

Let us illustrate the basic idea of the TIC method by a simple problem of the thermal diffusivity estimation. We consider the heat equation

$${\displaystyle \frac{\partial T}{\partial t}}=a{\displaystyle \frac{{\partial }^{2}T}{\partial x^2}},t>0,0<x<\infty .$$

(1)

Here x and t are spatial and temporal variables, respectively, $T(x,t)$ is the temperature field, a is the thermal diffusivity. This equation is accompanied with the initial condition

$$T(x,0)=0,0<x<\infty ,$$

(2)

the boundary conditions

$$T(0,t)=T_0\left(t\right),\underset{x\to \infty }{lim}T(x,t)=0,t>0,$$

(3)

and the additional internal condition

$$T(l,t)=T_l\left(t\right),t>0.$$

(4)

Here $T_0\left(t\right)$, $T_l\left(t\right)$ are known functions. Then the inverse coefficient problem is stated as follows: given the initial boundary value problem (1), (2), (3) and the additional condition (4), find the constant thermal diffusivity a.

The method of TIC is based on an integral transformation of the temperature field with respect to time. The Laplace transform

$$T^{*}(x,p)={\int }_0^{\infty }{e}^{-pt}T(x,t)dt$$

(5)

can be efficiently used for this purpose. The function $T_L^{*}\left(p\right)=T^{*}(L,p)$ is referred to as a time integral characteristic of the temperature field T at the point $x=L$.

The initial boundary value problem (1), (2), (3) after Laplace transformation takes the form

$$a{T}^{{*}^{\prime \prime }}(x,p)-p{T}^{*}(x,p)=0,0<x<\infty ,$$

(6)

$$T^{*}(0,p)=T_0^{*}\left(p\right),\underset{x\to \infty }{lim}{T}^{*}(x,p)=0,$$

(7)

and (4) gives

$$T^{*}(l,p)=T_l^{*}\left(p\right).$$

(8)

In (6), prime denotes differentiation with respect to x.

The solution of (6), (7) is

$$T^{*}(x,p)=T_0^{*}\left(p\right)e^{-\lambda x},\lambda =\sqrt{{\displaystyle \frac{p}{a}}}.$$

Then, by using (8), we find the following explicit representation of the thermal diffusivity via TICs of the temperature field:

$$a={\displaystyle \frac{p{l}^2}{{ln}^2\phi \left(p\right)}},\phi \left(p\right)={\displaystyle \frac{T_l^{*}\left(p\right)}{T_0^{*}\left(p\right)}}.$$

(9)

A main advantage of the described technique is that there is no necessity to perform the inverse Laplace transform. If the temperature functions $T_0\left(t\right)$ and $T_l\left(t\right)$ are known exactly, the representation (9) gives the exact value of a for all permitted values of p. However, in practice such functions are usually measured in an experiment with some errors. Then the Laplace parameter p should be considered as a regularization parameter and its value should be chosen in agreement with experimental errors. The explicit representation (9) permits to obtain a linear estimate of relative error for the thermal diffusivity as a function of p (see [13,14,15] for more details). The minimum of this function gives the optimal value of p.

3. The time-fractional dual-phase-lag heat equation

A TFDPL heat conduction model can be obtained similarly to the time-fractional Zener model in the theory of linear fractional viscoelasticity [25].

A general linear constitutive relation for the heat transfer by conduction in homogeneous and isotropic materials with memory is defined mathematically using a Riemann–Stieltjes integral as

$$\mathbf{q}\left(t\right)={\int }_{-\infty }^{\infty }r(t-s)d\mathbf{g}\left(s\right).$$

(10)

Here $\mathbf{q}\left(t\right)\equiv \mathbf{q}(\mathbf{x},t)$ is the heat flux vector, $\mathbf{g}\left(t\right)\equiv \mathbf{g}(\mathbf{x},t)=-\nabla T$ is the temperature gradient, and $r\left(t\right)$ is the heat conduction relaxation function which does not depend on the spatial coordinate $\mathbf{x}$.

In accordance to the physical principle of causality, the relaxation function $r\left(t\right)$ is zero for negative time. Hence, the constitutive relation (10) takes the form

$$\mathbf{q}\left(t\right)={\int }_{-\infty }^{t}r(t-s)d\mathbf{g}\left(s\right).$$

(11)

It is easy to see that this equation reduces to Fourier’s law $\mathbf{q}=-k\nabla T$ if $\mathbf{g}(-\infty )=0$ and $r\left(t\right)=k$ where k is the thermal conductivity which is a constant in time.

Assuming that $r\left(t\right)$ and $g^{\prime }\left(t\right)$ are continuous functions, we can rewrite (11) in a more convenient form

$$\mathbf{q}\left(t\right)={\int }_{-\infty }^{t}r(t-s){\mathbf{g}}^{\prime }\left(s\right)ds.$$

(12)

Let us now consider the case when the heat conduction relaxation function has the form

$$r\left(t\right)=r_0+r_1{E}_{\alpha }\left(-{\displaystyle \frac{t^{\alpha }}{{\tau }_q}}\right),\alpha \in (0,1),$$

(13)

where $r_0,r_1,{\tau }_q$ are constants and

$$E_{\alpha }\left(z\right)=\sum _{n=0}^{\infty }{\displaystyle \frac{z^n}{\Gamma (\alpha n+1)}}$$

(14)

is the Mittag-Leffler function (see, e.g., [17]). This function has the known property

$${ }_a^C{D}_t^{\alpha }{E}_{\alpha }\left[\lambda {(t-a)}^{\alpha }\right]=\lambda E_{\alpha }\left[\lambda {(t-a)}^{\alpha }\right],a\in \mathbb{R},\lambda \in \mathbb{C},$$

(15)

i.e. it is invariant with respect to the left-sided Caputo fractional derivative of order $\alpha$. This fractional derivative reads

$$\left({ }_a^C{D}_t^{\alpha }y\right)\left(t\right)\equiv \left({ }_a^{ }{I}_t^{1-\alpha }{y}^{\prime }\right)\left(t\right)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_a^t{\displaystyle \frac{y^{\prime }\left(s\right)}{{(t-s)}^{\alpha }}}ds,t>a,\alpha \in (0,1).$$

(16)

Here ${ }_a^{ }{I}_t^{1-\alpha }$ is the left-sided fractional integral operator of order $1-\alpha$. Letting $a=-\infty$ in (16), we obtain the Caputo fractional operator ${ }_{-\infty }^C{D}_t^{\alpha }$ for the whole axis $\mathbb{R}$. Applying this operator to both sides of (12), we get

$$\begin{array}{cc}\left({ }_{-\infty }^C{D}_t^{\alpha }\mathbf{q}\right)\left(t\right)\hfill & {\displaystyle ={ }_{-\infty }^C{D}_t^{\alpha }\left({\int }_{-\infty }^{t}r(t-s){\mathbf{g}}^{\prime }\left(s\right)ds\right)}\hfill \\ & {\displaystyle =r\left(0\right)\left({ }_{-\infty }^C{D}_t^{\alpha }\mathbf{g}\right)\left(t\right)+{ }_{-\infty }^{ }{I}_t^{1-\alpha }\left({\int }_{-\infty }^t{r}^{\prime }(t-s){\mathbf{g}}^{\prime }\left(s\right)ds\right)}\hfill \\ & {\displaystyle =r\left(0\right)\left({ }_{-\infty }^C{D}_t^{\alpha }\mathbf{g}\right)\left(t\right)+{\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_{-\infty }^t{\displaystyle \frac{1}{{(t-\xi )}^{\alpha }}}\left({\int }_{-\infty }^{\xi }{r}^{\prime }(\xi -s){\mathbf{g}}^{\prime }\left(s\right)ds\right)d\xi .}\hfill \end{array}$$

We assume here that $\mathbf{g}\left(t\right)$ is a <<sufficiently good>> function so that all integrals exist. Changing the order of integration in the last term of the above expression, we obtain

$$\begin{array}{c}\left({ }_{-\infty }^C{D}_t^{\alpha }\mathbf{q}\right)\left(t\right)=r\left(0\right)\left({ }_{-\infty }^C{D}_t^{\alpha }\mathbf{g}\right)\left(t\right)+{\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_{-\infty }^t{\mathbf{g}}^{\prime }\left(s\right)\left({\int }_s^t{\displaystyle \frac{r^{\prime }(\xi -s)}{{(t-\xi )}^{\alpha }}}d\xi \right)ds\hfill \\ \hfill =r\left(0\right)\left({ }_{-\infty }^C{D}_t^{\alpha }\mathbf{g}\right)\left(t\right)+{\int }_{-\infty }^t\left({ }_s^C{D}_t^{\alpha }r\right)(t-s){\mathbf{g}}^{\prime }\left(s\right)ds.\end{array}$$

(17)

Substituting $r\left(t\right)$ given by (13) into (17), and using (15), we have

$$\left({ }_{-\infty }^C{D}_t^{\alpha }\mathbf{q}\right)\left(t\right)=(r_0+r_1)\left({ }_{-\infty }^C{D}_t^{\alpha }\mathbf{g}\right)\left(t\right)-{\displaystyle \frac{r_1}{{\tau }_q}}{\int }_{-\infty }^t{E}_{\alpha }\left(-{\displaystyle \frac{{(t-s)}^{\alpha }}{{\tau }_q}}\right){\mathbf{g}}^{\prime }\left(s\right)ds.$$

(18)

On the other hand, the relation (12) with (13) takes the form

$$\mathbf{q}\left(t\right)=r_0\mathbf{g}\left(t\right)+r_1{\int }_{-\infty }^t{E}_{\alpha }\left(-{\displaystyle \frac{{(t-s)}^{\alpha }}{{\tau }_q}}\right){\mathbf{g}}^{\prime }\left(s\right)ds.$$

(19)

It is easy to see that integrals in the last terms of (18) and (19) coincide and therefore can be excluded. As a result, we obtain the time-fractional constitutive relation

$$\mathbf{q}+{\tau }_q{ }_{-\infty }^C{D}_t^{\alpha }\mathbf{q}=r_0\mathbf{g}+{\tau }_q(r_0+r_1)\left({ }_{-\infty }^C{D}_t^{\alpha }\mathbf{g}\right).$$

Using the definition of the function $\mathbf{g}$, this relation can be written as

$$\mathbf{q}+{\tau }_q{ }_{-\infty }^C{D}_t^{\alpha }\mathbf{q}=-k\left[\nabla T+{\tau }_T{ }_{-\infty }^C{D}_t^{\alpha }(\nabla T)\right].$$

(20)

Here $k=r_0$ is the thermal conductivity, ${\tau }_q$ is the fractional analogue of the thermal relaxation time (the so-called phase lag in the heat flux), and ${\tau }_T={\tau }_q(1+r_1/r_0)$ is the fractional analogue of the thermal retardation (the so-called temperature gradient phase lag).

The constitutive relation (20) describes the heat conduction in a material with full power-law memory. This relation is invariant with respect to translation in time and therefore the time origin can be arbitrarily chosen.

Let us now assume that there is not heat transfer in a material for time $t<0$. Then the relation (20) reduces to

$$\mathbf{q}+{\tau }_q{ }_0^C{D}_t^{\alpha }\mathbf{q}=-k\left[\nabla T+{\tau }_T{ }_0^C{D}_t^{\alpha }(\nabla T)\right].$$

(21)

Note that this relation is not invariant with respect to translation in time and the time origin is fixed in this case.

Now we can obtain TFDPL heat equation. The energy conservation law for a constant property material without heat sources can be written as

$$c\rho {\displaystyle \frac{\partial T}{\partial t}}+\nabla ·\mathbf{q}=0,$$

(22)

where c is the specific heat and $\rho$ is the density of material. Combining (21) and (22), after simple algebra, we get

$${\displaystyle \frac{\partial T}{\partial t}}+{\tau }_q{ }_0^C{D}_t^{\alpha +1}T=a\left[\Delta T+{\tau }_T{ }_0^C{D}_t^{\alpha }\left(\Delta T\right)\right],\alpha \in (0,1),$$

(23)

where $a=k/\left(c\rho \right)$ is the thermal diffusivity. TFDPL heat equation (22) models heat conduction in a material with power-law memory and constant temperature field for time $t<0$.

4. An Algorithm of TFDPL Model Parameters Estimation

Let us consider a one-dimensional case of TFDPL heat equation (23) in a half space, namely

$${\displaystyle \frac{\partial T}{\partial t}}+{\tau }_q{ }_0^C{D}_t^{\alpha +1}T=a\left[{\displaystyle \frac{{\partial }^{2}T}{\partial x^2}}+{\tau }_T{ }_0^C{D}_t^{\alpha }\left({\displaystyle \frac{{\partial }^{2}T}{\partial x^2}}\right)\right],t>0,0<x<\infty ,0<\alpha <1.$$

(24)

This is a time-fractional equation of order $\alpha +1\in (1,2)$ and therefore two initial conditions are needed for its unique solvability. We take them in the form

$$T(x,0)=T^0,{\displaystyle \frac{\partial T(x,t)}{\partial t}}{|}_{t=0}=0,0<x<\infty ,$$

(25)

where $T^0$ is a constant initial temperature.

We will also assume that (24) is accompanied by the boundary conditions

$$T(0,t)=T_0\left(t\right),\underset{x\to \infty }{lim}T(x,t)=T^0,t>0,$$

(26)

and by the additional internal condition

$$T(l,t)=T_l\left(t\right),l>0,t>0.$$

(27)

Here $T_0\left(t\right)$, $T_l\left(t\right)$ are known functions.

We will consider the following inverse problem: given the initial boundary value problem (24), (25), (26) and the additional condition (27), find the constants a, ${\tau }_T$, ${\tau }_q$, $\alpha$. For solving this problem, the method of TIC can be efficiently used.

For convenience, we introduce a new function $\theta (x,t)=T(x,t)-T^0$. Then the problem (24) –(27) takes the form

$${\displaystyle \frac{\partial \theta }{\partial t}}+{\tau }_q{ }_0^C{D}_t^{\alpha +1}\theta =a\left[{\displaystyle \frac{{\partial }^2\theta }{\partial x^2}}+{{\tau }_T}_0^C{D}_t^{\alpha }\left({\displaystyle \frac{{\partial }^2\theta }{\partial x^2}}\right)\right],t>0,0<x<\infty ,0<\alpha <1,$$

(28)

$$\theta (x,0)=0,{\displaystyle \frac{\partial \theta (x,t)}{\partial t}}{|}_{t=0}=0,0<x<\infty ,$$

(29)

$$\theta (0,t)={\theta }_0\left(t\right),\underset{x\to \infty }{lim}\theta (x,t)=0,t>0,$$

(30)

$$\theta (l,t)={\theta }_l\left(t\right),l>0,t>0.$$

(31)

Here the functions ${\theta }_0\left(t\right)=T_0\left(t\right)-T^0$ and ${\theta }_l\left(t\right)=T_l\left(t\right)-T^0$ are known.

The initial-boundary problem (28), (29), (30) after Laplace transform can be written as

$$a(1+{\tau }_T{p}^{\alpha }){\theta }^{{*}^{\prime \prime }}(x,p)-p(1+{\tau }_q{p}^{\alpha }){\theta }^{*}(x,p)=0,0<x<\infty ,0<\alpha <1,$$

(32)

$${\theta }^{*}(0,t)={\theta }_0^{*}\left(p\right),\underset{x\to \infty }{lim}{\theta }^{*}(x,p)=0,$$

(33)

and (31) gives

$${\theta }^{*}(l,p)={\theta }_l^{*}\left(p\right),l>0.$$

(34)

Here ${\theta }^{*}(x,p)$ denotes the Laplace transform of $\theta (x,t)$ which is defined by

$${\theta }^{*}(x,p)={\int }_0^{\infty }{e}^{-pt}\theta (x,t)dt.$$

In (32), prime denotes differentiation with respect to x.

The solution of (32), (33) is

$${\theta }^{*}(x,p)={\theta }_0^{*}\left(p\right)e^{-\lambda x},\lambda =\sqrt{{\displaystyle \frac{p(1+{\tau }_q{p}^{\alpha })}{a(1+{\tau }_T{p}^{\alpha })}}}.$$

Using the additional condition (34), we get the main equation for parameters estimation which can be written as

$$\Phi (a,b,{\tau }_q,\alpha ;p)=0,$$

(35)

where

$$\Phi (a,b,{\tau }_q,\alpha ;p)=(a+b{p}^{\alpha })\psi \left(p\right)-{\tau }_q{p}^{\alpha }-1,b=a{\tau }_t,\psi \left(p\right)={\displaystyle \frac{1}{p{l}^2}}{ln}^2{\displaystyle \frac{{\theta }_0^{*}\left(p\right)}{{\theta }_l^{*}\left(p\right)}}.$$

(36)

Note that the function $\Phi$ is linear with respect to a, b, ${\tau }_q$, and nonlinear with respect to $\alpha$.

As it was mentioned in Preliminaries section, the problem of finding the Laplace parameter p arises if the functions $T_0\left(t\right)$ and $T_l\left(t\right)$ are not known exactly. In the method of TIC, the Laplace parameter p is assumed to be real and positive. Therefore, it is naturally to assume that this parameter belongs to a finite interval $[p_{min},p_{max}]$$(0<p_{min}<p_{max}<\infty )$. A discussion of different approaches to estimation of $p_{min}$ and $p_{max}$ can be found in [14,15]. Then the considered inverse problem can be reduced to minimization problem

$$f(a,b,{\tau }_q,\alpha )\equiv {\int }_{p_{min}}^{p_{max}}{\Phi }^2(a,b,{\tau }_q,\alpha ;p)\to min.$$

(37)

This is a classical problem of finding a minimum of four variables function f. The physical constraints are $a>0$, $b>0$, ${\tau }_q>0$, and $\alpha \in (0,1)$. In general, this problem can be solved numerically by using different optimization software.

However, explicit TIC representations for desired parameters can be obtained in a special case when the order of fractional differentiation $\alpha$ is known. Let us consider (37) as unconstrained minimization problem. It is obvious that the function f defined by (37) is a quadratic function in three variables a, b and ${\tau }_q$. The necessary conditions for local optimality reads

$${\displaystyle \frac{\partial f}{\partial a}}=0,{\displaystyle \frac{\partial f}{\partial b}}=0,{\displaystyle \frac{\partial f}{\partial {\tau }_q}}=0.$$

These conditions give the system of linear equations

$$Az=B,$$

(38)

where

$$A=\left(\begin{array}{ccc}{I}_2^0& I_2^1& -I_1^1\\ I_2^1& I_2^2& -I_1^2\\ -I_1^1& -I_1^2& I_0^2\end{array}\right),B=\left(\begin{array}{c}{I}_1^0\\ I_1^1\\ -I_0^1\end{array}\right),z=\left(\begin{array}{c}a\\ b\\ {\tau }_q\end{array}\right),$$

and

$$I_m^n={\int }_{p_{min}}^{p_{max}}{\psi }^m\left(p\right)p^{\alpha n}dp.$$

Note that the matrix A is symmetric.

Using Cramer’s rule, the solution of (38) can be written in the explicit form

$$a={\displaystyle \frac{{\Delta }_1}{\Delta }},b={\displaystyle \frac{{\Delta }_2}{\Delta }},{\tau }_q={\displaystyle \frac{{\Delta }_3}{\Delta }},$$

(39)

where

$$\Delta =det\left(A\right),{\Delta }_i=det\left(A_i\right),i=1,2,3,$$

and $A_i$ is the matrix formed by replacing the i-th column of A by the column vector B. Thus, the explicit representations (39) permit to obtain the values of a, b and ${\tau }_q$ for a given value of $\alpha$.

The representations (39) can also be used in the case of unknown $\alpha$. Then we have

$${\displaystyle \frac{\partial f}{\partial \alpha }}=0$$

and obtain

$$ab{J}_2^1+b^2{J}_2^2-2b{\tau }_q{J}_1^2-(a{\tau }_q+b)J_1^1+{\tau }_q^2{J}_0^2+{\tau }_q{J}_0^1=0,$$

(40)

where

$$J_m^n\equiv J_m^n\left(\alpha \right)={\int }_{p_{min}}^{p_{max}}{\psi }^m\left(p\right)p^{\alpha n}lnpdp.$$

In (40), the parameters a, b and ${\tau }_q$ are the functions of $\alpha$ which are defined by (39). We thus obtain a single equation for $\alpha$. The equation (40) is nonlinear and quite complex. Therefore, it should be solved numerically.

As a result, we can state the following semi-explicit algorithm for parameters estimation in TFDPL heat conduction model:

• The Laplace transforms ${\theta }_0^{*}\left(p\right)$ and ${\theta }_l^{*}\left(p\right)$ are computed for a given functions ${\theta }_0\left(t\right)$ and ${\theta }_l\left(t\right)$, respectively.
• The function $\psi \left(p\right)$ from (36) is found.
• The order of fractional differentiation $\alpha$ is obtained by numerical solving of the nonlinear equation (40).
• The desired parameters a, b and ${\tau }_q$ are calculated by the explicit TIC representations given by (39).

It is necessary to note that the proposed algorithm is based on the unconstrained minimization problem. As a result, the order of fractional differentiation $\alpha$ which is obtained as the solution of (40) is not necessary belong to the interval $[0,1]$. Then the constrained minimization problem mentioned above should be considered and solved numerically. Note that usually this situation arises when the initial functions ${\theta }_0\left(t\right)$ and ${\theta }_l\left(t\right)$ have quite large errors (usually more than 5 %).

The considered problem of parameters estimation belongs to the class of inverse coefficient problems. Hence, it is an ill-posed problem in most cases. In the proposed algorithm, the stabilization of solution is achieved by integration with respect to the Laplace parameter p. However, numerical experiments show that the solution is stable only if $p_{max}/p_{min}=O\left({10}^k\right)$ with $k\ge 1$. If $p_{max}/p_{min}=O\left(1\right)$, the determinant $det\left(A\right)$ is close to zero and corresponding approximate solution is unstable. An additional regularization is needed in this case. For example, the Tikhonov regularization method can be used for this purpose.

5. An Example

To illustrate the above algorithm, let us consider the equation (28) with

$$a=1,{\tau }_T=8,{\tau }_q=0.2.$$

Hence, $b=8$ and

$$\psi \left(p\right)={\displaystyle \frac{1+0.2p^{\alpha }}{1+8p^{\alpha }}}.$$

(41)

Three different values of fractional order are considered: $\alpha =0.75$, $0.5$ and $0.25$. The graphs of the function $\psi \left(p\right)$ defined by (41) for these values of $\alpha$ are shown in Figure 1.

Denote by $\overline{T}(x,t)$ and $\tilde{T}(x,t)$ the exact and perturbed temperature fields, respectively. Let

$$|\overline{T}(x,t)-\tilde{T}(x,t)|\le {\Delta }_T,$$

where ${\Delta }_T$ is the upper error bound. Then it is easy to prove (see, e.g., [13]) that

$$|{\overline{T}}^{*}(x,p)-{\tilde{T}}^{*}(x,p)|\le {\displaystyle \frac{{\Delta }_T}{p}}.$$

Thus, the error of $T^{*}(x,p)$ is increased as the parameter p is decreased. The same is valid for the function $\psi \left(p\right)$.

To simulate experimental errors, we approximated the function $\psi \left(p\right)$ on the interval $[p_{min},p_{max}]$ by polynomials of different degrees with respect to a new dependent variable

$$P=lnp.$$

The values of $p_{min}=0.01$ and $p_{max}=10$ were used in all numerical calculations. Also, we utilized the relative error as a metric of accuracy. For a quantity q, it is defined by

$$\delta q={\displaystyle \frac{|\overline{q}-\tilde{q}|}{\overline{q}}}100\%,$$

where $\overline{q}$ and $\tilde{q}$ are exact and perturbed values of q, respectively.

Case 1: $\alpha =0.75$.

The following approximations of $\psi \left(p\right)$ were constructed:

$$\begin{array}{ccc}{\psi }_1\hfill & =\hfill & -2.0871605·{10}^{-4}{P}^3+1.5908660·{10}^{-2}{P}^2-7.6804434·{10}^{-2}P+0.14207071,\hfill \\ {\psi }_2\hfill & =\hfill & -4.3025293·{10}^{-4}{P}^4-1.0286478·{10}^{-3}{P}^3+2.0448222·{10}^{-2}{P}^2\hfill \\ \hfill & \hfill & -7.4697380·{10}^{-2}P+.13433121,\hfill \\ {\psi }_3\hfill & =\hfill & 5.4973602·{10}^{-6}{P}^6+1.1001135·{10}^{-4}{P}^5-4.0590709·{10}^{-4}{P}^4\hfill \\ \hfill & \hfill & -2.4193687·{10}^{-3}{P}^3+2.1332786·{10}^{-2}{P}^2-7.2620046·{10}^{-2}P+0.13333333.\hfill \end{array}$$

(42)

The graphs of relative errors for these functions are plotted in Figure 2. It can be seen that the maximum value of relative error is approximately equal to $6.5\%$ for ${\psi }_1$, $2\%$ for ${\psi }_2$, and $0.5\%$ for ${\psi }_3$.

Table 1 contains the results of parameters estimation by using the explicit expressions given in (39) for $\alpha =0.75$ and approximate functions ${\psi }_i$ from (42). Note that in this case $detA=O\left({10}^{-1}\right)$. It can be seen that the relative errors of a and b are of the same order of magnitude as the corresponding relative errors of functions ${\psi }_i$. The relative errors of ${\tau }_t$ and ${\tau }_q$ are also in good agreement with the relative errors of ${\psi }_3$ and ${\psi }_2$. However, in the case of quite large initial error when the function ${\psi }_1$ is used, the error level of ${\tau }_t$ and ${\tau }_q$ is highly increased.

Table 2 contains the results of parameters estimation for unknown $\alpha$. In this case the equation (40) was solved for each ${\psi }_i$ ($i=1,2,3$). As can be seen from the table, the accuracy of $\alpha$ identification is highly depends on error of initial data. The same is valid for other parameters.

Case 2: $\alpha =0.5$.

The following approximations of $\psi \left(p\right)$ were obtained:

$$\begin{array}{ccc}{\psi }_1\hfill & =\hfill & 7.9024062·{10}^{-3}{P}^2-5.0950384·{10}^{-2}P+0.13779823,\hfill \\ {\psi }_2\hfill & =\hfill & -4.0453570·{10}^{-4}{P}^3+8.4778168·{10}^{-3}{P}^2-4.8909305·{10}^{-2}P+0.13482677,\hfill \\ {\psi }_3\hfill & =\hfill & -1.1327318·{10}^{-4}{P}^4-5.8706810·{10}^{-4}{P}^3+9.6778542·{10}^{-3}{P}^2\hfill \\ \hfill & \hfill & -4.8760918·{10}^{-2}P+0.13333333.\hfill \end{array}$$

(43)

The graphs of relative errors for these functions are plotted in Figure 3. The maximum value of relative error is approximately equal to $5\%$ for ${\psi }_1$, $2.2\%$ for ${\psi }_2$, and $0.5\%$ for ${\psi }_3$.

The results of parameters estimation are given in Table 3 for $\alpha =0.5$, and in Table 4 for initially unknown $\alpha$. In general, the results of this case are close to previous one. However, the magnitudes of relative errors for estimated parameters are greater than those obtained previously. This is because $detA=O\left({10}^{-3}\right)$ in this case.

Case 3: $\alpha =0.25$.

In this case, the following approximations of $\psi \left(p\right)$ were constructed:

$$\begin{array}{ccc}{\psi }_1\hfill & =\hfill & -2.7262704·{10}^{-2}P+0.14595998;\hfill \\ {\psi }_2\hfill & =\hfill & 2.3647866·{10}^{-3}{P}^2-2.4888663·{10}^{-2}P+0.13385186;\hfill \\ {\psi }_3\hfill & =\hfill & -8.7191573·{10}^{-5}{P}^3+2.3025763·{10}^{-3}{P}^2-2.4102390·{10}^{-2}P+0.133395.\hfill \end{array}$$

(44)

The graphs of relative errors for the functions in (44) are plotted in Figure 4. The maximum value of relative error is approximately equal to $10\%$ for ${\psi }_1$, $1\%$ for ${\psi }_2$, and $0.1\%$ for ${\psi }_3$.

Table 5 contains the results of parameters estimation by using (39) for $\alpha =0.25$ and approximations (44). In this case we found $detA=O\left({10}^{-5}\right)$ and therefore the estimation results are highly sensitive to errors of initial data. It follows from the table that if we use ${\psi }_1$ function having relative error at most $\pm 10\%$, the proposed algorithm does not permit to identify ${\tau }_q$ with the chosen values of $p_{min}$ and $p_{max}$. However, we obtained reasonable values of all desired parameters for the function ${\psi }_3$. This demonstrates the stability of the algorithm.

Finally, Table 2 contains table 6 citatin, table 6 is not cited. the results of parameters estimation for unknown $\alpha$. It can be seen from the table that in all cases we have a high level of error, especially for ${\tau }_q$.

Thus, we can conclude that the proposed algorithm permits to estimate the parameters of TFDPL model with a reasonable accuracy if $detA=O\left({10}^k\right)$ with $k\ge -3$. For fixed values of $p_{min}$ and $p_{max}$, decreasing of $\alpha$ leads to decreasing of $detA$. Hence, the values of $p_{min}$ and $p_{max}$ should be chosen dependently on $\alpha$. Also, it should be pointed out that the parameters $\alpha$ and ${\tau }_q$ are more sensitive to error level in initial data than a and ${\tau }_T$.