Analytical and numerical treatment of the heat conduction equation obtained via time-fractional distributed-order heat conduction law

TL;DR

This paper develops analytical and numerical methods to solve a generalized heat conduction equation incorporating fractional and distributed-order laws, demonstrating their effectiveness through numerical examples.

Contribution

It introduces a combined analytical and numerical framework for solving distributed-order fractional heat conduction equations, extending classical models.

Findings

Analytical solutions obtained via Fourier and Laplace transforms.

Numerical solutions using finite difference schemes show good agreement with analytical results.

The methods effectively handle multi-term and power-type distributed-order laws.

Abstract

Generalization of the heat conduction equation is obtained by considering the system of equations consisting of the energy balance equation and fractional-order constitutive heat conduction law, assumed in the form of the distributed-order Cattaneo type. The Cauchy problem for system of energy balance equation and constitutive heat conduction law is treated analytically through Fourier and Laplace integral transform methods, as well as numerically by the method of finite differences through Adams-Bashforth and Gr\"{u}nwald-Letnikov schemes for approximation derivatives in temporal domain and leap frog scheme for spatial derivatives. Numerical examples, showing time evolution of temperature and heat flux spatial profiles, demonstrate applicability and good agreement of both methods in cases of multi-term and power-type distributed-order heat conduction laws.

Tables

Table 1. Table 1: Errors of numerically obtained solutions for temperature using Adams-Bashforth, centered with RAW filter, centered and Euler approximation schemes, with respect to analytically obtained response to Gaussian function in case of multi-term heat conduction law.

	$t_i$	Adams- Bashworth	Centered with RAW filter	Centered	Euler
${\delta }_{l^2}T\left(t_i\right)$	$0.01$	$2.760\times {10}^{-5}$	$2.768\times {10}^{-5}$	$2.777\times {10}^{-5}$	$1.829\times {10}^{-2}$
	$0.015$	$1.655\times {10}^{-5}$	$1.661\times {10}^{-5}$	$1.666\times {10}^{-5}$	$4.668\times {10}^{-2}$
	$0.02$	$1.108\times {10}^{-5}$	$1.113\times {10}^{-5}$	$1.115\times {10}^{-5}$	$7.782\times {10}^{-2}$
	$0.035$	$5.732\times {10}^{-6}$	$5.754\times {10}^{-6}$	$5.793\times {10}^{-6}$	$0.150$
	$0.05$	$5.167\times {10}^{-6}$	$5.174\times {10}^{-6}$	$1.020\times {10}^{-3}$	$0.203$
	$0.065$	$5.514\times {10}^{-6}$	$5.525\times {10}^{-6}$	$1.830\times {10}^6$	$0.240$
${\Delta }_{l^2}T\left(t_i\right)$	$0.01$	$2.089\times {10}^{-7}$	$2.094\times {10}^{-7}$	$2.101\times {10}^{-7}$	$1.384\times {10}^{-4}$
	$0.015$	$1.092\times {10}^{-7}$	$1.096\times {10}^{-7}$	$1.099\times {10}^{-7}$	$3.079\times {10}^{-4}$
	$0.02$	$6.401\times {10}^{-8}$	$6.429\times {10}^{-8}$	$6.442\times {10}^{-8}$	$4.495\times {10}^{-4}$
	$0.035$	$2.442\times {10}^{-8}$	$2.452\times {10}^{-8}$	$2.468\times {10}^{-8}$	$6.416\times {10}^{-4}$
	$0.05$	$1.797\times {10}^{-8}$	$1.800\times {10}^{-8}$	$3.547\times {10}^{-6}$	$7.051\times {10}^{-4}$
	$0.065$	$1.639\times {10}^{-8}$	$1.643\times {10}^{-8}$	$5.439\times {10}^3$	$7.144\times {10}^{-4}$
${\Delta }_{l^{\infty }}T$		$5.536\times {10}^{-5}$	$5.552\times {10}^{-5}$	$15.69$	$2.456\times {10}^{-3}$

Table 2. Table 2: Errors of numerically obtained solutions for heat flux using Adams-Bashforth, centered with RAW filter, centered and Euler approximation schemes, with respect to analytically obtained response to Gaussian function in case of multi-term heat conduction law.

	$t_i$	Adams- Bashworth	Centered with RAW filter	Centered	Euler
${\delta }_{l^2}q\left(t_i\right)$	$0.01$	$3.842\times {10}^{-6}$	$3.894\times {10}^{-6}$	$3.830\times {10}^{-6}$	$2.247\times {10}^{-2}$
	$0.015$	$2.238\times {10}^{-6}$	$2.288\times {10}^{-6}$	$2.235\times {10}^{-6}$	$6.373\times {10}^{-2}$
	$0.02$	$1.580\times {10}^{-6}$	$1.616\times {10}^{-6}$	$1.575\times {10}^{-6}$	$0.116$
	$0.035$	$2.521\times {10}^{-6}$	$2.563\times {10}^{-6}$	$2.507\times {10}^{-6}$	$0.254$
	$0.05$	$3.361\times {10}^{-6}$	$3.432\times {10}^{-6}$	$4.037\times {10}^{-4}$	$0.341$
	$0.065$	$3.612\times {10}^{-6}$	$3.697\times {10}^{-6}$	$7.598\times {10}^5$	$0.401$
${\Delta }_{l^2}q\left(t_i\right)$	$0.01$	$1.458\times {10}^{-8}$	$1.478\times {10}^{-8}$	$1.454\times {10}^{-8}$	$8.528\times {10}^{-5}$
	$0.015$	$1.152\times {10}^{-8}$	$1.178\times {10}^{-8}$	$1.151\times {10}^{-8}$	$3.282\times {10}^{-4}$
	$0.02$	$9.128\times {10}^{-9}$	$9.334\times {10}^{-9}$	$9.101\times {10}^{-9}$	$6.708\times {10}^{-4}$
	$0.035$	$1.370\times {10}^{-8}$	$1.393\times {10}^{-8}$	$1.363\times {10}^{-8}$	$1.378\times {10}^{-3}$
	$0.05$	$1.491\times {10}^{-8}$	$1.522\times {10}^{-8}$	$1.790\times {10}^{-6}$	$1.514\times {10}^{-3}$
	$0.065$	$1.306\times {10}^{-8}$	$1.336\times {10}^{-8}$	$2.746\times {10}^3$	$1.448\times {10}^{-3}$
${\Delta }_{l^{\infty }}q$		$1.426\times {10}^{-5}$	$1.437\times {10}^{-5}$	$11.19$	$3.614\times {10}^{-3}$

Table 3. Table 3: Errors of numerically obtained solutions for temperature using Adams-Bashforth, centered with RAW filter, centered and Euler approximation schemes, with respect to analytically obtained response to Gaussian function in case of power-type distributed-order heat conduction law.

	$t_i$	Adams- Bashworth	Centered with RAW filter	Centered	Euler
${\delta }_{l^2}T\left(t_i\right)$	$0.01$	$6.577\times {10}^{-6}$	$6.613\times {10}^{-6}$	$6.603\times {10}^{-6}$	$3.284\times {10}^{-2}$
	$0.015$	$4.562\times {10}^{-6}$	$4.582\times {10}^{-6}$	$4.532\times {10}^{-6}$	$8.450\times {10}^{-2}$
	$0.02$	$8.527\times {10}^{-6}$	$8.558\times {10}^{-6}$	$8.427\times {10}^{-6}$	$0.136$
	$0.035$	$5.579\times {10}^{-6}$	$5.643\times {10}^{-6}$	$5.542\times {10}^{-6}$	$0.273$
	$0.05$	$9.058\times {10}^{-6}$	$9.173\times {10}^{-6}$	$9.112\times {10}^{-6}$	$0.394$
	$0.065$	$1.139\times {10}^{-5}$	$1.156\times {10}^{-5}$	$7.032\times {10}^{-5}$	$0.481$
