# Numerical Simulation of One-Dimensional Fractional Nonsteady Heat Transfer Model Based on the Second Kind Chebyshev Wavelet

Authors: Fuqiang Zhao; Jiaquan Xie; Qingxue Huang
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

In the current study, a numerical technique for solving one-dimensional fractional nonsteady heat transfer model is presented. We construct the second kind Chebyshev wavelet and then derive the operational matrix of fractional-order integration. The operational matrix of fractional-order integration is utilized to reduce the original problem to a system of linear algebraic equations, and then the numerical solutions obtained by our method are compared with those obtained by CAS wavelet method. Lastly, illustrated examples are included to demonstrate the validity and applicability of the technique.


---

## Body

## 1. Introduction

Fractional calculus is a branch of mathematics that deals with generalization of the well-known operations of differentiations to arbitrary orders. Many papers on fractional calculus have been published for the real-world applications in science and engineering such as viscoelasticity [1], bioengineering [2], biology [3], and more can be found in [4, 5]. Moreover fractional partial differential equations also are widely used in the areas of signal processing [6], mechanics [7], econometrics [8], fluid dynamics [9], and electromagnetics [10]. As the analytical solutions of fractional partial differential equations are not easy to derive, the scholars are committed to obtain their numerical solutions of these equations.In recent years, various numerical methods have been proposed for solving fractional diffusion equations, these methods include wavelets methods [11–17], Jacobi, Legendre, and Chebyshev polynomials methods [18–21], spectral methods [22, 23], finite element method [24], wavelet Galerkin method [25], and finite difference methods [26, 27]. In [28], a new matrix method is proposed to solve two-dimensional time-dependent diffusion equations with Dirichlet boundary conditions. In [29], the authors utilize the second kind Chebyshev wavelets to obtain the numerical solutions of the convection diffusion equations. Xie et al. use the Chebyshev operational matrix method to numerically solve one-dimensional fractional convection diffusion equations in [30]. In this paper, we apply the second kind Chebyshev wavelet method to obtain the numerical solutions of one-dimensional fractional nonsteady heat transfer model. The obtained numerical solutions by our method have been compared with those obtained by CAS wavelet method.The current paper is organized as follows: Section2 introduces the basic definitions of fractional calculus. In Section 3, the mathematical model of nonsteady heat transfer problem is proposed. Section 4 illustrates the second kind Chebyshev wavelets and their properties. In Section 5, we apply the second kind Chebyshev wavelet for solving fractional nonsteady heat transfer model. Numerical examples are presented to test the proposed method in Section 6. Finally, a conclusion is drawn in Section 7.


## 2. One-Dimensional Nonsteady Heat Transfer Model

For one infinite plate sample, as shown in Figure1, the height is δ, the upper surface and the edge are adiabatic, and the lower surface is contacted with the fluid, which its temperature is tw. The heat conductivity coefficient of the sample is λ, the density is ρ, and the specific heat capacity is cp. The initial temperature is t0, taking the origin of coordinates on the sample adiabatic surfaces, and the nonsteady heat transfer model with the initial-boundary condition can be defined as follows [31]: (1)∂t∂τ=λ∂2tρcp∂x2,τ=0,t=t0,x=0,∂t∂x=0,x=σ,t=tw.Obviously, when the sample density ρ, heat conductivity coefficient λ, specific heat capacity cp, and thickness δ are known, we can obtain the temperature distribution at any position x and any time τ, which is the nonsteady heat conduction model with constant temperature boundary condition. Based on the above-mentioned model, we give the fractional-order nonsteady heat transfer model of the following form:(2)∂T∂t=λ∂αTρcp∂xα+gx,t,0≤x≤1,t≥0,1<α≤2,with the initial condition:(3)Tx,0=fx,0≤x≤1,and the boundary conditions:(4)T0,t=g0t,T1,t=g1t,0≤t≤1,where gx,t denotes source term, fx is a given function, and g0t, g1t are continuous functions with first-order derivative.Figure 1
Nonsteady heat transfer model with constant temperature boundary condition.


