Analytical and numerical treatment of the heat conduction equation obtained via time-fractional distributed-order heat conduction law
aa r X i v : . [ m a t h . NA ] M a r Analytical and numerical treatment ofthe heat conduction equation obtained viatime-fractional distributed-order heat conduction law
Velibor ˇZeli ∗ , Duˇsan Zorica † September 4, 2018
Abstract
Generalization of the heat conduction equation is obtained by considering the system ofequations consisting of the energy balance equation and fractional-order constitutive heatconduction law, assumed in the form of the distributed-order Cattaneo type. The Cauchyproblem for system of energy balance equation and constitutive heat conduction law istreated analytically through Fourier and Laplace integral transform methods, as well asnumerically by the method of finite differences through Adams-Bashforth and Gr¨unwald-Letnikov schemes for approximation derivatives in temporal domain and leap frog schemefor spatial derivatives. Numerical examples, showing time evolution of temperature andheat flux spatial profiles, demonstrate applicability and good agreement of both methods incases of multi-term and power-type distributed-order heat conduction laws.
Keywords:
Cattaneo type heat conduction law, fractional distributed-order constitutiveequation, integral transforms, finite differences
Heat conduction in one-dimensional rigid material is considered on infinite spatial domain x ∈ R and for time t > . Generalization of the heat conduction equation is considered by treatingtwo different processes in a material. The first one is heating, described by the energy balanceequation ρc ∂∂t T ( x, t ) = − ∂∂x q ( x, t ) , (1)where ρ 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 bythe Cattaneo type time-fractional distributed-order heat conduction law Z φ ( γ ) c D γt q ( x, t ) d γ = − λ ∂∂x T ( x, t ) , (2) ∗ Linn´e FLOW Centre, KTH Mechanics, SE-100 44 Stockholm, Sweden, [email protected] † Mathematical Institute, Serbian Academy of Arts and Sciences, Kneza Mihaila 36, 11000 Beograd, Serbia,dusan [email protected] and Department of Physics, Faculty of Sciences, University of Novi Sad, Trg D.Obradovi´ca 3, 21000 Novi Sad, Serbia φ is the constitutive function or distribution, λ is the thermal conductivity and c D γt is theoperator of Caputo fractional differentiation of order γ ∈ (0 , , defined by c D γt y ( t ) = t − γ Γ (1 − γ ) ∗ t d y ( t )d t , see [25], with ∗ t denoting the convolution in time: f ( t ) ∗ t g ( t ) = R t f ( u ) g ( t − u ) d u. Rather than obtaining and solving a single heat conduction equation, the aim is to solvethe system of equations consisting of energy balance equation (1) and constitutive equation (2),subject to initial T ( x,
0) = T ( x ) , q ( x,
0) = 0 , x ∈ R , (3)and boundary conditions lim x →±∞ T ( x, t ) = 0 , lim x →±∞ q ( x, t ) = 0 , t > . (4)In particular, two special cases of the heat conduction law (2) will be examined: multi-term heatconduction law, obtained for the choice of constitutive distribution as φ ( γ ) = τ δ ( γ − α ) + N X ν =1 τ ν δ ( γ − α ν ) , ≤ α < . . . < α N < ,τ , τ , . . . , τ N > , (5)consisting of at least two terms, where δ is used to denote the Dirac δ -distribution, and power-type distributed-order heat conduction law, obtained for the choice of constitutive function as φ ( γ ) = τ γ , τ > . (6)The constitutive equation (2) represents the generalization of known heat conduction lawssuch as Fourier, Cattaneo, fractional Cattaneo, which are obtained by choosing the constitutivedistribution as φ ( γ ) = δ ( γ ) , φ ( γ ) = τ δ ( γ −
1) + δ ( γ ) , φ ( γ ) = τ δ ( γ − α ) + δ ( γ ) , (7)respectively. The approach of considering generalized heat conduction equation through systemof balance and constitutive equation is also adopted in [1] within the classical theory using theanalogy 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 Cattaneoheat conduction law are studied in [13]. Time and space fractional heat conduction of Cattaneotype 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 fractionalgeneralization 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 theoryare reviewed in [24], while in [3] there is a collection of heat conduction problems within thetheory of fractional calculus.By combining the energy balance equation (1) with the constitutive equation (2), whereconstitutive distributions are given by (7), the classical heat conduction, telegraph and fractionaltelegraph equations are obtained as ∂T∂t = D ∂ T∂x , τ ∂ T∂t + ∂T∂t = D ∂ T∂x and τ c D α +1 t T + ∂T∂t = D ∂ T∂x , (8)2espectively, with the thermal diffusivity D = λρc . In [35], an unconditionally stable differencescheme is developed for (8) having the source term, while the equation of the form τ ∂ T∂t + c D αt T = D ∂ T∂x , with α ∈ (1 , , is solved for the Dirichlet boundary conditions in [39]. The equation of the form (8) , with theadditional (source) term corresponding to the impulse laser penetration into material causingits heating, is solved analytically in [41], while in [34] the similar problem with time-exponentialdecay of laser heating was treated numerically by developing unconditionally stable compactdifference scheme.Classical heat conduction equation (8) is often generalized within the theory of fractionalcalculus by replacing the first order partial time derivative with the fractional one, obtainingdiffusion-wave equation c D αt T = D ∂ T∂x , (9)describing subdiffusion if α ∈ (0 ,
1) and superdiffusion if α ∈ (1 , . The pioneering work ondiffusion-wave equation is [27], followed by further explorations in [16, 28]. Multi-dimensionalvariants 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 termis numerically treated in [45], while an adaptive scheme was developed in [53]. Diffusion-waveequation (9) on bounded domain was numerically considered in [54] through implicit differencescheme, while the non-local variant of (9) is analyzed in [43] through explicit and implicit finitedifference schemes. In order to numerically study anomalous infiltration in porous media thenon-linear and variable-order version of (9) is used in [44].The fractional generalization of telegraph equation (8) can be found in various forms differentfrom (8) like τ c D αt T + c D βt T = D ∂ T∂x , or τ c D αt T + c D αt T = D ∂ T∂x , (10)with 0 < α < β < , or α ∈ (0 , . Similarly as in [5], where (10) was treated analyticallyon infinite, semi-infinite and finite domains, in [48] the problem on the semi-infinite domainfor (10) was treated for a special choice of the boundary conditions. By the use of analyticalmethods, different versions of telegraph equation are treated on unbounded and bounded domainsin [11, 12, 22], including the non-locality as in [49]. The non-local version of (10) on unboundeddomain, where the non-locality is expressed through the fractional Laplacian is analyticallytreated in [9, 10], while the non-local version of the telegraph equation (10) on infinite domainis treated in [40].Generalizing the fractional telegraph equation by adding terms containing fractional deriva-tives of different orders lead to the distributed-order diffusion-wave equation Z φ ( γ ) c D γt T d γ = D ∂ T∂x , (11)analytically analyzed in [6, 7, 29, 30, 31]. A compact difference scheme for (11) on boundeddomains, with source term included, is developed in [50], as well as in [37], where the similarequation is also numerically analyzed. In [32], equation (11) including the classical Laplacianis considered on a bounded multi-dimensional domain. Local, two-sided space-fractional, andRiesz space fractional variants of (11) are analyzed through the implicit finite difference schemesin [20, 21, 52]. Multi-term time-fractional diffusion type equation is considered in [42], while3he maximum principle and numerical method for the multi-term time-space heat conductionequation of fractional order is considered in [51]. The Cauchy initial value problem on the real axis ( x ∈ R , t >
0) for the time-fractionaldistributed-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 ofintegral transform methods: Fourier transform with respect to spatial coordinate and Laplacetransform with respect to time, as well as by the finite difference method: leap frog numericalscheme for spatial coordinate, along with Gr¨unwald-Letnikov and third-order Adams-Bashforthtemporal numerical schemes. The use of Adams-Bashforth scheme will prove to give more ac-curate and stable results when compared with centered, centered with RAW filter and Eulerschemes. Two cases of the constitutive equation (2) will be examined: multi-term heat conduc-tion law, with the constitutive distribution given by (5), and power-type distributed-order heatconduction law, with the constitutive function given by (6).Dimensionless quantities¯ x = xx ∗ , ¯ t = tt ∗ , x ∗ = s λρc t ∗ , ¯ T = T Θ − , ¯ q = q r t ∗ λρc , ¯ φ = φ ( t ∗ ) γ , where the time-scale t ∗ will be determined according to the choice of constitutive distribu-tion/function φ, see (15) below, and where the constant Θ represents the reference temperature,introduced into system of equations (1) and (2), with subsequent omittance of bars, yield thefollowing form of governing equations ∂∂t T ( x, t ) = − ∂∂x q ( x, t ) , x ∈ R , t > , (12) Z φ ( γ ) c D γt q ( x, t ) d γ = − ∂∂x T ( x, t ) , x ∈ R , t > , (13)while the constitutive distribution (5) and the constitutive function (6) become (0 ≤ α < . . . <α N < φ ( γ ) = δ ( γ − α ) + N X ν =1 τ ν δ ( γ − α ν ) and φ ( γ ) = 1 , (14)respectively. In the case of constitutive distribution (5), respectively constitutive function (6),the time-scales t ∗ = τ α , respectively t ∗ = τ , (15)yield ¯ τ ν = τ ν τ − ανα , ν = 1 , , . . . , N, where bar is omitted in (14), respectively ¯ τ = 1 . Governing equations (12) and (13), with constitutive distribution/function (14), are subjectto (dimensionless) initial and boundary conditions T ( x,
0) = T ( x ) , q ( x,
0) = 0 , x ∈ R , (16)lim x →±∞ T ( x, t ) = 0 , lim x →±∞ q ( x, t ) = 0 , t > . (17)4 .1 Analytical solution Governing equations (12) and (13), with initial (16) and boundary conditions (17), will be analyt-ically solved by the integral transform method. Application of the Fourier, ˆ f ( ξ ) = F [ f ( x )] ( ξ ) = R ∞−∞ f ( x ) e − i ξx d x, ξ ∈ R , and Laplace transform ˜ f ( s ) = L [ f ( t )] ( s ) = R ∞ f ( t ) e − st d t, Re s > , to system of equations (12) and (13), with (16) and (17) taken into account, yields s b ˜ T ( ξ, s ) − ˆ T ( ξ ) = − i ξ b ˜ q ( ξ, s ) , ξ ∈ R , Re s > , (18) b ˜ q ( ξ, s ) Φ ( s ) = − i ξ b ˜ T ( ξ, s ) , ξ ∈ R , Re s > , (19)where Φ ( s ) = Z φ ( γ ) s γ d γ, Re s > , in cases of constitutive distribution/function (14) takes the following formsΦ ( s ) = s α + N X ν =1 τ ν s α ν and Φ ( s ) = s − s , Re s > . (20)Solution to system of equations (18) and (19) with respect to b ˜ T and b ˜ q reads ( ξ ∈ R , Re s > b ˜ T ( ξ, s ) = ˆ T ( ξ ) Φ ( s ) ξ + s Φ ( s ) and b ˜ q ( ξ, s ) = − ˆ T ( ξ ) i ξξ + s Φ ( s ) . (21)Using the well-known Fourier inversion formula F − (cid:20) ξ + λ (cid:21) ( x ) = 12 √ λ e −| x |√ λ , x ∈ R , λ ∈ C \ ( −∞ , , (22)in (21), along with the Fourier transform of a derivative and convolution, one obtains T ( x, t ) = T ( x ) ∗ x P ( x, t ) and q ( x, t ) = T ( x ) ∗ x Q ( x, t ) , (23)after additional inversion of the Laplace transform, where ( x ∈ R , Re s > P ( x, s ) = 12 r Φ ( s ) s e −| x | √ s Φ( s ) , (24)˜ Q ( x, s ) = − p s Φ ( s ) dd x e −| x | √ s Φ( s ) = 12 e −| x | √ s Φ( s ) sgn x, (25)with sgn x = 2 H ( x ) − , x ∈ R , and H being the Heaviside function. The justification for usingthe Fourier inversion formula (22), as well as the argumentation that complex square root iswell-defined, is given in Appendix A.When the Laplace inversion formula f ( t ) = 12 π i Z c +i ∞ c − i ∞ ˜ f ( s ) e st d s, t > , is applied to ˜ P and ˜ Q, given by (24) and (25), then solution kernels P and Q, appearing in (23),are obtained for x ∈ R , t > , as P ( x, t ) = 14 π Z ∞ (cid:16)p Φ + ( p ) e − i | x | √ p Φ + ( p ) + p Φ − ( p ) e i | x | √ p Φ − ( p ) (cid:17) e − pt √ p d p, (26) Q ( x, t ) = sgn x π i Z ∞ (cid:16) e i | x | √ p Φ − ( p ) − e − i | x | √ p Φ + ( p ) (cid:17) e − pt d p, (27)5here Φ + ( p ) = Φ (cid:0) p e i π (cid:1) and Φ − ( p ) = Φ (cid:0) p e − i π (cid:1) , with Φ given by (20).In order to obtain functions P and Q, (26) and (27), the Cauchy integral theorem I Γ ˜ P ( x, s ) e st d s = 0 , (28)is used, where Γ = Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ is the closed contour shown in Figure 1,within ˜ P ( x, s ) e st , x ∈ R , t > , is an analytic function, since apart from s = 0 , there are no otherbranching points of ˜ P and ˜ Q. The contour integration will be performed for ˜ P only, since theFigure 1: Contour Γ.similar procedure and arguments are applicable for ˜ Q as well. This will be shown in the sequelby proving that Φ , given by (20), has no zeros in the principal branch, i.e., for arg s ∈ ( − π, π ) . The real and imaginary parts of Φ , given by (20) , after substitution s = ρ e i ϕ , ρ > ,ϕ ∈ ( − π, π ) , read Re Φ ( ρ, ϕ ) = ρ α cos ( α ϕ ) + N X ν =1 τ ν ρ α ν cos ( α ν ϕ ) , (29)Im Φ ( ρ, ϕ ) = ρ α sin ( α ϕ ) + N X ν =1 τ ν ρ α ν sin ( α ν ϕ ) . (30)Since, by (30), it holds that Im Φ (¯ s ) = − Im Φ ( s ) , where bar denotes the complex conjugation,it is sufficient to analyze function Φ for ϕ ∈ [0 , π ) only. If ϕ ∈ (0 , π ) , then Im Φ ( ρ, ϕ ) > , sincefor 0 ≤ α < . . . < α N < , it is valid that sin ( α ν ϕ ) > , ν = 0 , , . . . , N (except for possible α = 0 , which does affect that Im Φ ( ρ, ϕ ) > ρ, ϕ ) = 0 , only if there is asingle term with α = 0 in the constitutive distribution (14) . This, however is not the case. If ϕ = 0 , then, by (29), Re Φ ( ρ, ϕ ) > . Therefore, function Φ , (20) , has no zeros in the principalbranch.In the case of Φ , given by (20) , the corresponding real and imaginary parts areRe Φ ( ρ, ϕ ) = ln ρ ( ρ cos ϕ −
1) + ρϕ sin ϕ ln ρ + ϕ , (31)Im Φ ( ρ, ϕ ) = ρ ln ρ sin ϕ − ϕ ( ρ cos ϕ − ρ + ϕ . (32)6n order to determine zeros of Φ, both real and imaginary part of Φ have to be zero, yielding thesystem of equations ln ρ ( ρ cos ϕ −
1) + ρϕ sin ϕ = 0 ,ρ ln ρ sin ϕ − ϕ ( ρ cos ϕ −
1) = 0 . Combining the equations, one obtains (cid:0) ln ρ + ϕ (cid:1) ( ρ cos ϕ −
1) = 0 . If ln ρ + ϕ = 0 , then real and imaginary part of Φ are not well defined, see (31) and (32).Therefore, ρ cos ϕ − , implying ρϕ sin ϕ = 0 , whose solutions are ρ = 0 , ϕ = 0 , ϕ = π. Thefirst and third solution are in contradiction with ρ cos ϕ − ρ = 1 corresponds to thesecond one. By substituting s = ρ e i ϕ (cid:12)(cid:12) ϕ =0 in (20) one obtains ρ − ρ = 0 , while, on the otherhand, lim ρ → ρ − ρ = 1 , implying that ρ = 1 is not a zero of Φ . Thus, there are no zeros in theprincipal branch of function Φ , (20) , as well.The integrals along contours Γ , Γ , parametrized by s = p e i π , and Γ , parametrized by s = p e − i π , in the limit when R tends to infinity and r tends to zero, yield ( x ∈ R , t > R →∞ Z Γ ˜ P ( x, s ) e st d s = 2 π i P ( x, t ) , lim R →∞ ,r → Z Γ ˜ P ( x, s ) e st d s = −
12i lim R →∞ ,r → Z rR p Φ ( p e i π ) √ p e − i | x | √ p Φ( p e i π ) e − pt d p = 12i Z ∞ p Φ + ( p ) √ p e − i | x | √ p Φ + ( p ) e − pt d p, lim R →∞ ,r → Z Γ ˜ P ( x, s ) e st d s = 12i lim R →∞ ,r → Z Rr p Φ ( p e − i π ) √ p e i | x | √ p Φ( p e − i π ) e − pt d p = 12i Z ∞ p Φ − ( p ) √ p e i | x | √ p Φ − ( p ) e − pt d p, while the integrals along Γ , Γ , Γ , Γ and Γ are zero. Using aforementioned integrals in theCauchy integral theorem (28) yields solution kernel P in the form given by (26).The absolute value of an integral along contour Γ , parametrized by s = q +i R, with q ∈ (0 , c )and R tending to infinity, is estimated as (cid:12)(cid:12)(cid:12)(cid:12)Z Γ ˜ P ( x, s ) e st d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z c (cid:12)(cid:12)(cid:12)p Φ ( q + i R ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) √ q + i R (cid:12)(cid:12) (cid:12)(cid:12)(cid:12) e −| x |√ q +i R √ Φ( q +i R ) (cid:12)(cid:12)(cid:12) e qt d q. (33)Note that √ q + i R ∼ √ R e i π , for R → ∞ . For Φ given by (20) , according to (29) and (30) andhaving in mind 0 ≤ α < . . . < α N < , as R → ∞ , one hasRe Φ (cid:16) R, π (cid:17) ∼ R α N cos α N π (cid:16) R, π (cid:17) ∼ R α N sin α N π , i.e., p Φ ( q + i R ) ∼ R αN e i αN π , while for Φ given by (20) , according to (31) and (32), as R → ∞ , one hasRe Φ (cid:16) R, π (cid:17) ∼ π R ln R and Im Φ (cid:16) R, π (cid:17) ∼ R ln R , i.e., p Φ ( q + i R ) ∼ r R ln R e i π ,
7o that, as R → ∞ , in the first case of Φ , (33) becomes (cid:12)(cid:12)(cid:12)(cid:12)Z Γ ˜ P ( x, s ) e st d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z c R − αN e −| x | R αN cos ( αN ) π e qt d q ≤ , since cos (1+ α N ) π > , while in the second case of Φ , (33) becomes (cid:12)(cid:12)(cid:12)(cid:12)Z Γ ˜ P ( x, s ) e st d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z c √ ln R e qt d q ≤ , implying in both cases lim R →∞ R Γ ˜ P ( x, s ) e st d s = 0 . Similar arguments yield lim R →∞ R Γ ˜ P ( x, s ) e st d s =0 . On the contour Γ , parametrized by s = R e i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , with R → ∞ , the absolute valueof the corresponding integral is estimated as (cid:12)(cid:12)(cid:12)(cid:12)Z Γ ˜ P ( x, s ) e st d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z π π R (cid:12)(cid:12)(cid:12)p Φ (
R, ϕ ) (cid:12)(cid:12)(cid:12) √ R (cid:12)(cid:12)(cid:12) e −| x |√ R e i ϕ √ Φ( R,ϕ ) (cid:12)(cid:12)(cid:12) e Rt cos ϕ d ϕ. (34)For Φ given by (20) , according to (29) and (30), as R → ∞ , one hasRe Φ ( R, ϕ ) ∼ R α N cos ( α N ϕ ) and Im Φ ( R, ϕ ) ∼ R α N sin ( α N ϕ ) , i.e., p Φ (
R, ϕ ) ∼ R αN e i αN ϕ , while for Φ given by (20) , according to (31) and (32), as R → ∞ , one hasRe Φ ( R, ϕ ) ∼ R ln R cos ϕ and Im Φ ( R, ϕ ) ∼ R ln R sin ϕ, i.e., p Φ (
R, ϕ ) ∼ r R ln R e i ϕ , so that, as R → ∞ , in the first case of Φ , (34) becomes (cid:12)(cid:12)(cid:12)(cid:12)Z Γ ˜ P ( x, s ) e st d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z π π R αN e Rt cos ϕ −| x | R αN cos ( αN ) ϕ d ϕ ≤ , while in the second case of Φ , (34) becomes (cid:12)(cid:12)(cid:12)(cid:12)Z Γ ˜ P ( x, s ) e st d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z π π R √ ln R e Rt cos ϕ −| x | R √ ln R cos ϕ d ϕ ≤ . In both cases cos ϕ < R →∞ R Γ ˜ P ( x, s ) e st d s = 0 . Similar arguments yield lim R →∞ R Γ ˜ P ( x, s ) e st d s = 0 . The absolute value of the integral along the contour Γ , parametrized by s = r e i ϕ , ϕ ∈ ( − π, π ) , with r → , is estimated as (cid:12)(cid:12)(cid:12)(cid:12)Z Γ ˜ P ( x, s ) e st d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z π − π r (cid:12)(cid:12)(cid:12)p Φ ( r, ϕ ) (cid:12)(cid:12)(cid:12) √ r (cid:12)(cid:12)(cid:12) e −| x |√ r e i ϕ √ Φ( r,ϕ ) (cid:12)(cid:12)(cid:12) e rt cos ϕ d ϕ. (35)For Φ given by (20) , according to (29) and (30), as r → , one hasRe Φ ( r, ϕ ) ∼ r α cos ( α ϕ ) and Im Φ ( r, ϕ ) ∼ r α sin ( α ϕ ) , i.e., p Φ ( r, ϕ ) ∼ r α e i α ϕ , , according to (31) and (32), as r → , one hasRe Φ ( r, ϕ ) ∼ − r and Im Φ ( r, ϕ ) ∼ ϕ ln r , i.e., p Φ ( r, ϕ ) ∼ p | ln r | , so that, as r → , in the first case of Φ , (35) becomes (cid:12)(cid:12)(cid:12)(cid:12)Z Γ ˜ P ( x, s ) e st d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z π − π r α e rt cos ϕ −| x | r α cos (1+ α ϕ d ϕ ≤ , while in the second case of Φ , (35) becomes (cid:12)(cid:12)(cid:12)(cid:12)Z Γ ˜ P ( x, s ) e st d s (cid:12)(cid:12)(cid:12)(cid:12) ≤ Z π − π r r | ln r | e rt cos ϕ −| x | √ r | ln r | cos ϕ d ϕ ≤ , implying in both cases lim r → R Γ ˜ P ( x, s ) e st d s = 0 . The Cauchy problem consisting of the (dimensionless) energy balance equation (12) and (dimen-sionless) constitutive Cattaneo type time-fractional distributed-order heat conduction law (13),corresponding to the time-fractional distributed-order Cattaneo type heat conduction, subjectto initial (16) and boundary conditions (17), will be numerically solved, as the coupled system ofequations, through the finite difference method, where discretization takes place on the spatialand temporal domain, with ∆ x and ∆ t being their steps, respectively.The spatial derivatives will be approximated using the leap frog scheme ∂∂x y ( x, t ) (cid:12)(cid:12)(cid:12)(cid:12) x = j ∆ x, t = n ∆ t ≈ y nj +1 − y nj − x , (36)which is standardly used second order accuracy scheme. The numerical calculations have shownthat, unless using this scheme in discretization of all spatial derivatives appearing in the governingequations (12) and (13), the divergence occurs during the calculation of the solution. Thepossibility of using other schemes in discretization of all spatial derivatives in the governingequations was not considered.The third order accuracy Adams-Bashforth scheme for the first order differential equation ofthe type ∂∂t y ( x, t ) = f ( x, t ) , reads y n +1 j − y nj ∆ t = 112 (cid:0) f nj − f n − j + 5 f n − j (cid:1) , (37)and it will be used to approximate the energy balance equation (12), since it significantly dumpsthe computational mode for ∆ 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 andaccuracy while calculating the solution to governing equations, with the first order accuracyEuler scheme y n +1 j − y nj ∆ t = f nj , (38)9he second order accuracy centered scheme y n +1 j − y n − j t = f nj , (39)and centered scheme with RAW filter, which is of the third order accuracy, thus improved whencompared with the centered scheme. For the properties and implementation of the centeredscheme with RAW filter see [47].Although the Caputo fractional derivative appears in the constitutive equation (13), Gr¨unwald-Letnikov approximation of the Riemann-Liouville derivative will be used, due to the zero initialcondition (16) , implying the equivalence of the Caputo and Riemann-Liouville fractional deriva-tives, so that, for α ∈ (0 , c D γt y ( x, t )] x = j ∆ x, t = n ∆ t ≈ t ) γ y nj + 1(∆ t ) γ n X k =1 ω k ( γ ) y n − kj , ω k ( γ ) = ( − k (cid:18) γk (cid:19) , see [38]. When applying the Gr¨unwald-Letnikov approximation, in order to avoid calculation ofthe binomial coefficients, i.e., gamma functions for large arguments, recurrence relation ω k ( γ ) = (cid:18) − γk (cid:19) ω k − ( γ ) , ω ( γ ) = 1 , 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) inthe energy balance equation (12), the approximation of governing equations (12), (13) reads T n +1 j − T nj ∆ t = − q nj +1 − q nj − x − q n − j +1 − q n − j − x + 5 q n − j +1 − q n − j − x ! , Z φ ( γ ) t ) γ q nj + 1(∆ t ) γ n X k =1 ω k ( γ ) q n − kj ! d γ = − T nj +1 − T nj − x , with the initial conditions (16) yielding T j = ( T ) j , q j = 0 , (40)so that one obtains T n +1 j = T nj − ∆ t x × (cid:0) (cid:0) q nj +1 − q nj − (cid:1) − (cid:0) q n − j +1 − q n − j − (cid:1) + 5 (cid:0) q n − j +1 − q n − j − (cid:1)(cid:1) , n = 0 , , . . . , (41) q nj = − W T nj +1 − T nj − x + n X k =1 W k q n − kj ! , n = 1 , , . . . , (42)where W k = Z φ ( γ )(∆ t ) γ ω k ( γ ) d γ, k = 0 , , . . . , n. (43)Solution is found by solving the coupled system of equations (41), (42). In the first step( n = 0), the initial conditions (40) used in (41) imply T j = T j , where it is assumed that10 − j = q − j = 0 , as well. The obtained T j is used in the second step ( n = 1), along withthe initial condition (40) , in (42), and q j is calculated. In this step, T j is calculated as well,according to (41). The algorithm for marching in time consists in repeating the calculations asin the second step. The system of equations (41), (42) is such that, while marching in time,the spatial domain shrinks due to the application of the leap frog scheme in the infinite spatialdomain, since the scheme does not calculate values of solution in the outer points for each timestep. The scheme follows the rule j = 2 n − , . . . , J − (2 n − , n = 1 , , . . . , for the heat flux and j = 2 n, . . . , J − n, n = 1 , , . . . , for the temperature, where J is the last point in the domain inthe first step.If, instead of using the Adams-Bashforth scheme, one uses the Euler (38) and the centeredscheme (39), the discretization of the energy balance equation (12), instead of (41), yields T n +1 j = T nj − ∆ t x (cid:0) q nj +1 − q nj − (cid:1) , n = 0 , , . . . , (44) T n +1 j = T n − j − ∆ t ∆ x (cid:0) q nj +1 − q nj − (cid:1) , n = 1 , , . . . , (45)respectively. In the case of the centered scheme with RAW filter one should use (45) along withthe procedure given in [47]. Solution is found by solving the coupled system of equations: either(44), (42) in the case of Euler discretization, or (45), (42) in the case of using the centeredscheme. The coupled system of equations (44), (42) is solved by following the same procedureas described when the Adams-Bashforth scheme is used. When solving the system of equations(45), (42), in the first step ( n = 0), the Euler scheme is used in approximating the energy balanceequation (12), i.e., (44) is assumed, implying T j = T j , 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 (42), since the history of heat flux is taken intoaccount through the weights W k , k = 0 , , . . . , n. The relation (43) shows that weights aredependent on the constitutive model. Namely, in the case of constitutive distribution (14) ,corresponding to the multi-term time-fractional heat conduction law, weights take the followingform W k = 1(∆ t ) α ω k ( α ) + N X ν =1 τ ν (∆ t ) α ν ω k ( α ν ) , k = 0 , , . . . , n, (46)while in the case of constitutive distribution (14) , corresponding to the power-type distributed-order heat conduction law, weights are W k = Z t ) γ ω k ( γ ) d γ, k = 0 , , . . . , n. (47)It is clear that in the case of multi-term law no further approximation of weights (46) is needed,as opposed to the case of power-type distributed-order law, where the trapezoidal method forintegral calculation, used for approximating weights (47), yields W k = M − X m =0 t ) m +12 ∆ γ ω k (cid:18) m + 12 ∆ γ (cid:19) ∆ γ, k = 0 , , . . . , n, (48)where M = γ . Temperature and heat flux as solutions to system of energy balance equation (12) and fractionaldistributed-order heat conduction law (13), with initial (16) and boundary conditions (17), an-11lytically obtained by (23) and numerically by (41), (42), are plotted in cases of multi-term andpower-type distributed-order laws, represented by constitutive distribution and function (14).Further, analytical response (23) to the initial temperature assumed as the Gaussian functionand the solutions obtained using different numerical schemes for energy balance equation: (41),(44), and (45), along with the same approximation of the constitutive equation (42), are com-pared.
In order to show time evolution of temperature and heat flux spatial profiles, the initial temper-ature distribution is assumed as T ( x ) = T δ ( x ) , x ∈ R , where δ is the Dirac delta distribution, so that temperature and heat flux, given by (23), reduceto the solution kernels P and Q, (26) and (27).In the case of multi-term heat conduction law, the model parameters are taken as α = 0 ,α = 0 . , α = 0 . , α = 0 . , τ = 0 . , τ = 0 . , and τ = 0 . , while the amplitude of the initialtemperature distribution is T = 1 . Figures 2 and 3 show the spatial profiles of temperature and -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 x T ( x , t ) t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . Figure 2: Spatial profiles of temperature at different time instances, obtained analitically formulti-term heat conduction law, with Dirac distribution as initial condition.heat flux at different time instances. The temperature distribution is symmetric with respectto the vertical axis, while the heat flows away from the origin (where the initial Dirac deltatemperature distribution was introduced), therefore producing the antisymmetric character ofthe heat flux spatial profiles. One notices that spatial profiles of temperature have the similarbehavior as in the case of the wave equation with energy dissipation effects included, see forexample [4, Figures 3.2, 3.3, 3.6, and 3.7], unlike the heat conduction equations with fractionalCattaneo and Jeffreys heat conduction laws, see [3, Figures 7.3, 7.4, and 7.14]. Namely, astime passes, the peaks of both temperature and heat flux profiles propagate in space (which ischaracteristic for wave-like behavior) and decrease in height while increasing in width (which is12 x -60-40-200204060 q ( x , t ) t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . Figure 3: Spatial profiles of heat flux at different time instances, obtained analitically for multi-term heat conduction law, with Dirac distribution as initial condition.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 classicalCattaneo constitutive law, might be considered as the propagation of heat waves. The diffusivecharacteristics 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 timeinstances in the case of power-type distributed-order heat conduction law, with T = 1 . Likewisethe case of heat conduction with multi-term law, one might consider this type of heat conductionas the propagation of heat waves. When compared to the case of heat conduction with multi-termlaw, the wave-like character of temperature and heat flux is more prominent, since the peaks aremore 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.
In order to be able to initialize the numerical scheme for simultaneous calculation of temperatureand heat flux and to approximate the fundamental solution to heat conduction equation, whosespatial profiles are shown in Section 3.1, the initial temperature distribution is assumed as theGaussian function T ( x ) = T √ πε e − x ε , x ∈ R , (49)as the smooth approximation of the Dirac delta distribution, which is obtained as ε → . 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(49) as initial condition, with the ones obtained numerically through (41), (42), with weights(43) reducing to (46) in the case of multi-term heat conduction law and to (48) in the case ofpower-type distributed-order law. 13 x T ( x , t ) t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . Figure 4: Spatial profiles of temperature at different time instances, obtained analitically forpower-type distributed-order heat conduction law, with Dirac distribution as initial condition. -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 x -100-50050100 q ( x , t ) t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . Figure 5: Spatial profiles of heat flux at different time instances, obtained analitically for power-type distributed-order heat conduction law, with Dirac distribution as initial condition.14igures 6 and 7 present the comparison of temperature and heat flux spatial profiles obtained -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 x T ( x , t ) t = 0 t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . Figure 6: Comparison of spatial profiles of temperature at different time instances, obtained nu-merically (solid lines) and analytically (dots) for multi-term heat conduction law, with Gaussianfunction as initial condition.analytically and numerically in the case of multi-term heat conduction law, where T = 0 .
001 and ε = 0 . α = 0 , α = 0 . ,α = 0 . , α = 0 . , τ = 0 . , τ = 0 . , and τ = 0 . , while the time and space steps inthe numerical scheme are ∆ t = 0 . x = 0 . . From Figure 6 one sees that near theinitial time instant, the temperature profiles resemble the Gaussian function and, as the processevolves, the profiles begin to resemble the profiles of the fundamental solution from Figure 2, thusconcluding that the initial condition shapes the solution significantly in the beginning, while thecharacteristics of the process become prominent later on. This observation is supported by theheat 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, whenthe process becomes dominant.In the case of power-type distributed-order law, the scheme requires an additional discretiza-tion of the integral in weights (47), reducing them to (48), where the integral discretizationstep is ∆ γ = 0 . . The other parameters are as in the case of multi-term law. Comparison oftemperature and heat flux spatial profiles are shown in Figures 8 and 9. Again, as in the caseof multi-term law, near the initial time instant, the temperature profiles resemble the Gaussianfunction, see Figure 8, while later on, the profiles resemble the profiles of the fundamental solu-tion from Figure 4. Also, the peaks’ height of heat flux profiles increase while the initial conditionis dominant, see Figure 9, and decrease, as fundamental one does in Figure 5, when the processbecomes dominant. As opposed to the case of multi-term law, in this case the process begins toshape 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 stepsand parameters taken. 15 x -0.0050.0000.005 q ( x , t ) t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . Figure 7: Comparison of spatial profiles of heat flux at different time instances, obtained nu-merically (solid lines) and analytically (dots) for multi-term heat conduction law, with Gaussianfunction as initial condition. -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 x T ( x , t ) t = 0 t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . Figure 8: Comparison of spatial profiles of temperature at different time instances, obtainednumerically (solid lines) and analytically (dots) for power-type distributed-order heat conductionlaw, with Gaussian function as initial condition.16 x -0.010-0.0050.0000.0050.010 q ( x , t ) t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . Figure 9: Comparison of spatial profiles of heat flux at different time instances, obtained numer-ically (solid lines) and analytically (dots) for power-type distributed-order heat conduction law,with Gaussian function as initial condition.
The aim is to test the accuracy and stability of numerical schemes using different approximationsof the energy balance equation (12), while the constitutive heat conduction law (13) is approxi-mated by (42) in all cases, by comparing the corresponding solutions among themselves and withthe 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 isused in the Gr¨unwald-Letnikov form, while the approximations of the energy balance equationare obtained by using the following schemes: Adams-Bashforth (41), Euler (44), centered (45),and centered with RAW filter. The model parameters, as well as the parameters appearing inthe Gaussian function, along with the discretization steps and time instances are throughout thissection kept the same as in Section 3.2.Figures 10 and 11 present comparison of temperature and heat flux spatial profiles obtainedanalytically according to (23) and numerically according to Euler (44), centered with RAW filter,and centered (45) scheme used in the energy balance equation, along with the multi-term heatconduction law approximated by (42), (46), as the responses to the Gaussian function. In thecase of Euler approximation, the numerical scheme (44), (42) seems to be stable, however it yieldsinaccurate 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 (45), (42) yields accurate results for bothtemperature and heat flux, however it becomes unstable by exhibiting high frequency oscillationshaving small amplitudes at t = 0 .
05 and having large amplitudes at t = 0 . , thus not depictedon Figures 10 and 11. This is due to the existence of the computational mode that may not bedamped in the three level schemes, see [33]. Having the accuracy also improved, centered schemewith 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 (45), (42)that uses Adams-Bashforth approximation. Since the RAW filter averages values of the centered17 x T ( x , t ) −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 x T ( x , t ) −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 x T ( x , t ) Figure 10: Comparison of analytically obtained temperature profile (dots with dashed line) withnumerical one (solid line) for Euler (top left), centered with RAW filter (top right) and centered(bottom) scheme, as response to Gaussian function in case of multi-term heat conduction law. − . − . − . . . . . x − . . . q ( x , t ) − . − . − . . . . . x − . . . q ( x , t ) − . − . − . . . . . x − . . . q ( x , t ) Figure 11: Comparison of analytically obtained heat flux profile (dots with dashed line) withnumerical one (solid line) for Euler (top left), centered with RAW filter (top right) and centered(bottom) scheme, as response to Gaussian function in case of multi-term heat conduction law.18cheme, it damps the amplitudes of oscillations, thus yielding the stable solution.Tables 1 and 2 contain absolute and relative errors of numerically obtained solutions forTable 1: Errors of numerically obtained solutions for temperature using Adams-Bashforth, cen-tered with RAW filter, centered and Euler approximation schemes, with respect to analyticallyobtained response to Gaussian function in case of multi-term heat conduction law. t i Adams-Bashworth Centered withRAW filter Centered Euler δ l T ( t i ) 0 .
01 2 . × − . × − . × − . × − .
015 1 . × − . × − . × − . × − .
02 1 . × − . × − . × − . × − .
035 5 . × − . × − . × − . .
05 5 . × − . × − . × − . .
065 5 . × − . × − . × . l T ( t i ) 0 .
01 2 . × − . × − . × − . × − .
015 1 . × − . × − . × − . × − .
02 6 . × − . × − . × − . × − .
035 2 . × − . × − . × − . × − .
05 1 . × − . × − . × − . × − .
065 1 . × − . × − . × . × − ∆ l ∞ T . × − . × − .
69 2 . × − temperature T ns and heat flux q ns using Adams-Bashforth (41), centered with RAW filter, cen-tered (45) and Euler (44) approximation scheme in the energy balance equation, along with themulti-term heat conduction law approximated by (42), (46), with respect to analytically obtainedresponses T as and q as to the initial Gaussian temperature distribution. The errors are calculatedby using l and l ∞ norms according to∆ l u ( t i ) = k u as ( · , t i ) − u ns ( · , t i ) k l , δ l u ( t i ) = ∆ l u ( t i ) k u as ( · , t i ) k l , ∆ l ∞ u = k u as − u ns k l ∞ , (50)where k u ( · , t i ) k l = vuut j max − j min j max X j = j min ( u ( j ∆ x, t i )) , k u k l ∞ = max j min ≤ j ≤ j max , ≤ i ≤ u ( j ∆ x, t i ) , with j min and j max being dependant on the iteration step, as described in Section 2.2. For alltime instances considered, the relative and absolute l as well as l ∞ errors for both temperatureand heat flux, produced by the numerical scheme (41), (42), that uses the Adams-Bashforthapproximation of the constitutive law, are smaller than the corresponding errors produced bythe numerical scheme that uses the centered approximation with RAW filter. This, along withthe reduced computational time in the case of using the Adams-Bashforth scheme implies itsadvantage. The scheme (45), (42), that uses the centered approximation of the constitutiveequation, for time instances t = 0 . , . . . , . , as opposed to the temperature, produces for theheat flux smaller values of the relative and absolute l errors than the Adams-Bashforth schemeand centered scheme with RAW filter. The relative and absolute l errors, for both temperatureand heat flux, increase when there are high frequency oscillations with small amplitudes ( t = 0 . t i Adams-Bashworth Centered withRAW filter Centered Euler δ l q ( t i ) 0 .
01 3 . × − . × − . × − . × − .
015 2 . × − . × − . × − . × − .
02 1 . × − . × − . × − . .
035 2 . × − . × − . × − . .
05 3 . × − . × − . × − . .
065 3 . × − . × − . × . l q ( t i ) 0 .
01 1 . × − . × − . × − . × − .
015 1 . × − . × − . × − . × − .
02 9 . × − . × − . × − . × − .
035 1 . × − . × − . × − . × − .
05 1 . × − . × − . × − . × − .
065 1 . × − . × − . × . × − ∆ l ∞ q . × − . × − .
19 3 . × − and they become significantly large when the amplitudes increase ( t = 0 . l ∞ error is large dueto the high frequency oscillations with large amplitudes. The scheme (44), (42), that uses theEuler approximation of the constitutive equation, produces the solution with the lowest accuracy,except for the case when (45), (42) exhibits high frequency oscillations, and the increasing errorin each time instant.Figures 12 and 13 present comparison of temperature and heat flux spatial profiles obtainedanalytically according to (23) and numerically according to Euler (44), centered with RAWfilter, and centered (45) scheme used in the energy balance equation, along with the power-type distributed-order heat conduction law approximated by (42), (48), as the responses to theGaussian function. Similarly as in the case of multi-term constitutive law, the use of Eulerapproximation gives solutions that seem to be stable, but inaccurate; the centered scheme alsobecomes unstable, but at the time instance later than the one in the case of multi-term law; thecentered scheme with RAW filter seems to yield accurate and stable results for the parametersand time interval considered.Tables 3 and 4 contain absolute and relative errors, calculated by (50), of numerically ob-tained solutions for temperature T ns and heat flux q ns using Adams-Bashforth (41), centered withRAW filter, centered (45) and Euler (44) approximation schemes in the energy balance equation,along with the power-type distributed-order heat conduction law approximated by (42), (48),with respect to analytically obtained response T as and q as to the initial Gaussian temperaturedistribution. Similarly as in the case of multi-term law, the use of Adams-Bashforth approxima-tion produces smaller relative and absolute l and l ∞ errors when compared with the centeredscheme with RAW filter; the Euler scheme also produces the solution with the lowest accuracy,having relative and absolute l errors increasing with time. When compared with the Adams-Bashforth scheme and centered scheme with RAW filter, the use of centered scheme producessmaller values of the relative and absolute l errors for both temperature and heat flux, exceptat time instances t = 0 .
01 and t = 0 .
05 for temperature and at t = 0 .
065 for both tempera-20 x T ( x , t ) −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 x T ( x , t ) −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 x T ( x , t ) Figure 12: Comparison of analytically obtained temperature profile (dots with dashed line) withnumerical one (solid line) for Euler (top left), centered with RAW filter (top right) and centered(bottom) scheme, as response to Gaussian function in case of power-type distributed-order heatconduction law. − . − . − . . . . . x − . − . . . . q ( x , t ) − . − . − . . . . . x − . − . . . . q ( x , t ) − . − . − . . . . . x − . − . . . . q ( x , t ) Figure 13: Comparison of analytically obtained heat flux profile (dots with dashed line) withnumerical one (solid line) for Euler (top left), centered with RAW filter (top right) and centered(bottom) scheme, as response to Gaussian function in case of power-type distributed-order heatconduction law. 21able 3: Errors of numerically obtained solutions for temperature using Adams-Bashforth, cen-tered with RAW filter, centered and Euler approximation schemes, with respect to analyticallyobtained response to Gaussian function in case of power-type distributed-order heat conductionlaw. t i Adams-Bashworth Centered withRAW filter Centered Euler δ l T ( t i ) 0 .
01 6 . × − . × − . × − . × − .
015 4 . × − . × − . × − . × − .
02 8 . × − . × − . × − . .
035 5 . × − . × − . × − . .
05 9 . × − . × − . × − . .
065 1 . × − . × − . × − . l T ( t i ) 0 .
01 4 . × − . × − . × − . × − .
015 2 . × − . × − . × − . × − .
02 4 . × − . × − . × − . × − .
035 2 . × − . × − . × − . × − .
05 2 . × − . × − . × − . × − .
065 3 . × − . × − . × − . × − ∆ l ∞ T . × − . × − . × − . × − Table 4: Errors of numerically obtained solutions for heat flux using Adams-Bashforth, centeredwith RAW filter, centered and Euler approximation schemes, with respect to analytically ob-tained response to Gaussian function in case of power-type distributed-order heat conductionlaw. t i Adams-Bashworth Centered withRAW filter Centered Euler δ l q ( t i ) 0 .
01 5 . × − . × − . × − . × − .
015 3 . × − . × − . × − . × − .
02 1 . × − . × − . × − . .
035 8 . × − . × − . × − . .
05 1 . × − . × − . × − . .
065 1 . × − . × − . × − . l q ( t i ) 0 .
01 4 . × − . × − . × − . × − .
015 4 . × − . × − . × − . × − .
02 1 . × − . × − . × − . × − .
035 9 . × − . × − . × − . × − .
05 1 . × − . × − . × − . × − .
065 1 . × − . × − . × − . × − ∆ l ∞ q . × − . × − . × − . × − l error of temperature between this scheme and the scheme thatuses Adams-Bashforth approximation. This difference is not as prominent as in the case of heatflux. In the case of temperature, the scheme that uses Adams-Bashforth approximation has thesmallest l ∞ error, as opposed to the case of heat flux, where the smallest l ∞ error is for thescheme that uses centered approximation. The classical heat conduction equation is generalized by considering the system of equationsconsisting 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 distri-bution/function given by (5) and (6), respectively. The Cauchy initial value problem on the realaxis is considered by subjecting governing equations (1) and (2) to initial and boundary condi-tions (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 respectto spatial coordinate and Laplace transform with respect to time, as well as by the finite differ-ence method: leap frog numerical scheme for spatial coordinate, along with Gr¨unwald-Letnikovand third-order Adams-Bashforth temporal numerical schemes. The analytical solution (23) isobtained as a convolution of initial temperature distribution with the solution kernels, givenby (26) and (27), while the numerical solution is obtained through (41) and (42), with weights(43), reducing to (46) and (48) in the cases of multi-term and power-type distributed-order heatconduction 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 factthat the system of governing equations is coupled.The response to the initial Dirac delta distribution yielded temperature and heat flux spatialprofiles having the similar form as in the case of the telegraph equation, i.e., wave equation withenergy dissipation effects included, see Figures 2 - 5, thus describing the propagation of heatwaves, as opposed to the case of the heat conduction equations with fractional Cattaneo andJeffreys heat conduction laws, having purely diffusive character.Good agreement between analytical and numerical methods in cases of multi-term, see Fig-ures 6 and 7, and power-type distributed-order heat conduction laws, see Figures 8 and 9, isfound by comparing temperature and heat flux profiles, obtained analytically by convolving thesolution kernels with the Gaussian function as initial condition and numerically through (41),(42), showing applicability of the joint use of Adams-Bashforth approximation of the energy bal-ance equation, leap frog scheme for spatial derivatives, and Gr¨unwald-Letnikov approximationof the fractional derivative.Justification for the use of Adams-Bashforth scheme in approximating the energy balanceequation is found by comparing the absolute and relative l and l ∞ errors (obtained with re-spect to the analytical solutions) of temperature and heat flux, produced by using the followingschemes: Adams-Bashforth (41), Euler (44), centered (45), and centered with RAW filter, whilethe heat conduction law is in all cases approximated by (42). Namely, the Euler scheme provedto 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 thecentered scheme with RAW filter gave stable solutions requiring longer computational time and23aving higher values of all errors. Therefore, the use of Adams-Bashforth scheme proved to bethe best choice, both because of its accuracy and stability when compared with Euler, centeredand centered with RAW filter scheme. A Justification for applicability of the Fourier inversionformula (22)
In order to shown that the Fourier inversion formula (22) applies, one has to prove that s Φ ( s ) ∈ C \ ( −∞ ,
0] for Re s > . More precisely, it will be shown that arg ( s Φ ( s )) ∈ [0 , π ) , for arg s ∈ (cid:2) , π (cid:1) (and arg ( s Φ ( s )) ∈ ( − π, , for arg s ∈ (cid:0) π , (cid:1) ) implying that the complex square root, p s Φ ( s ) , for Re s > , is well-defined.In the case of constitutive distribution φ, given by (14) , the substitution s = ρ e i ϕ , ρ > ,ϕ ∈ (cid:2) − π , π (cid:3) , implies that the real and imaginary parts of function s Φ ( s ) , Re s > , where Φ isgiven by (20) , readRe ( s Φ ( s )) | s = ρ e i ϕ = ρ α +1 cos (( α + 1) ϕ ) + N X ν =1 τ ν ρ α ν +1 cos (( α ν + 1) ϕ ) , (51)Im ( s Φ ( s )) | s = ρ e i ϕ = ρ α +1 sin (( α + 1) ϕ ) + N X ν =1 τ ν ρ α ν +1 sin (( α ν + 1) ϕ ) . (52)Since, by (52), it holds that Im (¯ s Φ (¯ s )) = − Im ( s Φ ( s )) , where bar denotes the complex conjuga-tion, it is sufficient to analyze function s Φ ( s ) , Re s > , for ϕ ∈ (cid:2) , π (cid:3) only. If ϕ ∈ (cid:0) , π (cid:3) , thenIm ( s Φ ( s )) > , since for 0 ≤ α < . . . < α N < , it is valid that implying sin (( α ν + 1) ϕ ) > ,ν = 0 , , . . . , N. If ϕ = 0 , then, by (51), Re ( s Φ ( s )) > . Therefore, s Φ ( s ) ∈ C \ ( −∞ ,
0] and theFourier inversion formula (22) applies, as well as that arg ( s Φ ( s )) ∈ (0 , π ) , for ϕ ∈ (cid:0) , π (cid:1) . In the case of constitutive function φ, given by (14) , it will be shown using the argumentprinciple that function ψ ( s ) = s Φ ( s ) + ξ , s ∈ C , (53)where Φ is given by (20) , has no zeros in the right complex half-plane (Re s >
0) for any ξ ∈ R , implying the applicability of the Fourier inversion formula (22). Moreover, it will also be shownthat arg ( ψ ( s )) ∈ (0 , π ) , for ϕ ∈ (cid:0) , π (cid:1) . By substituting s = ρ e i ϕ , ρ > , ϕ ∈ (cid:2) − π , π (cid:3) , in (53),the real and imaginary parts of function ψ are obtained asRe ψ ( ρ, ϕ ) = ρ ln ρ ( ρ cos (2 ϕ ) − cos ϕ ) + ϕ ( ρ sin (2 ϕ ) − sin ϕ )ln ρ + ϕ + ξ , (54)Im ψ ( ρ, ϕ ) = ρ ln ρ ( ρ sin (2 ϕ ) − sin ϕ ) − ϕ ( ρ cos (2 ϕ ) − cos ϕ )ln ρ + ϕ . (55)Similarly as in the previous case, Im ( ψ (¯ s )) = − Im ( ψ ( s )) , thus is sufficient to analyze function ψ only in the upper right complex quarter-plane. In order to apply the argument principle, thecontour Γ = γ ∪ γ ∪ γ ∪ γ , shown in Figure 14, is used. The contour γ is parameterized by s = x , x ∈ ( r, R ) , with r → R → ∞ , so that function ψ, (53), reads ψ ( x ) = x ( x − x + ξ > , since x − x are of the same sign for x ∈ (0 , ∞ ) . Moreover, ψ ( x ) → ξ as x → ψ ( x ) → ∞ as x → ∞ . e s Im s g g g g Figure 14: Contour Γ.The contour γ is parametrized by s = R e i ϕ , ϕ ∈ (cid:0) , π (cid:1) , with R → ∞ . For R sufficiently large,by (54) and (55), it is obtainedRe ψ ( R, ϕ ) ∼ R ln R cos (2 ϕ ) + ξ and Im ψ ( R, ϕ ) ∼ R ln R sin (2 ϕ ) > . In particular, Re ψ ( R, → ∞ and Re ψ (cid:16) R, π (cid:17) → −∞ , as R → ∞ , Im ψ ( R,
0) = 0 and Im ψ (cid:16) R, π (cid:17) → ∞ , as R → ∞ . Along γ , which is parametrized by s = ρ e i π , ρ ∈ ( r, R ), with r → R → ∞ , by (54) and(55), the real and imaginary parts of ψ readRe ψ (cid:16) ρ, π (cid:17) = − ρ ρ ln ρ + π ln ρ + π + ξ and Im ψ (cid:16) ρ, π (cid:17) = ρ π ρ − ln ρ ln ρ + π > , since ρ > ln ρ, for all ρ ∈ (0 , ∞ ) . Also,Re ψ (cid:16) ρ, π (cid:17) → ξ , as ρ → ψ (cid:16) ρ, π (cid:17) → −∞ , as ρ → ∞ , Im ψ (cid:16) ρ, π (cid:17) → , as ρ → ψ (cid:16) ρ, π (cid:17) → ∞ , as ρ → ∞ . The last part of the contour Γ is the arc γ , parametrized by s = r e i ϕ , ϕ ∈ (cid:0) , π (cid:1) , with r → r sufficiently small, yieldRe ψ ( r, ϕ ) = − r ln r cos ϕ + ξ > ψ ( r, ϕ ) = − r ln r sin ϕ > , as well as Re ψ ( r, → ξ and Re ψ (cid:16) r, π (cid:17) → ξ , as r → , Im ψ ( r,
0) = 0 and Im ψ (cid:16) r, π (cid:17) → , as r → . Summing up, it is evident that as the complex variable s changes along the contour Γ , with r tending to zero and R tending to infinity, the imaginary part of function ψ, (53), stays non-negative implying that arg ( ψ ( s )) ∈ (0 , π ) , for ϕ ∈ (cid:0) , π (cid:1) . This ensures the applicability of theFourier inversion formula (22) and that the complex square root of ψ is well-defined.25 cknowledgment This work is supported by the Serbian Ministry of Science, Education and Technological Devel-opment under grant 174005, as well as by the Provincial Government of Vojvodina under grant114 − − References [1] J. Alvarez-Ramirez, G. Fernandez-Anaya, F. J. Valdes-Parada, and J. A. Ochoa-Tapia. Ahigh-order extension for the Cattaneo’s diffusion equation.
Physica A , 368:345–354, 2006.[2] T. M. Atanackovic, S. Konjik, Lj. Oparnica, and D. Zorica. The Cattaneo type space-timefractional heat conduction equation.
Continuum Mechanics and Thermodynamics , 24:293–311, 2012.[3] T. M. Atanackovic, S. Pilipovic, B. Stankovic, and D. Zorica.
Fractional Calculus withApplications in Mechanics: Vibrations and Diffusion Processes . Wiley-ISTE, London, 2014.[4] T. M. Atanackovic, S. Pilipovic, B. Stankovic, and D. Zorica.
Fractional Calculus withApplications in Mechanics: Wave Propagation, Impact and Variational Principles . Wiley-ISTE, London, 2014.[5] T. M. Atanackovic, S. Pilipovic, and D. Zorica. Diffusion wave equation with two fractionalderivatives of different order.
Journal of Physics A: Mathematical and Theoretical , 40:5319–5333, 2007.[6] T. M. Atanackovic, S. Pilipovic, and D. Zorica. Time distributed-order diffusion-wave equa-tion. I. Volterra type equation.
Proceedings of the Royal Society A: Mathematical, Physicaland Engineering Sciences , 465:1869–1891, 2009.[7] T. M. Atanackovic, S. Pilipovic, and D. Zorica. Time distributed-order diffusion-wave equa-tion. II. Applications of the Laplace and Fourier transformations.
Proceedings of the RoyalSociety A: Mathematical, Physical and Engineering Sciences , 465:1893–1917, 2009.[8] G. Borino, M. Di Paola, and M. Zingales. A non-local model of fractional heat conductionin rigid bodies.
European Physical Journal Special Topics , 193:173–184, 2011.[9] R. F. Camargo, E. Capelas de Oliveira, and J. Vaz Jr. On the generalized Mittag-Lefflerfunction and its application in a fractional telegraph equation.
Mathematical Physics Anal-ysis and Geometry , 15:1–16, 2012.[10] R. F. Camargo, A. O. Chiacchio, and E. Capelas de Oliveira. Differentiation to fractionalorders and the fractional telegraph equation.
Journal of Mathematical Physics , 49:033505–1–12, 2008.[11] R. C. Cascaval, E. C. Eckstein, C. L. Frota, and J. A. Goldstein. Fractional telegraphequations.
Journal of Mathematical Analysis and Applications , 276:145–159, 2002.[12] J. Chen, F. Liu, and V. Anh. Analytical solution for the time-fractional telegraph equationby the method of separating variables.
Journal of Mathematical Analysis and Applications ,338:1364–1377, 2008. 2613] A. Compte and R. Metzler. The generalized Cattaneo equation for the description of anoma-lous transport processes.
Journal of Physics A: Mathematical and General , 30:7277–7289,1997.[14] D. R. Durran. The third-order Adams-Bashforth method: An attractive alternative toleapfrog time differencing.
Monthly Weather Review , 119:702–720, 1991.[15] G. Fernandez-Anaya, F. J. Valdes-Parada, and J. Alvarez-Ramirez. On generalized fractionalCattaneo’s equations.
Physica A , 390:4198–4202, 2011.[16] R. Gorenflo, Y. Luchko, and F. Mainardi. Wright functions as scale-invariant solutions of thediffusion-wave equation.
Journal of Computational and Applied Mathematics , 118:175–191,2000.[17] A. Hanyga. Multidimensional solutions of space-fractional diffusion equations.
Proceedingsof the Royal Society A: Mathematical, Physical and Engineering Sciences , 457:2993–3005,2001.[18] A. Hanyga. Multi-dimensional solutions of space-time-fractional diffusion equations.
Pro-ceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , 458:429–450, 2002.[19] A. Hanyga. Multidimensional solutions of time-fractional diffusion-wave equations.
Proceed-ings of the Royal Society A: Mathematical, Physical and Engineering Sciences , 458:933–957,2002.[20] X. Hu, F. Liu, V. Anh, and I. Turner. A numerical investigation of the time distributed-orderdiffusion model.
ANZIAM , 55 (EMAC2013):C464–C478, 2014.[21] X. Hu, F. Liu, I. Turner, and V. Anh. An implicit numerical method of a new timedistributed-order and two-sided space-fractional advection-dispersion equation.
NumericalAlgorithms , 72:393–407, 2016.[22] F. Huang. Analytical solution for the time-fractional telegraph equation.
Journal of AppliedMathematics , 2009:890158–1–9, 2009.[23] J. D. Hunter. Matplotlib: A 2D graphics environment.
Computing in Science & Engineering ,9:90–95, 2007.[24] D. D. Joseph and L. Preziosi. Heat waves.
Reviews of Modern Physics , 61:41–73, 1989.[25] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo.
Theory and Applications of FractionalDifferential Equations . Elsevier B.V., Amsterdam, 2006.[26] L. Liu, L. Zheng, F. Liu, and X. Zhang. An improved heat conduction model with Riesz frac-tional Cattaneo-Christov flux.
International Journal of Heat and Mass Transfer , 103:1191–1197, 2016.[27] F. Mainardi. Fractional relaxation-oscillation and fractional diffusion-wave phenomena.
Chaos, Solitons & Fractals , 7:1461–1477, 1996.[28] F. Mainardi, Y. Luchko, and G. Pagnini. The fundamental solution of the space-timefractional diffusion equation.
Fractional Calculus and Applied Analysis , 4:153–192, 2001.2729] F. Mainardi, A. Mura, R. Gorenflo, and M. Stojanovic. The two forms of fractional relax-ation of distributed order.
Journal of Vibration and Control , 13:1249–1268, 2007.[30] F. Mainardi, G. Pagnini, and R. Gorenflo. Some aspects of fractional diffusion equations ofsingle and distributed order.
Applied Mathematics and Computation , 187:295–305, 2007.[31] F. Mainardi, G. Pagnini, A. Mura, and R. Gorenflo. Time-fractional diffusion of distributedorder.
Journal of Vibration and Control , 14:1267–1290, 2008.[32] M. M. Meerschaert, E. Nane, and P. Vellaisamy. Distributed-order fractional diffusions onbounded domains.
Journal of Mathematical Analysis and Applications , 379:216–228, 2011.[33] F. Mesinger and A. Arakawa.
Numerical methods used in atmospheric models , volume 1 of
GARP publications series No. 17 . World Meteorological Organization, International Councilof Scientific Unions, Geneve, 1976.[34] T. N. Mishra and K. N. Rai. Numerical solution of FSPL heat conduction equation foranalysis of thermal propagation.
Applied Mathematics and Computation , 273:1006–1017,2016.[35] R. K. Mohanty. An unconditionally stable difference scheme for the one-space-dimensionallinear hyperbolic equation.
Applied Mathematics Letters , 17:101–105, 2004.[36] M. S. Mongiov´ı and M. Zingales. A non-local model of thermal energy transport: Thefractional temperature equation.
International Journal of Heat and Mass Transfer , 67:593–601, 2013.[37] M. L. Morgado and M. Rebelo. Numerical approximation of distributed order reaction-diffusion equations.
Journal of Computational and Applied Mathematics , 275:216–227, 2015.[38] I. Podlubny.
Fractional Differential Equations . Academic Press, San Diego, 1999.[39] H. Qi and X. Guo. Transient fractional heat conduction with generalized Cattaneo model.
International Journal of Heat and Mass Transfer , 76:535–539, 2014.[40] H. Qi and X. Jiang. Solutions of the space-time fractional Cattaneo diffusion equation.
Physica A , 390:1876–1883, 2011.[41] HT. Qi, HY. Xu, and XW. Guo. The Cattaneo-type time fractional heat conduction equationfor laser heating.
Computers and Mathematics with Applications , 66:824–831, 2013.[42] M. R. Rapai´c and Z. D. Jeliˇci´c. Optimal control of a class of fractional heat diffusion systems.
Nonlinear Dynamics , 62:39–51, 2010.[43] S. Shen, F. Liu, and V. Anh. Numerical approximations and solution techniques for thespace-time Riesz-Caputo fractional advection-diffusion equation.
Numerical Algorithms ,56:383–403, 2011.[44] S. Shen, F. Liu, Q. Liu, and V. Anh. Numerical simulation of anomalous infiltration inporous media.
Numerical Algorithms , 68:443–454, 2015.[45] ZZ. Sun and X. Wu. A fully discrete difference scheme for a diffusion-wave system.
AppliedNumerical Mathematics , 56:193–209, 2006.[46] M. ˇZecov´a and J. Terp´ak. Heat conduction modeling by using fractional-order derivatives.
Applied Mathematics and Computation , 257:365–373, 2015.2847] P. D. Williams. A proposed modification to the Robert-Asselin time filter.
Monthly WeatherReview , 137:2538–2546, 2009.[48] HY. Xu, HT. Qi, and XY. Jiang. Fractional Cattaneo heat equation in a semi-infinitemedium.
Chinese Physics B , 22:014401–1–6, 2013.[49] S. Yakubovich and M. M. Rodrigues. Fundamental solutions of the fractional two-parametertelegraph equation.
Integral Transforms and Special Functions , 23:509–519, 2012.[50] H. Ye, F. Liu, and V. Anh. Compact difference scheme for distributed-order time-fractionaldiffusion-wave equation on bounded domains.
Journal of Computational Physics , 298:652–660, 2015.[51] H. Ye, F. Liu, V. Anh, and I. Turner. Maximum principle and numerical method for themulti-term time-space Riesz-Caputo fractional differential equations.
Applied Mathematicsand Computation , 227:531–540, 2014.[52] H. Ye, F. Liu, V. Anh, and I. Turner. Numerical analysis for the time distributed-order andRiesz space fractional diffusions on bounded domains.
IMA Journal of Applied Mathematics ,80:825–838, 2015.[53] S. B. Yuste and J. Quintana-Murillo. Fast, accurate and robust adaptive finite differencemethods for fractional diffusion equations.
Numerical Algorithms , 71:207–228, 2016.[54] P. Zhuang and F. Liu. Implicit difference approximation for the time fractional diffusionequation.
Journal of Applied Mathematics and Computing , 22:87–99, 2006.[55] M. Zingales. Fractional-order theory of heat transport in rigid bodies.