${\Delta }_{l^2}T\left(t_i\right)$	$0.01$	$4.654\times {10}^{-8}$	$4.679\times {10}^{-8}$	$4.672\times {10}^{-8}$	$2.324\times {10}^{-4}$
	$0.015$	$2.690\times {10}^{-8}$	$2.701\times {10}^{-8}$	$2.671\times {10}^{-8}$	$4.982\times {10}^{-4}$
	$0.02$	$4.309\times {10}^{-8}$	$4.325\times {10}^{-8}$	$4.259\times {10}^{-8}$	$6.849\times {10}^{-4}$
	$0.035$	$2.129\times {10}^{-8}$	$2.154\times {10}^{-8}$	$2.115\times {10}^{-8}$	$1.040\times {10}^{-3}$
	$0.05$	$2.906\times {10}^{-8}$	$2.943\times {10}^{-8}$	$2.923\times {10}^{-8}$	$1.264\times {10}^{-3}$
	$0.065$	$3.155\times {10}^{-8}$	$3.203\times {10}^{-8}$	$1.948\times {10}^{-7}$	$1.334\times {10}^{-3}$
${\Delta }_{l^{\infty }}T$		$2.267\times {10}^{-5}$	$2.272\times {10}^{-5}$	$8.842\times {10}^{-5}$	$2.699\times {10}^{-3}$

Table 4. Table 4: Errors of numerically obtained solutions for heat flux using Adams-Bashforth, centered with RAW filter, centered and Euler approximation schemes, with respect to analytically obtained response to Gaussian function in case of power-type distributed-order heat conduction law.

	$t_i$	Adams- Bashworth	Centered with RAW filter	Centered	Euler
${\delta }_{l^2}q\left(t_i\right)$	$0.01$	$5.647\times {10}^{-6}$	$5.731\times {10}^{-6}$	$5.615\times {10}^{-6}$	$3.167\times {10}^{-2}$
	$0.015$	$3.488\times {10}^{-6}$	$3.579\times {10}^{-6}$	$3.478\times {10}^{-6}$	$9.435\times {10}^{-2}$
	$0.02$	$1.126\times {10}^{-5}$	$1.139\times {10}^{-5}$	$1.124\times {10}^{-5}$	$0.171$
	$0.035$	$8.189\times {10}^{-6}$	$8.355\times {10}^{-6}$	$8.157\times {10}^{-6}$	$0.360$
	$0.05$	$1.218\times {10}^{-5}$	$1.246\times {10}^{-5}$	$1.215\times {10}^{-5}$	$0.496$
	$0.065$	$1.481\times {10}^{-5}$	$1.517\times {10}^{-5}$	$1.989\times {10}^{-5}$	$0.595$
${\Delta }_{l^2}q\left(t_i\right)$	$0.01$	$4.683\times {10}^{-8}$	$4.753\times {10}^{-8}$	$4.657\times {10}^{-8}$	$2.627\times {10}^{-4}$
	$0.015$	$4.060\times {10}^{-8}$	$4.166\times {10}^{-8}$	$4.049\times {10}^{-8}$	$1.098\times {10}^{-3}$
	$0.02$	$1.466\times {10}^{-7}$	$1.483\times {10}^{-7}$	$1.463\times {10}^{-7}$	$2.228\times {10}^{-3}$
	$0.035$	$9.660\times {10}^{-8}$	$9.855\times {10}^{-8}$	$9.622\times {10}^{-8}$	$4.243\times {10}^{-3}$
	$0.05$	$1.189\times {10}^{-7}$	$1.216\times {10}^{-7}$	$1.186\times {10}^{-7}$	$4.838\times {10}^{-3}$
	$0.065$	$1.205\times {10}^{-7}$	$1.234\times {10}^{-7}$	$1.619\times {10}^{-7}$	$4.840\times {10}^{-3}$
${\Delta }_{l^{\infty }}q$		$4.375\times {10}^{-5}$	$4.399\times {10}^{-5}$	$4.368\times {10}^{-5}$	$6.074\times {10}^{-3}$

Equations

$$\rho c\frac{\partial }{\partial t}T\left(x,t\right)=-\frac{\partial }{\partial x}q\left(x,t\right),$$

$${\int }_0^1\varphi \left(\gamma \right){}_0^c{\mathrm{D}}_t^{\gamma }q\left(x,t\right)\mathrm{d}\gamma =-\lambda \frac{\partial }{\partial x}T\left(x,t\right),$$

$${}_0^c{\mathrm{D}}_t^{\gamma }y\left(t\right)=\frac{t^{-\gamma }}{\Gamma \left(1-\gamma \right)}{\ast }_t\frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t},$$

$$T\left(x,0\right)=T_0\left(x\right),q\left(x,0\right)=0,x\in \mathbb{R},$$

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

$$\phi\left(\gamma\right)=\tau_{0}\,\delta\left(\gamma-\alpha_{0}\right)+\sum_{\nu=1}^{N}\tau_{\nu}\,\delta\left(\gamma-\alpha_{\nu}\right),\;\;\begin{array}[]{c}0\leq\alpha_{0}<\ldots<\alpha_{N}<1,\\ \tau_{0},\tau_{1},\ldots,\tau_{N}>0,\end{array}$$

$$\varphi \left(\gamma \right)={\tau }^{\gamma },\tau >0.$$

$$\varphi \left(\gamma \right)=\delta \left(\gamma \right),\varphi \left(\gamma \right)=\tau \delta \left(\gamma -1\right)+\delta \left(\gamma \right),\varphi \left(\gamma \right)=\tau \delta \left(\gamma -\alpha \right)+\delta \left(\gamma \right),$$

$$\frac{\partial T}{\partial t}=\mathcal{D}\frac{{\partial }^{2}T}{\partial x^2},\tau \frac{{\partial }^{2}T}{\partial t^2}+\frac{\partial T}{\partial t}=\mathcal{D}\frac{{\partial }^{2}T}{\partial x^2}\text{and}\tau {}_0^c{\mathrm{D}}_t^{\alpha +1}T+\frac{\partial T}{\partial t}=\mathcal{D}\frac{{\partial }^{2}T}{\partial x^2},$$

$$\tau \frac{{\partial }^{2}T}{\partial t^2}+{}_0^c{\mathrm{D}}_t^{\alpha }T=\mathcal{D}\frac{{\partial }^{2}T}{\partial x^2},\text{with}\alpha \in \left(1,2\right),$$

$${}_0^c{\mathrm{D}}_t^{\alpha }T=\mathcal{D}\frac{{\partial }^{2}T}{\partial x^2},$$

$$\tau {}_0^c{\mathrm{D}}_t^{\alpha }T+{}_0^c{\mathrm{D}}_t^{\beta }T=\mathcal{D}\frac{{\partial }^{2}T}{\partial x^2},\text{or}\tau {}_0^c{\mathrm{D}}_t^{2\alpha }T+{}_0^c{\mathrm{D}}_t^{\alpha }T=\mathcal{D}\frac{{\partial }^{2}T}{\partial x^2},$$

$${\int }_0^2\varphi \left(\gamma \right){}_0^c{\mathrm{D}}_t^{\gamma }T\mathrm{d}\gamma =\mathcal{D}\frac{{\partial }^{2}T}{\partial x^2},$$

$$\bar{x}=\frac{x}{x^{\ast }},\bar{t}=\frac{t}{t^{\ast }},x^{\ast }=\sqrt{\frac{\lambda }{\rho c}{t}^{\ast }},\bar{T}=\frac{T}{{\Theta }_0}-1,\bar{q}=q\frac{1}{{\Theta }_0}\sqrt{\frac{t^{\ast }}{\lambda \rho c}},\bar{\varphi }=\frac{\varphi }{{\left(t^{\ast }\right)}^{\gamma }},$$