## 3. Preliminaries of the Fractional Calculus

In this section, we give some necessary definitions and mathematical preliminaries on fractional calculus which will be used further in this paper.Definition 1.
The Riemann-Liouville fractional integral operatorIαα>0 of a function ft is defined as follows [4]:(5)Iαft=1Γα∫0tt-τα-1fτdτ,α>0,α∈R+.Some properties of the operator Iα are as follows:(6)IαIβft=Iα+βft,α>0,β>0,(7)Iαtγ=Γ1+γΓ1+γ+αtα+γ,γ>-1.Definition 2.
The Caputo fractional derivative0Dtα of a function ft is defined as follows [4]:(8)0Dtαft=1Γn-α∫0tfnτt-τn-α+1dτ,n-1<α≤n,n∈N.
Some properties of the Caputo fractional derivative are as follows:(9)0Dtαtβ=Γ1+βΓ1+β-αtβ-α,0<α<β+1,β>-1,IαDαft=ft-∑k=0n-1fk0+tkk!,n-1<α≤n,n∈N.


## 4. The Second Kind Chebyshev Wavelet and Its Operational Matrix of Fractional Integration

### 4.1. The Second Kind Chebyshev Wavelet and Its Properties

The second kind Chebyshev waveletψnmt=ψk,n,m,t has four arguments, n=1,2,…,2k-1,k∈N∗. They are defined on the interval 0,1 as follows [19]:(10)ψnmt=2k/2U~m2kt-2n+1,n-12k-1≤t<n2k-1,0,o.w.with(11)U~mt=2πUmt,m=0,1,2,…,M-1.Here Umt are the second kind Chebyshev polynomials which are orthogonal with respect to the weight function wt=1-t2 and satisfy the following recursive formula:(12)U0t=1,U1t=2t,Um+1t=2tUmt-Um-1t,m=1,2,….A functionft defined over 0,1 may be expanded in terms of the second kind Chebyshev wavelet as follows:(13)ft≃∑n=12k-1∑m=0M-1cnmψnmt=CTΨt,where(14)cnm=ft,ψnmtωn=∫01ωntψnmtdt,and the weight function wnt=w2kt-2n+1. Moreover, C and Ψt are m^=2k-1M column vectors given by(15)C=c10,c11,…,c1M-1,c20,c21,…,c2M-1,…,c2k-10,…,c2k-1M-1T,Ψt=ψ10,ψ11,…,ψ1M-1,ψ20,ψ21,…,ψ2M-1,…,ψ2k-10,…,ψ2k-1M-1T.Take the collocation points as follows:(16)ti=2i-12kM,i=1,2,…,2k-1M,m^=2k-1M.We define the second kind Chebyshev wavelet matrix Φm^×m^ as(17)Φm^×m^=Ψ12m^,Ψ32m^,…,Ψ2m^-12m^.An arbitrary function of two variablesTx,t defined over 0,1×0,1 may be expanded into Chebyshev wavelets basis as follows:(18)Tx,t≃∑i=1m^∑j=1m^dijψixψjt=ΨTxDΨt,where D=dijm^×m^ and dij=ψix,Tx,t,ψjt.The following theorem discusses the convergence and accuracy estimation of the proposed method.Theorem 3.
Letft be a second-order derivative square-integrable function defined over 0,1 with bounded second-order derivative, satisfying f′′t≤B for some constants B; then(1)
ft can be expanded as an infinite sum of the second kind Chebyshev wavelets and the series converge to ft uniformly, that is,(19)ft=∑n=0∞∑m∈Z∞cnmψnmt,wherecnm=ft,ψnmtLω20,1.(2)
(20)σf,k,M<πB23∑n=2k-1+1∞1n5∑m=M∞1m-141/2,whereσf,k,M=∫01ft-∑n=12k-1∑m=0M-1cnmψnmt2ωntdt1/2.