$$\frac{\partial }{\partial t}T\left(x,t\right)=-\frac{\partial }{\partial x}q\left(x,t\right),x\in \mathbb{R},t>0,$$

$${\int }_0^1\varphi \left(\gamma \right){}_0^c{\mathrm{D}}_t^{\gamma }q\left(x,t\right)\mathrm{d}\gamma =-\frac{\partial }{\partial x}T\left(x,t\right),x\in \mathbb{R},t>0,$$

$$\varphi \left(\gamma \right)=\delta \left(\gamma -{\alpha }_0\right)+\sum _{\nu =1}^N{\tau }_{\nu }\delta \left(\gamma -{\alpha }_{\nu }\right)\text{and}\varphi \left(\gamma \right)=1,$$

$$t^{\ast }={\tau }_0^{\frac{1}{{\alpha }_0}},\text{respectively}{t}^{\ast }=\tau ,$$

$$T\left(x,0\right)=T_0\left(x\right),q\left(x,0\right)=0,x\in \mathbb{R},$$

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

$$s\hat{\tilde{T}}\left(\xi ,s\right)-{\hat{T}}_0\left(\xi \right)=-\mathrm{i}\xi \hat{\tilde{q}}\left(\xi ,s\right),\xi \in \mathbb{R},\mathrm{Re}s>0,$$

$$\hat{\tilde{q}}\left(\xi ,s\right)\Phi \left(s\right)=-\mathrm{i}\xi \hat{\tilde{T}}\left(\xi ,s\right),\xi \in \mathbb{R},\mathrm{Re}s>0,$$

$$\Phi \left(s\right)={\int }_0^1\varphi \left(\gamma \right)s^{\gamma }\mathrm{d}\gamma ,\mathrm{Re}s>0,$$

$$\Phi \left(s\right)=s^{{\alpha }_0}+\sum _{\nu =1}^N{\tau }_{\nu }{s}^{{\alpha }_{\nu }}\text{and}\Phi \left(s\right)=\frac{s-1}{\ln s},\mathrm{Re}s>0.$$

$$\hat{\tilde{T}}\left(\xi ,s\right)={\hat{T}}_0\left(\xi \right)\frac{\Phi \left(s\right)}{{\xi }^2+s\Phi \left(s\right)}\text{and}\hat{\tilde{q}}\left(\xi ,s\right)=-{\hat{T}}_0\left(\xi \right)\frac{\mathrm{i}\xi }{{\xi }^2+s\Phi \left(s\right)}.$$

$${\mathcal{F}}^{-1}\left[\frac{1}{{\xi }^2+\lambda }\right]\left(x\right)=\frac{1}{2\sqrt{\lambda }}{\mathrm{e}}^{-|x|\sqrt{\lambda }},x\in \mathbb{R},\lambda \in \mathbb{C}\backslash \left(-\infty ,0\right],$$

$$T\left(x,t\right)=T_0\left(x\right){\ast }_{x}P\left(x,t\right)\text{and}q\left(x,t\right)=T_0\left(x\right){\ast }_{x}Q\left(x,t\right),$$

$$\tilde{P}\left(x,s\right)$$

$$\tilde{Q}\left(x,s\right)$$

$$f\left(t\right)=\frac{1}{2\pi \mathrm{i}}{\int }_{c-\mathrm{i}\infty }^{c+\mathrm{i}\infty }\tilde{f}\left(s\right){\mathrm{e}}^{st}\mathrm{d}s,t>0,$$

$$P\left(x,t\right)$$

$$Q\left(x,t\right)$$

$${\Phi }^{+}\left(p\right)=\Phi \left(p{\mathrm{e}}^{\mathrm{i}\pi }\right)\text{and}{\Phi }^{-}\left(p\right)=\Phi \left(p{\mathrm{e}}^{-\mathrm{i}\pi }\right),$$

Full text

Keywords: Cattaneo type heat conduction law, fractional distributed-order constitutive equation, integral transforms, finite differences

1 Introduction

Heat conduction in one-dimensional rigid material is considered on infinite spatial domain $x\in\mathbb{R}$ and for time $t>0.$ Generalization of the heat conduction equation is considered by treating two different processes in a material. The first one is heating, described by the energy balance equation

[TABLE]

where $\rho$ is used to denote the material density, $c$ is the specific heat capacity, while $T$ and $q$ denote temperature and heat flux respectively. The second one is heat conduction, described by the Cattaneo type time-fractional distributed-order heat conduction law

[TABLE]

where $\phi$ is the constitutive function or distribution, $\lambda$ is the thermal conductivity and ${}_{0}^{c}\mathrm{D}_{t}^{\gamma}$ is the operator of Caputo fractional differentiation of order $\gamma\in\left(0,1\right),$ defined by

[TABLE]

see [25], with $\ast_{t}$ denoting the convolution in time: $f\left(t\right)\ast_{t}g\left(t\right)=\int_{0}^{t}f\left(u\right)g\left(t-u\right)\mathrm{d}u.$

Rather than obtaining and solving a single heat conduction equation, the aim is to solve the system of equations consisting of energy balance equation (1) and constitutive equation (2), subject to initial

[TABLE]

and boundary conditions

[TABLE]

In particular, two special cases of the heat conduction law (2) will be examined: multi-term heat conduction law, obtained for the choice of constitutive distribution as

[TABLE]

consisting of at least two terms, where $\delta$ is used to denote the Dirac $\delta$-distribution, and power-type distributed-order heat conduction law, obtained for the choice of constitutive function as

[TABLE]

The constitutive equation (2) represents the generalization of known heat conduction laws such as Fourier, Cattaneo, fractional Cattaneo, which are obtained by choosing the constitutive distribution as

[TABLE]

respectively. The approach of considering generalized heat conduction equation through system of balance and constitutive equation is also adopted in [1] within the classical theory using the analogy with circuits and extending the results within the theory of fractional calculus in [15]. Anomalous transport processes through space and time fractional generalizations of the Cattaneo heat conduction law are studied in [13]. Time and space fractional heat conduction of Cattaneo type is studied, analytically on infinite domain in [2] and with physical justification for non-locality introduction in [8, 36, 55]. Heat conduction problem with the Riesz space fractional generalization of the Cattaneo-Christov heat conduction model is numerically treated in [26]. Heat conduction problems with different heat conduction laws in terms of the classical theory are reviewed in [24], while in [3] there is a collection of heat conduction problems within the theory of fractional calculus.

By combining the energy balance equation (1) with the constitutive equation (2), where constitutive distributions are given by (7), the classical heat conduction, telegraph and fractional telegraph equations are obtained as

[TABLE]

respectively, with the thermal diffusivity $\mathcal{D=}\frac{\lambda}{\rho c}.$ In [35], an unconditionally stable difference scheme is developed for (8)2 having the source term, while the equation of the form

[TABLE]

is solved for the Dirichlet boundary conditions in [39]. The equation of the form (8)3, with the additional (source) term corresponding to the impulse laser penetration into material causing its heating, is solved analytically in [41], while in [34] the similar problem with time-exponential decay of laser heating was treated numerically by developing unconditionally stable compact difference scheme.

Classical heat conduction equation (8)1 is often generalized within the theory of fractional calculus by replacing the first order partial time derivative with the fractional one, obtaining diffusion-wave equation

[TABLE]

describing subdiffusion if $\alpha\in\left(0,1\right)$ and superdiffusion if $\alpha\in\left(1,2\right).$ The pioneering work on diffusion-wave equation is [27], followed by further explorations in [16, 28]. Multi-dimensional variants of (9) are studied in [17, 18, 19] for space-fractional, space-time-fractional and time-fractional cases. Dirichlet problem for (9), using several finite difference schemes, such as explicit, implicit and Crank-Nicolson schemes is considered in [46]. Similar problem with the source term is numerically treated in [45], while an adaptive scheme was developed in [53]. Diffusion-wave equation (9) on bounded domain was numerically considered in [54] through implicit difference scheme, while the non-local variant of (9) is analyzed in [43] through explicit and implicit finite difference schemes. In order to numerically study anomalous infiltration in porous media the non-linear and variable-order version of (9) is used in [44].

The fractional generalization of telegraph equation (8)2 can be found in various forms different from (8)3 like

[TABLE]

with $0<\alpha<\beta<2,$ or $\alpha\in\left(0,1\right).$ Similarly as in [5], where (10)1 was treated analytically on infinite, semi-infinite and finite domains, in [48] the problem on the semi-infinite domain for (10)1 was treated for a special choice of the boundary conditions. By the use of analytical methods, different versions of telegraph equation are treated on unbounded and bounded domains in [11, 12, 22], including the non-locality as in [49]. The non-local version of (10)1 on unbounded domain, where the non-locality is expressed through the fractional Laplacian is analytically treated in [9, 10], while the non-local version of the telegraph equation (10)2 on infinite domain is treated in [40].

Generalizing the fractional telegraph equation by adding terms containing fractional derivatives of different orders lead to the distributed-order diffusion-wave equation

[TABLE]

analytically analyzed in [6, 7, 29, 30, 31]. A compact difference scheme for (11) on bounded domains, with source term included, is developed in [50], as well as in [37], where the similar equation is also numerically analyzed. In [32], equation (11) including the classical Laplacian is considered on a bounded multi-dimensional domain. Local, two-sided space-fractional, and Riesz space fractional variants of (11) are analyzed through the implicit finite difference schemes in [20, 21, 52]. Multi-term time-fractional diffusion type equation is considered in [42], while the maximum principle and numerical method for the multi-term time-space heat conduction equation of fractional order is considered in [51].

2 Solution to Cauchy problem

The Cauchy initial value problem on the real axis ($x\in\mathbb{R},$ $t>0$) for the time-fractional distributed-order Cattaneo type heat conduction, i.e., the system of energy balance equation (1) and constitutive Cattaneo type time-fractional distributed-order heat conduction law (2), subject to initial (3) and boundary conditions (4), will be analytically solved by the means of integral transform methods: Fourier transform with respect to spatial coordinate and Laplace transform with respect to time, as well as by the finite difference method: leap frog numerical scheme for spatial coordinate, along with Grünwald-Letnikov and third-order Adams-Bashforth temporal numerical schemes. The use of Adams-Bashforth scheme will prove to give more accurate and stable results when compared with centered, centered with RAW filter and Euler schemes. Two cases of the constitutive equation (2) will be examined: multi-term heat conduction law, with the constitutive distribution given by (5), and power-type distributed-order heat conduction law, with the constitutive function given by (6).

Dimensionless quantities

[TABLE]

where the time-scale $t^{\ast}$ will be determined according to the choice of constitutive distribution/function $\phi,$ see (15) below, and where the constant $\Theta_{0}$ represents the reference temperature, introduced into system of equations (1) and (2), with subsequent omittance of bars, yield the following form of governing equations

[TABLE]

while the constitutive distribution (5) and the constitutive function (6) become ($0\leq\alpha_{0}<\ldots<\alpha_{N}<1$)

[TABLE]

respectively. In the case of constitutive distribution (5), respectively constitutive function (6), the time-scales

[TABLE]

yield $\bar{\tau}_{\nu}=\tau_{\nu}\,\tau_{0}^{-\frac{\alpha_{\nu}}{\alpha_{0}}},$ $\nu=1,2,\ldots,N,$ where bar is omitted in (14), respectively $\bar{\tau}=1.$

Governing equations (12) and (13), with constitutive distribution/function (14), are subject to (dimensionless) initial and boundary conditions

[TABLE]

2.1 Analytical solution

Governing equations (12) and (13), with initial (16) and boundary conditions (17), will be analytically solved by the integral transform method. Application of the Fourier, $\hat{f}\left(\xi\right)=\mathcal{F}\left[f\left(x\right)\right]\left(\xi\right)=\int_{-\infty}^{\infty}f\left(x\right)\mathrm{e}^{-\mathrm{i}\xi x}\mathrm{d}x,$ $\xi\in\mathbb{R},$ and Laplace transform $\tilde{f}\left(s\right)=\mathcal{L}\left[f\left(t\right)\right]\left(s\right)=\int_{0}^{\infty}f\left(t\right)\mathrm{e}^{-st}\mathrm{d}t,$ $\mathop{\mathrm{R}e}s>0,$ to system of equations (12) and (13), with (16) and (17) taken into account, yields

[TABLE]

where

[TABLE]

in cases of constitutive distribution/function (14) takes the following forms

[TABLE]

Solution to system of equations (18) and (19) with respect to $\widehat{\tilde{T}}$ and $\widehat{\tilde{q}}$ reads ($\xi\in\mathbb{R},$ $\mathop{\mathrm{R}e}s>0$)

[TABLE]

Using the well-known Fourier inversion formula

[TABLE]

in (21), along with the Fourier transform of a derivative and convolution, one obtains

[TABLE]

after additional inversion of the Laplace transform, where ($x\in\mathbb{R},$ $\mathop{\mathrm{R}e}s>0$)

[TABLE]

with $\mathop{\mathrm{s}gn}x=2H\left(x\right)-1,$ $x\in\mathbb{R},$ and $H$ being the Heaviside function. The justification for using the Fourier inversion formula (22), as well as the argumentation that complex square root is well-defined, is given in Appendix A.

When the Laplace inversion formula

[TABLE]

is applied to $\tilde{P}$ and $\tilde{Q},$ given by (24) and (25), then solution kernels $P$ and $Q,$ appearing in (23), are obtained for $x\in\mathbb{R},$ $t>0,$ as

[TABLE]

where

[TABLE]

with $\Phi$ given by (20).

In order to obtain functions $P$ and $Q,$ (26) and (27), the Cauchy integral theorem

[TABLE]

is used, where $\Gamma=\Gamma_{1}\cup\Gamma_{2}\cup\Gamma_{3}\cup\Gamma_{4}\cup\Gamma_{5}\cup\Gamma_{6}\cup\Gamma_{7}\cup\Gamma_{8}$ is the closed contour shown in Figure 1, within $\tilde{P}\left(x,s\right)\mathrm{e}^{st},$ $x\in\mathbb{R},$ $t>0,$ is an analytic function, since apart from $s=0,$ there are no other branching points of $\tilde{P}$ and $\tilde{Q}.$

The contour integration will be performed for $\tilde{P}$ only, since the similar procedure and arguments are applicable for $\tilde{Q}$ as well. This will be shown in the sequel by proving that $\Phi,$ given by (20), has no zeros in the principal branch, i.e., for $\arg s\in\left(-\pi,\pi\right).$

The real and imaginary parts of $\Phi,$ given by (20)1, after substitution $s=\rho\,\mathrm{e}^{\mathrm{i}\varphi},$ $\rho>0,$ $\varphi\in\left(-\pi,\pi\right),$ read

[TABLE]

Since, by (30), it holds that $\mathop{\mathrm{I}m}\Phi\left(\bar{s}\right)=-\mathop{\mathrm{I}m}\Phi\left(s\right),$ where bar denotes the complex conjugation, it is sufficient to analyze function $\Phi$ for $\varphi\in\left[0,\pi\right)$ only. If $\varphi\in\left(0,\pi\right),$ then $\mathop{\mathrm{I}m}\Phi\left(\rho,\varphi\right)>0,$ since for $0\leq\alpha_{0}<\ldots<\alpha_{N}<1,$ it is valid that $\sin\left(\alpha_{\nu}\varphi\right)>0,$ $\nu=0,1,\ldots,N$ (except for possible $\alpha_{0}=0,$ which does affect that $\mathop{\mathrm{I}m}\Phi\left(\rho,\varphi\right)>0$). Note that $\mathop{\mathrm{I}m}\Phi\left(\rho,\varphi\right)=0,$ only if there is a single term with $\alpha_{0}=0$ in the constitutive distribution (14)1. This, however is not the case. If $\varphi=0,$ then, by (29), $\mathop{\mathrm{R}e}\Phi\left(\rho,\varphi\right)>0.$ Therefore, function $\Phi,$ (20)1, has no zeros in the principal branch.