### 4.2. Operational Matrix of Fractional Integration

On the interval0,1, we defined a m^ – set of block-pulse functions (BPFs) as(21)bit=1,im^≤t<i+1m^,0,o.w.i=0,1,2,…,m^-1.The functions bit are disjoint and orthogonal:(22)bitbjt=0,i≠j,bit,i=j,∫01bisbjsds=0,i≠j,1m,i=j.Similarly, the second kind Chebyshev wavelet may be expanded into an m^-term block-pulse functions as(23)Ψt=Φm^×m^Bm^t.Kilicman has given the block-pulse functions operational matrix of fractional integrationFα of following form:(24)IαBm^t≈FαBm^t,where (25)Bm^t=b0t,b1t,…,bm^-1tT,Fα=1m^α1Γα+21ξ1ξ2ξ3⋯ξm^-101ξ1ξ2⋯ξm^-2001ξ1⋯ξm^-3⋮⋮⋱⋱⋮00⋯01ξ1000⋯01.Next, we derive the second kind Chebyshev wavelet operational matrix of fractional integration. Let(26)IαΨt=Pm^×m^αΨt,where Pm^×m^α is called the second kind Chebyshev wavelet operational matrix of fractional integration and it can be given by(27)Pm^×m^α=Φm^×m^FαΦm^×m^-1.For More details, see [29].


## 4.1. The Second Kind Chebyshev Wavelet and Its Properties

The second kind Chebyshev waveletψnmt=ψk,n,m,t has four arguments, n=1,2,…,2k-1,k∈N∗. They are defined on the interval 0,1 as follows [19]:(10)ψnmt=2k/2U~m2kt-2n+1,n-12k-1≤t<n2k-1,0,o.w.with(11)U~mt=2πUmt,m=0,1,2,…,M-1.Here Umt are the second kind Chebyshev polynomials which are orthogonal with respect to the weight function wt=1-t2 and satisfy the following recursive formula:(12)U0t=1,U1t=2t,Um+1t=2tUmt-Um-1t,m=1,2,….A functionft defined over 0,1 may be expanded in terms of the second kind Chebyshev wavelet as follows:(13)ft≃∑n=12k-1∑m=0M-1cnmψnmt=CTΨt,where(14)cnm=ft,ψnmtωn=∫01ωntψnmtdt,and the weight function wnt=w2kt-2n+1. Moreover, C and Ψt are m^=2k-1M column vectors given by(15)C=c10,c11,…,c1M-1,c20,c21,…,c2M-1,…,c2k-10,…,c2k-1M-1T,Ψt=ψ10,ψ11,…,ψ1M-1,ψ20,ψ21,…,ψ2M-1,…,ψ2k-10,…,ψ2k-1M-1T.Take the collocation points as follows:(16)ti=2i-12kM,i=1,2,…,2k-1M,m^=2k-1M.We define the second kind Chebyshev wavelet matrix Φm^×m^ as(17)Φm^×m^=Ψ12m^,Ψ32m^,…,Ψ2m^-12m^.An arbitrary function of two variablesTx,t defined over 0,1×0,1 may be expanded into Chebyshev wavelets basis as follows:(18)Tx,t≃∑i=1m^∑j=1m^dijψixψjt=ΨTxDΨt,where D=dijm^×m^ and dij=ψix,Tx,t,ψjt.The following theorem discusses the convergence and accuracy estimation of the proposed method.Theorem 3.
Letft be a second-order derivative square-integrable function defined over 0,1 with bounded second-order derivative, satisfying f′′t≤B for some constants B; then(1)
ft can be expanded as an infinite sum of the second kind Chebyshev wavelets and the series converge to ft uniformly, that is,(19)ft=∑n=0∞∑m∈Z∞cnmψnmt,wherecnm=ft,ψnmtLω20,1.(2)
(20)σf,k,M<πB23∑n=2k-1+1∞1n5∑m=M∞1m-141/2,whereσf,k,M=∫01ft-∑n=12k-1∑m=0M-1cnmψnmt2ωntdt1/2.