In the case of $\Phi,$ given by (20)${}_{2},$ the corresponding real and imaginary parts are

[TABLE]

In order to determine zeros of $\Phi$, both real and imaginary part of $\Phi$ have to be zero, yielding the system of equations

[TABLE]

Combining the equations, one obtains

[TABLE]

If $\ln^{2}\rho+\varphi^{2}=0,$ then real and imaginary part of $\Phi$ are not well defined, see (31) and (32). Therefore, $\rho\cos\varphi-1=0,$ implying $\rho\varphi\sin\varphi=0,$ whose solutions are $\rho=0,$ $\varphi=0,$ $\varphi=\pi.$ The first and third solution are in contradiction with $\rho\cos\varphi-1=0$, while $\rho=1$ corresponds to the second one. By substituting $s=\left.\rho\,\mathrm{e}^{\mathrm{i}\varphi}\right|_{\varphi=0}$ in (20)2 one obtains $\frac{\rho-1}{\ln\rho}=0,$ while, on the other hand, $\lim_{\rho\rightarrow 1}\frac{\rho-1}{\ln\rho}=1,$ implying that $\rho=1$ is not a zero of $\Phi.$ Thus, there are no zeros in the principal branch of function $\Phi,$ (20)2, as well.

The integrals along contours $\Gamma_{1},$ $\Gamma_{4},$ parametrized by $s=p\,\mathrm{e}^{\mathrm{i}\pi},$ and $\Gamma_{6},$ parametrized by $s=p\,\mathrm{e}^{-\mathrm{i}\pi},$ in the limit when $R$ tends to infinity and $r$ tends to zero, yield ($x\in\mathbb{R},$ $t>0$)

[TABLE]

while the integrals along $\Gamma_{2},$ $\Gamma_{3},$ $\Gamma_{5},$ $\Gamma_{7}$ and $\Gamma_{8}$ are zero. Using aforementioned integrals in the Cauchy integral theorem (28) yields solution kernel $P$ in the form given by (26).

The absolute value of an integral along contour $\Gamma_{2},$ parametrized by $s=q+\mathrm{i}R,$ with $q\in\left(0,c\right)$ and $R$ tending to infinity, is estimated as

[TABLE]

Note that $\sqrt{q+\mathrm{i}R}\sim\sqrt{R}\mathrm{e}^{\mathrm{i}\frac{\pi}{4}},$ for $R\rightarrow\infty.$ For $\Phi$ given by (20)1, according to (29) and (30) and having in mind $0\leq\alpha_{0}<\ldots<\alpha_{N}<1,$ as $R\rightarrow\infty,$ one has

[TABLE]

while for $\Phi$ given by (20)2, according to (31) and (32), as $R\rightarrow\infty,$ one has

[TABLE]

so that, as $R\rightarrow\infty,$ in the first case of $\Phi,$ (33) becomes

[TABLE]

since $\cos\frac{\left(1+\alpha_{N}\right)\pi}{4}>0,$ while in the second case of $\Phi,$ (33) becomes

[TABLE]

implying in both cases $\lim_{R\rightarrow\infty}\int_{\Gamma_{2}}\tilde{P}\left(x,s\right)\mathrm{e}^{st}\mathrm{d}s=0.$ Similar arguments yield $\lim_{R\rightarrow\infty}\int_{\Gamma_{8}}\tilde{P}\left(x,s\right)\mathrm{e}^{st}\mathrm{d}s=0.$

On the contour $\Gamma_{3},$ parametrized by $s=R\,\mathrm{e}^{\mathrm{i}\varphi}$, $\varphi\in\left(\frac{\pi}{2},\pi\right),$ with $R\rightarrow\infty,$ the absolute value of the corresponding integral is estimated as

[TABLE]

For $\Phi$ given by (20)1, according to (29) and (30), as $R\rightarrow\infty,$ one has

[TABLE]

while for $\Phi$ given by (20)2, according to (31) and (32), as $R\rightarrow\infty,$ one has

[TABLE]

so that, as $R\rightarrow\infty,$ in the first case of $\Phi,$ (34) becomes

[TABLE]

while in the second case of $\Phi,$ (34) becomes

[TABLE]

In both cases $\cos\varphi<0$ and the first term in the exponential is of the highest order, implying $\lim_{R\rightarrow\infty}\int_{\Gamma_{3}}\tilde{P}\left(x,s\right)\mathrm{e}^{st}\mathrm{d}s=0.$ Similar arguments yield $\lim_{R\rightarrow\infty}\int_{\Gamma_{7}}\tilde{P}\left(x,s\right)\mathrm{e}^{st}\mathrm{d}s=0.$

The absolute value of the integral along the contour $\Gamma_{5},$ parametrized by $s=r\,\mathrm{e}^{\mathrm{i}\varphi}$, $\varphi\in\left(-\pi,\pi\right),$ with $r\rightarrow 0,$ is estimated as

[TABLE]

For $\Phi$ given by (20)1, according to (29) and (30), as $r\rightarrow 0,$ one has

[TABLE]

while for $\Phi$ given by (20)2, according to (31) and (32), as $r\rightarrow 0,$ one has

[TABLE]

so that, as $r\rightarrow 0,$ in the first case of $\Phi,$ (35) becomes

[TABLE]

while in the second case of $\Phi,$ (35) becomes

[TABLE]

implying in both cases $\lim_{r\rightarrow 0}\int_{\Gamma_{5}}\tilde{P}\left(x,s\right)\mathrm{e}^{st}\mathrm{d}s=0.$

2.2 Numerical solution through finite difference scheme

The Cauchy problem consisting of the (dimensionless) energy balance equation (12) and (dimensionless) constitutive Cattaneo type time-fractional distributed-order heat conduction law (13), corresponding to the time-fractional distributed-order Cattaneo type heat conduction, subject to initial (16) and boundary conditions (17), will be numerically solved, as the coupled system of equations, through the finite difference method, where discretization takes place on the spatial and temporal domain, with $\Delta x$ and $\Delta t$ being their steps, respectively.

The spatial derivatives will be approximated using the leap frog scheme

[TABLE]

which is standardly used second order accuracy scheme. The numerical calculations have shown that, unless using this scheme in discretization of all spatial derivatives appearing in the governing equations (12) and (13), the divergence occurs during the calculation of the solution. The possibility of using other schemes in discretization of all spatial derivatives in the governing equations was not considered.

The third order accuracy Adams-Bashforth scheme for the first order differential equation of the type

[TABLE]

reads

[TABLE]

and it will be used to approximate the energy balance equation (12), since it significantly dumps the computational mode for $\Delta t$ small enough, see [14]. The use of the third order Adams-Bashforth scheme will be justified in Section 3.3 by comparing its performance in stability and accuracy while calculating the solution to governing equations, with the first order accuracy Euler scheme

[TABLE]

the second order accuracy centered scheme

[TABLE]

and centered scheme with RAW filter, which is of the third order accuracy, thus improved when compared with the centered scheme. For the properties and implementation of the centered scheme with RAW filter see [47].

Although the Caputo fractional derivative appears in the constitutive equation (13), Grünwald-Letnikov approximation of the Riemann-Liouville derivative will be used, due to the zero initial condition (16)2, implying the equivalence of the Caputo and Riemann-Liouville fractional derivatives, so that, for $\alpha\in\left(0,1\right)$

[TABLE]

see [38]. When applying the Grünwald-Letnikov approximation, in order to avoid calculation of the binomial coefficients, i.e., gamma functions for large arguments, recurrence relation

[TABLE]

is adopted. Additionally, the integral appearing in the distributed-order constitutive law (13) will be approximated using the trapezoidal method.

Employing the leap frog scheme (36) in spatial domain and Adams-Bashforth scheme (37) in the energy balance equation (12), the approximation of governing equations (12), (13) reads

[TABLE]

with the initial conditions (16) yielding

[TABLE]

so that one obtains

[TABLE]

where

[TABLE]

Solution is found by solving the coupled system of equations (42), (43). In the first step ($n=0$), the initial conditions (40) used in (42) imply $T_{j}^{1}=T_{j}^{0},$ where it is assumed that $q_{j}^{-1}=q_{j}^{-2}=0,$ as well. The obtained $T_{j}^{1}$ is used in the second step ($n=1$), along with the initial condition (40)${}_{2},$ in (43), and $q_{j}^{1}$ is calculated. In this step, $T_{j}^{2}$ is calculated as well, according to (42). The algorithm for marching in time consists in repeating the calculations as in the second step. The system of equations (42), (43) is such that, while marching in time, the spatial domain shrinks due to the application of the leap frog scheme in the infinite spatial domain, since the scheme does not calculate values of solution in the outer points for each time step. The scheme follows the rule $j=2n-1,\ldots,J-\left(2n-1\right),$ $n=1,2,\ldots,$ for the heat flux and $j=2n,\ldots,J-2n,$ $n=1,2,\ldots,$ for the temperature, where $J$ is the last point in the domain in the first step.

If, instead of using the Adams-Bashforth scheme, one uses the Euler (38) and the centered scheme (39), the discretization of the energy balance equation (12), instead of (42), yields

[TABLE]

respectively. In the case of the centered scheme with RAW filter one should use (48) along with the procedure given in [47]. Solution is found by solving the coupled system of equations: either (46), (43) in the case of Euler discretization, or (48), (43) in the case of using the centered scheme. The coupled system of equations (46), (43) is solved by following the same procedure as described when the Adams-Bashforth scheme is used. When solving the system of equations (48), (43), in the first step ($n=0$), the Euler scheme is used in approximating the energy balance equation (12), i.e., (46) is assumed, implying $T_{j}^{1}=T_{j}^{0},$ according to the initial conditions (40). Following steps are same as already described for the case of using the Adams-Bashforth scheme.

The memory effects are clearly visible in (43), since the history of heat flux is taken into account through the weights $W_{k},$ $k=0,1,\ldots,n.$ The relation (44) shows that weights are dependent on the constitutive model. Namely, in the case of constitutive distribution (14)1, corresponding to the multi-term time-fractional heat conduction law, weights take the following form

[TABLE]

while in the case of constitutive distribution (14)2, corresponding to the power-type distributed-order heat conduction law, weights are

[TABLE]

It is clear that in the case of multi-term law no further approximation of weights (49) is needed, as opposed to the case of power-type distributed-order law, where the trapezoidal method for integral calculation, used for approximating weights (50), yields

[TABLE]

where $M=\frac{1}{\Delta\gamma}.$

3 Results

Temperature and heat flux as solutions to system of energy balance equation (12) and fractional distributed-order heat conduction law (13), with initial (16) and boundary conditions (17), analytically obtained by (23) and numerically by (42), (43), are plotted in cases of multi-term and power-type distributed-order laws, represented by constitutive distribution and function (14). Further, analytical response (23) to the initial temperature assumed as the Gaussian function and the solutions obtained using different numerical schemes for energy balance equation: (42), (46), and (48), along with the same approximation of the constitutive equation (43), are compared.

3.1 Analytically obtained temperature and heat flux

In order to show time evolution of temperature and heat flux spatial profiles, the initial temperature distribution is assumed as

[TABLE]

where $\delta$ is the Dirac delta distribution, so that temperature and heat flux, given by (23), reduce to the solution kernels $P$ and $Q,$ (26) and (27).

In the case of multi-term heat conduction law, the model parameters are taken as $\alpha_{0}=0,$ $\alpha_{1}=0.25,$ $\alpha_{2}=0.5,$ $\alpha_{3}=0.75,$ $\tau_{1}=0.4,$ $\tau_{2}=0.6,$ and $\tau_{3}=0.8,$ while the amplitude of the initial temperature distribution is $T_{0}=1.$ Figures 2

and 3

show the spatial profiles of temperature and heat flux at different time instances. The temperature distribution is symmetric with respect to the vertical axis, while the heat flows away from the origin (where the initial Dirac delta temperature distribution was introduced), therefore producing the antisymmetric character of the heat flux spatial profiles. One notices that spatial profiles of temperature have the similar behavior as in the case of the wave equation with energy dissipation effects included, see for example [4, Figures 3.2, 3.3, 3.6, and 3.7], unlike the heat conduction equations with fractional Cattaneo and Jeffreys heat conduction laws, see [3, Figures 7.3, 7.4, and 7.14]. Namely, as time passes, the peaks of both temperature and heat flux profiles propagate in space (which is characteristic for wave-like behavior) and decrease in height while increasing in width (which is characteristic for diffusion-like behavior in which case the peaks do not propagate). Therefore, heat conduction with multi-term heat conduction law, similarly as in the case of the classical Cattaneo constitutive law, might be considered as the propagation of heat waves. The diffusive characteristics of the process can be observed through the decrease of peaks’ height during time, as well as by the increase of their width.

Figures 4

and 5

show the spatial profiles of temperature and heat flux at different time instances in the case of power-type distributed-order heat conduction law, with $T_{0}=1.$ Likewise the case of heat conduction with multi-term law, one might consider this type of heat conduction as the propagation of heat waves. When compared to the case of heat conduction with multi-term law, the wave-like character of temperature and heat flux is more prominent, since the peaks are more localized (higher and narrower). However, the diffusion-like character is also noticeable, since the height of the peaks decreases while their width increases as time passes.

3.2 Numerically obtained temperature and heat flux

In order to be able to initialize the numerical scheme for simultaneous calculation of temperature and heat flux and to approximate the fundamental solution to heat conduction equation, whose spatial profiles are shown in Section 3.1, the initial temperature distribution is assumed as the Gaussian function

[TABLE]

as the smooth approximation of the Dirac delta distribution, which is obtained as $\varepsilon\rightarrow 0.$

The aim is to compare temperature and heat flux profiles, obtained analytically through (23) by convolving the solution kernels $P$ and $Q,$ (26) and (27), with the Gaussian function (52) as initial condition, with the ones obtained numerically through (42), (43), with weights (44) reducing to (49) in the case of multi-term heat conduction law and to (51) in the case of power-type distributed-order law.

Figures 6

and 7

present the comparison of temperature and heat flux spatial profiles obtained analytically and numerically in the case of multi-term heat conduction law, where $T_{0}=0.001$ and $\varepsilon=0.0005$ in the Gaussian function (52), the model parameters as taken as $\alpha_{0}=0,$ $\alpha_{1}=0.25,$ $\alpha_{2}=0.5,$ $\alpha_{3}=0.75,$ $\tau_{1}=0.4,$ $\tau_{2}=0.6,$ and $\tau_{3}=0.8,$ while the time and space steps in the numerical scheme are $\Delta t=0.0001$ and $\Delta x=0.001.$ From Figure 6 one sees that near the initial time instant, the temperature profiles resemble the Gaussian function and, as the process evolves, the profiles begin to resemble the profiles of the fundamental solution from Figure 2, thus concluding that the initial condition shapes the solution significantly in the beginning, while the characteristics of the process become prominent later on. This observation is supported by the heat flux profiles, whose peak height increases while the Gaussian function shapes the profiles, see Figure 7, and decreases, as expected from fundamental solution profiles on Figure 3, when the process becomes dominant.

In the case of power-type distributed-order law, the scheme requires an additional discretization of the integral in weights (50), reducing them to (51), where the integral discretization step is $\Delta\gamma=0.005.$ The other parameters are as in the case of multi-term law. Comparison of temperature and heat flux spatial profiles are shown in Figures 8