## 4.2. Operational Matrix of Fractional Integration

On the interval0,1, we defined a m^ – set of block-pulse functions (BPFs) as(21)bit=1,im^≤t<i+1m^,0,o.w.i=0,1,2,…,m^-1.The functions bit are disjoint and orthogonal:(22)bitbjt=0,i≠j,bit,i=j,∫01bisbjsds=0,i≠j,1m,i=j.Similarly, the second kind Chebyshev wavelet may be expanded into an m^-term block-pulse functions as(23)Ψt=Φm^×m^Bm^t.Kilicman has given the block-pulse functions operational matrix of fractional integrationFα of following form:(24)IαBm^t≈FαBm^t,where (25)Bm^t=b0t,b1t,…,bm^-1tT,Fα=1m^α1Γα+21ξ1ξ2ξ3⋯ξm^-101ξ1ξ2⋯ξm^-2001ξ1⋯ξm^-3⋮⋮⋱⋱⋮00⋯01ξ1000⋯01.Next, we derive the second kind Chebyshev wavelet operational matrix of fractional integration. Let(26)IαΨt=Pm^×m^αΨt,where Pm^×m^α is called the second kind Chebyshev wavelet operational matrix of fractional integration and it can be given by(27)Pm^×m^α=Φm^×m^FαΦm^×m^-1.For More details, see [29].


## 5. Numerical Implementation

In this section, we use the second kind Chebyshev wavelets method for numerically solving the nonsteady fractional-order heat transfer model with initial-boundary conditions. In order to solve this problem, we assume(28)∂3T∂t∂x2=ΨTxDΨt,where D=dijm^×m^ is an unknown matrix which should be determined, and Ψ· is the vector defined in (15). By integrating (28) from 0 to t, we obtain(29)∂2T∂x2=∂2T∂x2t=0+ΨTxDPm^×m^Ψt.Making use of the initial condition (3) enables one to put (29) in the following form:(30)∂2T∂x2=f′′x+ΨTxDPm^×m^Ψt. Then we have(31)∂αT∂xα=Ix2-α∂2T∂x2=Ix2-α∂2T∂x2t=0+ΨTxDPm^×m^Ψt=Ix2-αf′′x+ΨTxPm^×m^2-αTDPm^×m^Ψt.By integrating (30) two times from 0 to x, we obtain(32)Tx,t=T0,t+x∂T∂xx=0+fx-f0-xf′0+ΨTxPm^×m^2TDPm^×m^Ψt,and, by putting x=1 in (32), we get(33)Tx,t=T0,t+xHt+fx-f0-xf′0+ΨTxPm^×m^2TDPm^×m^Ψt,where(34)Ht=T1,t-T0,t+f0+f′0-f1-ΨT1Pm^×m^2TDPm^×m^Ψt.By one time differentiation of (33) with respect to t, we obtain(35)∂T∂t=T′0,t+xH′t+ΨTxPm^×m^2TDPm^×m^Ψt,where(36)H′t=T′1,t-T′0,t-ΨT1Pm^×m^2TDPm^×m^Ψt.Now by substituting (31) and (35) into (2) and combining (4) and taking collocation points xi=2i-1/m^, tj=2j-1/m^, i,j=1,2,3,…,m^, we obtain the following linear system of algebraic equations:(37)T′0,tj+xiT′1,tj-T′0,tj-ΨT1Pm^×m^2TDPm^×m^Ψtj+ΨTxiPm^×m^2TDPm^×m^Ψtj=aIx2-vf′′xi+aΨTxiPm^×m^2-vTDPm^×m^Ψtj+gxi,tj,i,j,=1,2,3,…,m^.By solving this system to determine D, we can get the numerical solution of this problem by substituting D into (33).


## 6. Numerical Simulations