and 9.

Again, as in the case of multi-term law, near the initial time instant, the temperature profiles resemble the Gaussian function, see Figure 8, while later on, the profiles resemble the profiles of the fundamental solution from Figure 4. Also, the peaks’ height of heat flux profiles increase while the initial condition is dominant, see Figure 9, and decrease, as fundamental one does in Figure 5, when the process becomes dominant. As opposed to the case of multi-term law, in this case the process begins to shape the profiles for smaller times and the wave-like character is also more prominent.

Good agreement between analytical and numerical solution is evident from all Figures 6 - 9. The numerical scheme seems to be stable for selected time length, with the discretization steps and parameters taken.

3.3 Comparison of numerically obtained solutions

The aim is to test the accuracy and stability of numerical schemes using different approximations of the energy balance equation (12), while the constitutive heat conduction law (13) is approximated by (43) in all cases, by comparing the corresponding solutions among themselves and with the analytically obtained response (23) to the Gaussian function as initial temperature. Namely, the spatial derivatives are approximated by the leap frog scheme, the fractional derivative is used in the Grünwald-Letnikov form, while the approximations of the energy balance equation are obtained by using the following schemes: Adams-Bashforth (42), Euler (46), centered (48), and centered with RAW filter. The model parameters, as well as the parameters appearing in the Gaussian function, along with the discretization steps and time instances are throughout this section kept the same as in Section 3.2.

Figures 10

and 11

present comparison of temperature and heat flux spatial profiles obtained analytically according to (23) and numerically according to Euler (46), centered with RAW filter, and centered (48) scheme used in the energy balance equation, along with the multi-term heat conduction law approximated by (43), (49), as the responses to the Gaussian function. In the case of Euler approximation, the numerical scheme (46), (43) seems to be stable, however it yields inaccurate results for both temperature and heat flux when compared to the analytical solution, presumably due to the first order accuracy of the Euler scheme. Contrary to the previous case, when centered approximation is used, numerical scheme (48), (43) yields accurate results for both temperature and heat flux, however it becomes unstable by exhibiting high frequency oscillations having small amplitudes at $t=0.05$ and having large amplitudes at $t=0.065,$ thus not depicted on Figures 10 and 11. This is due to the existence of the computational mode that may not be damped in the three level schemes, see [33]. Having the accuracy also improved, centered scheme with RAW filter seems to yield stable results for the parameters and time interval considered, however it consumes more computational time when compared to numerical scheme (48), (43) that uses Adams-Bashforth approximation. Since the RAW filter averages values of the centered scheme, it damps the amplitudes of oscillations, thus yielding the stable solution.

Tables 1

and 2

contain absolute and relative errors of numerically obtained solutions for temperature $T_{ns}$ and heat flux $q_{ns}$ using Adams-Bashforth (42), centered with RAW filter, centered (48) and Euler (46) approximation scheme in the energy balance equation, along with the multi-term heat conduction law approximated by (43), (49), with respect to analytically obtained responses $T_{as}$ and $q_{as}$ to the initial Gaussian temperature distribution. The errors are calculated by using $l^{2}$ and $l^{\infty}$ norms according to

[TABLE]

where

[TABLE]

with $j_{\min}$ and $j_{\max}$ being dependant on the iteration step, as described in Section 2.2. For all time instances considered, the relative and absolute $l^{2}$ as well as $l^{\infty}$ errors for both temperature and heat flux, produced by the numerical scheme (42), (43), that uses the Adams-Bashforth approximation of the constitutive law, are smaller than the corresponding errors produced by the numerical scheme that uses the centered approximation with RAW filter. This, along with the reduced computational time in the case of using the Adams-Bashforth scheme implies its advantage. The scheme (48), (43), that uses the centered approximation of the constitutive equation, for time instances $t=0.01,\ldots,0.035,$ as opposed to the temperature, produces for the heat flux smaller values of the relative and absolute $l^{2}$ errors than the Adams-Bashforth scheme and centered scheme with RAW filter. The relative and absolute $l^{2}$ errors, for both temperature and heat flux, increase when there are high frequency oscillations with small amplitudes ($t=0.05$) and they become significantly large when the amplitudes increase ($t=0.065$), which is also evident from Figures 10 and 11. When using the centered scheme, the $l^{\infty}$ error is large due to the high frequency oscillations with large amplitudes. The scheme (46), (43), that uses the Euler approximation of the constitutive equation, produces the solution with the lowest accuracy, except for the case when (48), (43) exhibits high frequency oscillations, and the increasing error in each time instant.

Figures 12

and 13

present comparison of temperature and heat flux spatial profiles obtained analytically according to (23) and numerically according to Euler (46), centered with RAW filter, and centered (48) scheme used in the energy balance equation, along with the power-type distributed-order heat conduction law approximated by (43), (51), as the responses to the Gaussian function. Similarly as in the case of multi-term constitutive law, the use of Euler approximation gives solutions that seem to be stable, but inaccurate; the centered scheme also becomes unstable, but at the time instance later than the one in the case of multi-term law; the centered scheme with RAW filter seems to yield accurate and stable results for the parameters and time interval considered.

Tables 3

and 4

contain absolute and relative errors, calculated by (53), of numerically obtained solutions for temperature $T_{ns}$ and heat flux $q_{ns}$ using Adams-Bashforth (42), centered with RAW filter, centered (48) and Euler (46) approximation schemes in the energy balance equation, along with the power-type distributed-order heat conduction law approximated by (43), (51), with respect to analytically obtained response $T_{as}$ and $q_{as}$ to the initial Gaussian temperature distribution. Similarly as in the case of multi-term law, the use of Adams-Bashforth approximation produces smaller relative and absolute $l^{2}$ and $l^{\infty}$ errors when compared with the centered scheme with RAW filter; the Euler scheme also produces the solution with the lowest accuracy, having relative and absolute $l^{2}$ errors increasing with time. When compared with the Adams-Bashforth scheme and centered scheme with RAW filter, the use of centered scheme produces smaller values of the relative and absolute $l^{2}$ errors for both temperature and heat flux, except at time instances $t=0.01$ and $t=0.05$ for temperature and at $t=0.065$ for both temperature and heat flux. As one expects, the instability of scheme (48), (43) is not so clearly visible as in the case of multi-term law, since the solution has high frequency oscillations with small amplitudes, see Figures 12 and 13. However, the instability is implied by the difference of one order of magnitude in relative $l^{2}$ error of temperature between this scheme and the scheme that uses Adams-Bashforth approximation. This difference is not as prominent as in the case of heat flux. In the case of temperature, the scheme that uses Adams-Bashforth approximation has the smallest $l^{\infty}$ error, as opposed to the case of heat flux, where the smallest $l^{\infty}$ error is for the scheme that uses centered approximation.

4 Conclusion

The classical heat conduction equation is generalized by considering the system of equations consisting of the energy balance equation (1) and Cattaneo type time-fractional distributed-order constitutive heat conduction law (2). Two cases of the constitutive equation are examined: multi-term and power-type distributed-order heat conduction laws, with the constitutive distribution/function given by (5) and (6), respectively. The Cauchy initial value problem on the real axis is considered by subjecting governing equations (1) and (2) to initial and boundary conditions (3) and (4). Corresponding dimensionless system of equations (12) and (13), with (16) and (17), is solved analytically through integral transform methods: Fourier transform with respect to spatial coordinate and Laplace transform with respect to time, as well as by the finite difference method: leap frog numerical scheme for spatial coordinate, along with Grünwald-Letnikov and third-order Adams-Bashforth temporal numerical schemes. The analytical solution (23) is obtained as a convolution of initial temperature distribution with the solution kernels, given by (26) and (27), while the numerical solution is obtained through (42) and (43), with weights (44), reducing to (49) and (51) in the cases of multi-term and power-type distributed-order heat conduction laws, respectively. Note that solutions for temperature and heat flux naturally arise, since the scheme requires both temperature and heat flux for marching in time, due to the fact that the system of governing equations is coupled.