In this section, we use the proposed method to solve the initial-boundary problem of nonsteady heat transfer equations. The following numerical examples are given to show the effectiveness and practicability of the proposed method and the results have been compared with the analytical solution.Example 4.
Consider the following fractional-order nonsteady heat transfer model:(38)∂T∂t=λ∂1.5Tρc-p∂x1.5+gx,t,0≤x≤1,t≥0,where the parameters ρ=7500, c-p=0.795, λ=800, and gx,t=xx-12t-1-0.302793571044498x0.5tt-1 with initial-boundary condition Tx,0=T0,t=T1,t=0. The analytical solution of this problem is Tx,t=xtx-1t-1. The graph of the analytical solution is shown in Figure 2. The graphs of the numerical solutions when k=M=3, k=M=4, k=M=5 are shown in Figures 3–5. From Examples 4, 6, and 7, it can be concluded that the numerical solutions approximate to the analytical solution for a given value k, as M increases, or, for a given value M, as k increases.Figure 2
Analytical solution.Figure 3
Numerical solution withk=M=3.Figure 4
Numerical solution withk=M=4.Figure 5
Numerical solution withk=M=5.Example 5.
Consider the following fractional-order nonsteady heat transfer equation:(39)∂T∂t=∂1.8T∂x1.8+2x2t-2x0.2t2Γ1.2,0≤x≤1,t≥0,with initial-boundary condition Tx,0=T0,t=0, T1,t=t2. The analytical solution of this problem is Tx,t=x2t2. When k=M=3, k=M=4, k=M=5, the numerical solutions obtained by our method and those obtained by CAS wavelet method at some values of x, t are listed in Table 1.Table 1
The numerical solutions obtained by our method and those obtained by CAS wavelet method whenk=M=3, k=M=4, k=M=5.
t
x
Anal. Sol.
k
=
M
=
3
k
=
M
=
4
k
=
M
=
5
Our method
CAS wavelet
Our method
CAS wavelet
Our method
CAS wavelet
0.2
0.3
0.0036000
0.00362673
0.01527126
0.00360257
0.00471281
0.00360019
0.00382719
0.6
0.0144000
0.01445390
0.03638127
0.01440370
0.01673180
0.01440048
0.01492319
0.9
0.0324000
0.03248217
0.07371928
0.03240631
0.03631963
0.03240060
0.03273187
0.5
0.3
0.0225000
0.02253176
0.04872121
0.02250487
0.02826189
0.02250046
0.02293819
0.6
0.0900000
0.09061074
0.12739812
0.09006721
0.09537428
0.09000059
0.09072347
0.9
0.2025000
0.20257431
0.25873179
0.20250850
0.20736183
0.20250074
0.20301829
0.8
0.3
0.0576000
0.05765362
0.09381981
0.05760489
0.60121872
0.05760062
0.05830218
0.6
0.2304000
0.23048904
0.28237189
0.23040790
0.23833829
0.23040074
0.23138192
0.9
0.5184000
0.51851904
0.60381038
0.51841027
0.52478172
0.51840112
0.51953785Example 6.
We consider the following second-order nonsteady heat transfer model:(40)∂T∂t=2∂2T∂x2+3sin⁡x-sin⁡t-2cos⁡t,0≤x≤1,t>0,in such a way that Tx,0=sin⁡x+1, T0,t=cos⁡t, T1,t=sin⁡1+cos⁡t. The analytical solution of the system is Tx,t=sin⁡x+cos⁡t. The absolute errors between the numerical and analytical solutions obtained by our method and CAS wavelet method at some values of x, t when k=3, M=3,M=4,M=5 are shown in Table 2. Table 2 shows that our method has a better approximation than CAS wavelet method.Table 2
The absolute errors obtained by our method and CAS wavelet method whenM=3, M=4, M=5.
x
,
t
Anal. Sol.
k
=
3, M=3
k
=
3, M=4
k
=
3, M=5
Our method
CAS wavelet
Our method
CAS wavelet
Our method
CAS wavelet
(
0,0
)
1.00000000
1.627162e-4
2.381923e-2
8.719295
e
-
6
8.737819
e
-
4
2.319280
e
-
6
2.648278
e
-
4
(
0.1,0.1
)
1.09483758
1.738173e-4
2.731899e-2
5.371912
e
-
6
6.271928
e
-
4
2.842802
e
-
6
3.748217
e
-
4
(
0.2,0.2
)
1.17873590
2.371827e-4
3.759289e-2
2.361827
e
-
5
3.271929
e
-
3
4.830209
e
-
6
4.684278
e
-
4
(
0.3,0.3
)
1.25085669
2.731872e-4
4.542767e-2
4.731872
e
-
5
4.281912
e
-
3
5.371982
e
-
6
6.472938
e
-
4
(
0.4,0.4
)
1.31047933
3.261772e-4
5.251757e-2
5.219289
e
-
5
5.381018
e
-
3
7.381928
e
-
7
7.863982
e
-
4
(
0.5,0.5
)
1.35700810
8.271985e-5
4.378391e-2
6.319288
e
-
5
6.379843
e
-
3
6.238299
e
-
6
7.635176
e
-
4
(
0.6,0.6
)
1.38997808
4.268278e-4
8.373456e-3
5.738273
e
-
5
5.792808
e
-
3
8.302930
e
-
6
8.368386
e
-
4
(
0.7,0.7
)
1.40905987
4.791982e-4
7.371928e-2
7.382093
e
-
5
7.728732
e
-
3
1.983100
e
-
5
9.673817
e
-
4
(
0.8,0.8
)
1.41406280
5.281928e-4
6.367643e-2
8.382938
e
-
5
8.732763
e
-
3
9.381098
e
-
6
2.371927
e
-
3
(
0.9,0.9
)
1.40493687
6.782916e-4
7.371892e-2
9.381982
e
-
5
9.738273
e
-
3
2.381983
e
-
5
4.281988
e
-
3
(
1.0,1.0
)
1.38177329
9.381928e-4
8.263828e-2
9.983787
e
-
5
9.425146
e
-
3
3.313910
e
-
5
8.871999
e
-
4Example 7.
Consider the following second-order nonsteady heat transfer model:(41)∂T∂t=λ∂2Tρcp∂x2+gx,t,0≤x≤1,t≥0,where the parameters ρ=7500, c-p=0.795, λ=1000, and gx,t=-πsin⁡πxsin⁡πt+0.167714884696017π2sin⁡πxcos⁡πt, in such a way that Tx,0=sin⁡πx, T0,t=T1,t=0. The analytical solution of this problem is Tx,t=sin⁡πxcos⁡πt. The graphs of the analytical and numerical solutions, when M=3, k=3,4,5, are shown in Figures 6–9.Figure 6
Analytical solution.Figure 7
Numerical solution withk=3.Figure 8
Numerical solution withk=4.Figure 9
Numerical solution withk=5.Example 8.
Consider (41), with α=2,1.9,1.8,1.7; the numerical solutions when k=M=4 at t=0.3,0.6,0.95 are shown in Figure 10. This example is introduced to verify the robustness of the proposed method; when the fractional order gradually approaches to 2, the numerical solutions are in agreement with the analytical solution.Figure 10
The numerical solutions withα=2,1.9,1.8,1.7 when k=M=4.


## 7. Conclusions

This paper presents a numerical technique for approximating solutions of one-dimensional fractional nonsteady heat transfer model by combining the second kind Chebyshev wavelet with its operational matrix of fractional-order integration. In the proposed method, a small number of grid points guarantee the necessary accuracy. The main advantage of wavelet method for solving the kinds of equations is that, after dispersing the coefficients, matrix of algebraic equations is sparse. The solution is convenient, even though the size of increment may be large. Several examples are given to demonstrate the powerfulness of the proposed method.



---

*Source: 2658124-2017-12-11.xml*