The response to the initial Dirac delta distribution yielded temperature and heat flux spatial profiles having the similar form as in the case of the telegraph equation, i.e., wave equation with energy dissipation effects included, see Figures 2 - 5, thus describing the propagation of heat waves, as opposed to the case of the heat conduction equations with fractional Cattaneo and Jeffreys heat conduction laws, having purely diffusive character.

Good agreement between analytical and numerical methods in cases of multi-term, see Figures 6 and 7, and power-type distributed-order heat conduction laws, see Figures 8 and 9, is found by comparing temperature and heat flux profiles, obtained analytically by convolving the solution kernels with the Gaussian function as initial condition and numerically through (42), (43), showing applicability of the joint use of Adams-Bashforth approximation of the energy balance equation, leap frog scheme for spatial derivatives, and Grünwald-Letnikov approximation of the fractional derivative.

Justification for the use of Adams-Bashforth scheme in approximating the energy balance equation is found by comparing the absolute and relative $l^{2}$ and $l^{\infty}$ errors (obtained with respect to the analytical solutions) of temperature and heat flux, produced by using the following schemes: Adams-Bashforth (42), Euler (46), centered (48), and centered with RAW filter, while the heat conduction law is in all cases approximated by (43). Namely, the Euler scheme proved to give stable, but inaccurate solutions, contrary to the centered scheme that yielded unstable, but the most accurate solutions for heat flux within the time domain of its stability, while the centered scheme with RAW filter gave stable solutions requiring longer computational time and having higher values of all errors. Therefore, the use of Adams-Bashforth scheme proved to be the best choice, both because of its accuracy and stability when compared with Euler, centered and centered with RAW filter scheme.

Appendix A Justification for applicability of the Fourier inversion formula (22)

In order to shown that the Fourier inversion formula (22) applies, one has to prove that $s\Phi\left(s\right)\in\mathbb{C}\backslash\left(-\infty,0\right]$ for $\mathop{\mathrm{R}e}s>0.$ More precisely, it will be shown that $\arg\left(s\Phi\left(s\right)\right)\in\left[0,\pi\right),$ for $\arg s\in\left[0,\frac{\pi}{2}\right)$ (and $\arg\left(s\Phi\left(s\right)\right)\in\left(-\pi,0\right),$ for $\arg s\in\left(\frac{\pi}{2},0\right)$) implying that the complex square root, $\sqrt{s\Phi\left(s\right)},$ for $\mathop{\mathrm{R}e}s>0,$ is well-defined.

In the case of constitutive distribution $\phi,$ given by (14)1, the substitution $s=\rho\,\mathrm{e}^{\mathrm{i}\varphi},$ $\rho>0,$ $\varphi\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right],$ implies that the real and imaginary parts of function $s\Phi\left(s\right),$ $\mathop{\mathrm{R}e}s>0,$ where $\Phi$ is given by (20)1, read

[TABLE]

Since, by (55), it holds that $\mathop{\mathrm{I}m}\left(\bar{s}\Phi\left(\bar{s}\right)\right)=-\mathop{\mathrm{I}m}\left(s\Phi\left(s\right)\right),$ where bar denotes the complex conjugation, it is sufficient to analyze function $s\Phi\left(s\right),$ $\mathop{\mathrm{R}e}s>0,$ for $\varphi\in\left[0,\frac{\pi}{2}\right]$ only. If $\varphi\in\left(0,\frac{\pi}{2}\right],$ then $\mathop{\mathrm{I}m}\left(s\Phi\left(s\right)\right)>0,$ since for $0\leq\alpha_{0}<\ldots<\alpha_{N}<1,$ it is valid that implying $\sin\left(\left(\alpha_{\nu}+1\right)\varphi\right)>0,$ $\nu=0,1,\ldots,N.$ If $\varphi=0,$ then, by (54), $\mathop{\mathrm{R}e}\left(s\Phi\left(s\right)\right)>0.$ Therefore, $s\Phi\left(s\right)\in\mathbb{C}\backslash\left(-\infty,0\right]$ and the Fourier inversion formula (22) applies, as well as that $\arg\left(s\Phi\left(s\right)\right)\in\left(0,\pi\right),$ for $\varphi\in\left(0,\frac{\pi}{2}\right).$

In the case of constitutive function $\phi,$ given by (14)${}_{2},$ it will be shown using the argument principle that function

[TABLE]

where $\Phi$ is given by (20)2, has no zeros in the right complex half-plane ($\mathop{\mathrm{R}e}s>0$) for any $\xi\in\mathbb{R},$ implying the applicability of the Fourier inversion formula (22). Moreover, it will also be shown that $\arg\left(\psi\left(s\right)\right)\in\left(0,\pi\right),$ for $\varphi\in\left(0,\frac{\pi}{2}\right).$ By substituting $s=\rho\,\mathrm{e}^{\mathrm{i}\varphi},$ $\rho>0,$ $\varphi\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right],$ in (56), the real and imaginary parts of function $\psi$ are obtained as

[TABLE]

Similarly as in the previous case, $\mathop{\mathrm{I}m}\left(\psi\left(\bar{s}\right)\right)=-\mathop{\mathrm{I}m}\left(\psi\left(s\right)\right),$ thus is sufficient to analyze function $\psi$ only in the upper right complex quarter-plane. In order to apply the argument principle, the contour $\Gamma=\gamma_{1}\cup\gamma_{2}\cup\gamma_{3}\cup\gamma_{4},$ shown in Figure 14, is used.

The contour $\gamma_{1}$ is parameterized by $s=x$, $x\in(r,R),$ with $r\rightarrow 0$ and $R\rightarrow\infty$, so that function $\psi,$ (56), reads

[TABLE]

since $x-1$ and $\ln x$ are of the same sign for $x\in\left(0,\infty\right).$ Moreover,

[TABLE]

The contour $\gamma_{2}$ is parametrized by $s=R\,\mathrm{e}^{\mathrm{i}\varphi}$, $\varphi\in\left(0,\frac{\pi}{2}\right),$ with $R\rightarrow\infty$. For $R$ sufficiently large, by (57) and (58), it is obtained

[TABLE]

In particular,

[TABLE]

Along $\gamma_{3}$, which is parametrized by $s=\rho\,\mathrm{e}^{\mathrm{i}\frac{\pi}{2}}$, $\rho\in\left(r,R\right)$, with $r\rightarrow 0$ and $R\rightarrow\infty$, by (57) and (58), the real and imaginary parts of $\psi$ read

[TABLE]

since $\rho>\ln\rho,$ for all $\rho\in\left(0,\infty\right).$ Also,

[TABLE]

The last part of the contour $\Gamma$ is the arc $\gamma_{4}$, parametrized by $s=r\,\mathrm{e}^{\mathrm{i}\varphi}$, $\varphi\in\left(0,\frac{\pi}{2}\right)$, with $r\rightarrow 0$. Again, (57) and (58), for $r$ sufficiently small, yield

[TABLE]

as well as

[TABLE]

Summing up, it is evident that as the complex variable $s$ changes along the contour $\Gamma,$ with $r$ tending to zero and $R$ tending to infinity, the imaginary part of function $\psi,$ (56), stays non-negative implying that $\arg\left(\psi\left(s\right)\right)\in\left(0,\pi\right),$ for $\varphi\in\left(0,\frac{\pi}{2}\right).$ This ensures the applicability of the Fourier inversion formula (22) and that the complex square root of $\psi$ is well-defined.

Acknowledgment

This work is supported by the Serbian Ministry of Science, Education and Technological Development under grant $174005$, as well as by the Provincial Government of Vojvodina under grant $114-451-2098$.

Figures presenting plots of solutions have been produced using Matplotlib, [23].