Crank-Nicolson Scheme for Space Fractional Heat
Conduction Equation with Mixed Boundary
Condition
Rafał Brociek
Institute of Mathematics
Silesian University of Technology
Kaszubska 23, 44-100 Gliwice, Poland
Email: rafal.brociek@polsl.pl
Abstract—This paper presents Crank-Nicolson scheme for numerical solution of space fractional diffusion equation with
space fractional heat conduction equation, formulated with boundary condition of the first kind, and in paper [12] au-
Riemann-Liouville fractional derivative. Dirichlet and Robin thors presents finite difference method for two-dimensional
boundary condition will be considered. To illustrate the accu- fractional dispersion equation. In both papers, as the fractional
racy of described method some computational examples will be derivative, the Riemann-Liouville derivative was used.
presented as well.
In this paper we present numerical solution of space frac-
I. I NTRODUCTION tional heat conduction equation. To the equation the Dirichlet
and Robin boundary condition was added. In order to solve the
In recent years the applications of mathematical models equation, the Crank-Nicolson scheme was used. To illustrate
using the fractional order derivatives are very popular in the accuracy of described method some computational exam-
technical science. Different types of phenomena in physics, ples will be presented as well. The aim to achieve numerical
biology, viscoelasticity, heat transfer, electrical engineering, solution of considered model is application it to solve inverse
control theory, fluid and continuum mechanics can be modeled problems. Numerical solution of described model is called
by using the fractional order derivatives [6], [23], [10], solution of direct problem. In the process of solving inverse
[19], [33]. For example in paper [20] mathematical models of problem, it is required to solve the direct problem many times.
supercapacitors are considered. Authors investigated models More about inverse problems of fractional order can be found
based on electrical circuits in the form of RC ladder networks. in papers [2], [3], [4], [5]. In these papers, swarm optimization
These models are described using fractional order differential algorithms are introduced. More about intelligent algorithms
equations. Often we are not able to solve these models in and its application can be found in [24], [25], [26], [27], [28].
an analytical way, so it is important to develop approximate
methods of solving differential equations of fractional order.
Murio in paper [21] presents the implicit finite differ- II. F ORMULATION OF THE P ROBLEM
ence approximation for the time fractional diffusion equations
We discuss the following space fractional heat conduction
with homogeneous Dirichlet boundary conditions, formulated
equation
by Caputo derivative. In paper [1], implicit finite difference
method was used to solve time fractional heat conduction equa- ∂u(x, t) ∂ α u(x, t)
tion with Neumann and Robin boundary conditions. In these c% = λ(x) + g(x, t), (1)
case also, the Caputo derivative was used. In paper [30] authors ∂t ∂xα
describes approximated solution of the fractional equation with defined in area D = {(x, t) : x ∈ [a, b], t ∈ [0, T )}, where
Dirichlet and Neumann boundary conditions. λ(x) is thermal conductivity coefficient, c and % denote the
Paper [7] describes numerical method for diffusion equa- specific heat and density. We assume that λ(x) > 0 and α ∈
tion with Dirichlet boundary conditions. To discretization (1, 2]. The initial condition is also posed
fractional derivative authors used finite difference method
and Kansa method. Paper [11] presents numerical solution u(x, 0) = f (x), x ∈ [a, b], (2)
of model describe by fractional differential equation with
space fractional derivative and Dirichlet boundary conditions. as well as the homogeneous Dirichlet and Robin boundary
Fractional derivative used in this model was Riemann-Liouville conditions
derivative. Authors used finite volume method.
u(a, t) = q(t), t ∈ [0, T ), (3)
Also Meerschaert and Tadjeran dealt with fractional dif- ∂u(b, t)
ferential equations [13], [12], [14], [35]. Paper [35] describes −λ(x) = h(t)(u(b, t) − u∞ ), t ∈ [0, T ), (4)
∂x
Copyright c 2016 held by the author. where h describes the heat transfer coefficient and u∞ is the
ambient temperature.
41
Fractional derivative with respect to space, which occurs in
equation (1), will be the Riemann-Liouville fractional deriva- 1
tive [23], [8] determined as follows (IN −1 −A1 )U k+1 = (IN −1 +A2 )U k +Gk+ 2 ∆t, k = 0, 1, 2, . . . ,
(11)
Zx
∂ α u(x, t) 1 ∂n where
= u(s, t)(x − s)n−1−α ds, (5)
∂tα Γ(n − α) ∂xn
a
U k = [U1k , U2k , . . . , UN
k T
−1 ] ,
where α ∈ (n − 1, n] and Γ(·) is the Gamma function [34]. In
case of α ∈ (1, 2), the equation (1) describes super-diffusive h g k+ 12 ∆t k+ 1
process [15], [16], and for α = 2, we get the classical heat k+ 12 g 2 ∆t
G ∆t = 1
, . . . , N −2 ,
conduction equation. c% c%
k+ 21
III. N UMERICAL S OLUTION gN −1 ∆t λN −1 ∆tu∞ ∆xhk ∆xhk+1 i
+ + ,
c% 2(∆x)α λN + ∆xhk λN + ∆xhk+1
In this section we describe numerical solution of equation
(1) using the finite difference method. Let N, M ∈ N be
IN −1 is the identity matrix of size N −1×N −1, and matrices
the size of grid in space and time, respectively. We denote
A1 and A2 are defined as follows
grid steps ∆x = (b−a)N and ∆t = T /M . Therefore, we get
following grid
i = 1, 2, . . . , N − 1, j = 1, 2, . . . , N − 1,
λ ∆t
2c%(∆x)α ωα,i−j+1
i
j ≤ i − 1,
S = (xi , tk ), xi = i ∆x, tk = k ∆t, i = 0, 1, . . . , N,
(6)
λi ∆t ω
j = i ∧ i 6= N − 1,
k = 0, 1, 2, . . . , M .
2c%(∆x)
α α,1
λi ∆t λN
a1ij = 2c%(∆x)α (ωα,1 + λN +∆xhk+1 ) j = i = N − 1,
We assume the following notation λi = λ(xi ), gik =
g(xi , tk ), fi = f (xi ), hk = h(tk ). Values of approximate λi ∆t
2c%(∆x)α ωα,0 j = i + 1,
function in points (xi , tk ), we denoted by Uik .
0 j > i + 1,
In order to approximate fractional derivative (5), we used
right-shifted Grünwald formula [17]
λ ∆t
2c%(∆x)α ωα,i−j+1
i
j ≤ i − 1,
α N
∂ u(x, t) 1 1 Γ(j − α)
X λi ∆t
= lim α u(x−(j −1)r, t),
α ωα,1 j = i ∧ i 6= N − 1,
2c%(∆x)
α
∂x Γ(−α) N →∞ r j=0 Γ(j + 1)
(7) a2ij = λi ∆t λN
2c%(∆x)α (ωα,1 + λN +∆xhk ) j = i = N − 1,
where N ∈ N and r = x−a λi ∆t α ωα,0 j = i + 1,
N . We denote
2c%(∆x)
0 j > i + 1,
Γ(j − α)
ωα,j = . (8)
Γ(−α)Γ(j + 1) Solving systems of equations defined by (11), we obtain
approximate values of temperatures in points of grid (6).
Discretize the equation (1) and using approximation of
fractional derivative, we get In paper [35] authors presents proof of unconditionally
stability of presented method in case of homogeneous Dirichlet
boundary condition. Doing similar proof it can be proven that
Uik+1 − Uik λi
i+1
X
k+1
i+1
X
k
described in this paper method is unconditionally stable.
= ωα,j U i−j+1 + ωα,j Ui−j+1
∆t 2c%(∆x)α j=0 j=0 IV. E XPERIMENTAL RESULTS
k+ 1
gi 2 In this section we presents numerical examples of described
+ .
c% method.
(9)
Example 1: Let consider equation (1), defined in area
Now, we approximate boundary condition of the third kind,
so we obtain D = {(x, t) : x ∈ [0, 1], t ∈ [0, 1]},
with following data
k k ∞
k λN UN −1 + ∆xh u
UN = . (10) 1
λN + ∆xhk α = 1.7, λ(x) = Γ(2.2)x2.8 , c = % = 1, u∞ = 50,
6
Having regard to both boundary condition, the equation (9) 0.550901e−t
may be written in matrix form g(x, t) = −(1 + x)x3 e−t , h(t) = − .
e−t − 50
42
To equation we added an initial condition Figure 2 presents exact, approximate solutions and errors
for time t = 1.
u(x, 0) = x3 , x ∈ [0, 1]. a) b)
0.006
0.35
0.005
Exact solution of this problem is function 0.30
0.25
0.004
uH x,1 L
error
0.20
0.003
3 −t 0.15
u(x, t) = x e . 0.002
0.10
0.001
0.05
0.000 0.00
Maximal and average absolute errors will be defined by 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x x
formulas
Fig. 2. Distribution of errors for approximate solution (a) and approximate
solution (points), exact solution (solid line) (b) in time t = 1 (N = M = 100)
∆max = max |uki − Uik |, (example 1)
0≤i≤N
1≤k≤M
Example 2: Let consider again equation (1) with following
N XM
1 X data
∆avg = |uki − Uik |.
(N + 1)(M + 1) i=0 1
k=0 α = 1.9, λ(x) = , c = % = 1, u∞ = 100, a = 1, b = 2,
2
Table I shows errors of approximate solution for different grids. 32.7273 − 32.7273x
g(x, t) = −2t(x−1)+0.0525569(t2 −1) +
TABLE I. G RIDS , MAXIMAL ABSOLUTE ERRORS ∆max AND AVERAGE (x − 1)1.9
ABSOLUTE ERRORS ∆avg ( EXAMPLE 1)
18.1818 15.5455 − 31.0909x + 15.5455x2
+ + ,
Grid N × M ∆max ∆avg (x − 1)0.9 (x − 1)2.9
10 × 10 6.24590 · 10−2 1.23295 · 10−2
6.20448 · 10−2 1.26218 · 10−2 1
10 × 50 (1 − t2 )
20 × 20 3.09731 · 10−2 5.58958 · 10−3 h(t) = 2 .
50 × 50 1.23473 · 10−2 2.08623 · 10−3 99 + t2
100 × 100 6.16827 · 10−3 1.01756 · 10−3
100 × 200 6.16490 · 10−3 1.01786 · 10−3
and with initial condition
100 × 300 6.16429 · 10−3 1.01808 · 10−3
200 × 100 3.08626 · 10−3 5.02667 · 10−4 u(x, 0) = x − 1, t ∈ [0, 1].
300 × 100 2.05917 · 10−3 3.34012 · 10−4
Exact solution of this problem is given by function
With increase in the first dimension N of the grid, a
fixed second dimension M , errors decrease. For example, for
M = 100 and N = 100, 200, 300 absolute errors not exceed u(x, t) = (x − 1)(1 − t2 ).
6.17 · 10−3 , 3.09 · 10−3 , 2.06 · 10−3 . Increasing the dimension
of the grid with respect to time for a fixed dimension N Table II presents errors of approximate solution for different
insignificantly impact on errors of approximate solution. grids.
Distribution of errors in points of the grid is presented on TABLE II. G RIDS , MAXIMAL ABSOLUTE ERRORS ∆max AND
figure 1. AVERAGE ABSOLUTE ERRORS ∆avg ( EXAMPLE 2)
Grid N × M ∆max ∆avg
10 × 10 7.14165 · 10−4 4.40379 · 10−4
10 × 50 6.25414 · 10−4 3.82786 · 10−4
20 × 20 4.16933 · 10−4 2.54149 · 10−4
50 × 50 2.11043 · 10−4 1.22171 · 10−4
100 × 100 1.2352 · 10−4 6.88743 · 10−5
100 × 200 1.23129 · 10−4 6.84358 · 10−5
100 × 300 1.23059 · 10−4 6.83813 · 10−5
200 × 100 7.64731 · 10−5 3.8691 · 10−5
300 × 100 5.90598 · 10−5 2.74809 · 10−5
Similarly as in example 1, reducing the grid step ∆x results
in a significant reduction in average errors δavg and maximum
errors ∆max of approximate solution. Also, the density of
the grid in the time domain results in a slight decrease in
errors. Distribution of errors in the area of D, for example 2
is presented in the figure 3.
We also present distribution of errors for approximate
solution in time t = 1 (figure 4).
Fig. 1. Distribution of errors (N = M = 100) (example 1) Example 3: Now, we take the following data
43
V. C ONCLUSION
In this paper numerical solution of space heat conduction
equation is presented. To the equation Robin and Neumann
boundary conditions were added. Author used Crank-Nicolson
scheme, and right-shifted Grünwald formula for approximation
fractional Riemann-Liouville derivative. To illustrate the accu-
racy of presented method three examples are presented. In each
example results are good. With an increase in the density of
the grid, errors of approximate solution decreases.
R EFERENCES
[1] R. Brociek: Implicite finite difference method for time fractional heat
equation with mixed boundary conditions. Zesz. Nauk. PŚ., Mat.
Stosow. 4 (2014), 72–85.
[2] R. Brociek, D. Słota, R. Wituła: Reconstruction of the thermal conduc-
tivity coefficient in the time fractional diffusion equation, in: Advances
in Modelling and Control of Non-integer-Order Systems, K.J. Latawiec,
Fig. 3. Distribution of errors (N = M = 100) (example 2) M. Łukaniszyn, R. Stanisławski, eds., Lecture Notes in Electrical
Engineering, 320, Springer, 239–247 (2015).
a) b) [3] R. Brociek, D. Słota: Reconstruction of the boundary condition for
0.00008
0 the heat conduction equation of fractional order, Thermal Science 19,
-0.00002
35–42 (2015).
0.00006
[4] R. Brociek, D. Słota: Application of intelligent algorithm to solve
uH x,1 L
error
0.00004
-0.00004
the fractional heat conduction inverse problem, in: Information and
-0.00006
Software Technologies, G. Dregvaite, R. Damasevicius, eds., Commu-
0.00002
nications in Computer and Information Science, 538, Springer, 356–366
0
-0.00008
(2015).
1.0 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 1.6 1.8 2.0
x x [5] R. Brociek, D. Słota: Application and comparison of intelligent
algorithms to solve the fractional heat conduction inverse problem. Inf.
Fig. 4. Distribution of errors for approximate solution (a) and approximate Technol. Control T. 45 nr 2 (2016), 184–194.
solution (points), exact solution (solid line) (b) in time t = 1 (N = M = 100)
(example 2) [6] A. Carpinteri, F. Mainardi: Fractal and Fractional Calculus in Contin-
uum Mechanics. Springer, New York (1997)
[7] W. Chen, L. Ye, H. Sun: Fractional diffusion equations by the kansa
method. Comput. Math. Appl. 59 (2010), 1614–1620.
[8] K. Diethelm: The Analysis of Fractional Differential Equations.
α = 1.6, λ(x) = 1, c = % = 1, u∞ = 120, a = 1, b = 2, Springer, Berlin (2010)
[9] G.H. Gao, Z.A. Sun: A compact finite difference scheme for the
5.15228 fractional sub-diffusion equations. J. Comput. Phys. 230 (2011), 586–
g(x, t) = − +
(x − 1)0.6 595.
[10] J. Klafter, S. Lim, R. Metzler: Fractional dynamics. Resent advances.
t3 (0.0751374 + 1.87843(x − 1) + 0.375687x) World Scientific, New Jersey (2012)
+ +3xt2 (1−x),
(x − 1)0.6 [11] F. Liu, P. Zhuang, I. Turner, K. Burrage, V. Anh: A new fractional
finite volume method for solving the fractional diffusion equation. Appl.
Math. Modelling 38 (2014), 3871–3878.
8 − 3t3
h(t) = , f (x) = 8(x − 1). [12] M.M. Meerschaert, H.P. Scheffler, Ch. Tadjeran: Finite difference
112 + 2t3 methods for two-dimensional fractional dispersion equation. J. Comput.
Phys. 211 (2006), 249–261.
Exact solution of this problem is represented by function [13] M.M. Meerschaert, Ch. Tadjeran: Finite difference approximations for
fractional advection-dispersion flow equations. J. Comput. Appl. Math.
172 (2006), 65–77.
1
u(x, t) = 2(x − 1)(4 − xt3 ). [14] M.M. Meerschaert, Ch. Tadjeran: Finite difference approximations for
2 two-sided space-fractional partial differential equations. Appl. Numer.
Math. 56 (2006), 80–90.
Just as it did in the previous examples, we will examine [15] M.M. Meerschaert, D.A. Benson, H.P. Scheffler, P. Becker-Kern: Gov-
the errors depending on the density of the grid. In table III erning equations and solutions of anomalous random walk limits. Phys.
errors of approximate solutions are presented. Rev. E 66 (2002), 102R–105R.
[16] R. Metzler, J. Klafter: The restaurant at the end of the random
TABLE III. G RIDS , MAXIMAL ABSOLUTE ERRORS ∆max AND walk: recent developments in the description of anomalous transport
AVERAGE ABSOLUTE ERRORS ∆avg ( EXAMPLE 3) by fractional dynamics. J. Phys. Rev. A 37 (2004), R161–R208.
[17] K. Miller, B. Ross: An Introduction to the Fractional Calculus and
Grid N × M ∆max ∆avg
Fractional Differential. Wiley, New York, 1993.
10 × 10 1.2654 · 10−1 4.8953 · 10−2
50 × 50 2.48109 · 10−2 1.2146 · 10−2 [18] W. Mitkowski: Approximation of fractional diffusion-wave equation.
70 × 70 1.77869 · 10−2 8.93166 · 10−3 Acta Mechanica et Aautomatica 5, no. 2 (2011), 65–68.
100 × 100 1.25078 · 10−2 6.41806 · 10−3 [19] W. Mitkowski, J. Kacprzyk, J. Baranowski: Advances in the Theory
150 × 150 8.38441 · 10−3 4.38665 · 10−3 and Applications of Non-integer Order Systems. Springer Inter. Publ.,
Cham (2013)
44
[20] W. Mitkowski, P. Skruch: Fractional-order models of the supercapaci-
tors in the form of RC ladder networks. Bulletin of The Polish Academy
of Sciences Technical Sciences, Vol. 61, No. 3, 2013. pp. 581-587.
[21] D. Murio: Implicit finite difference approximation for time fractional
diffusion equations. Comput. Math. Appl. 56 (2008) 1138–1145
[22] A. Obra̧czka: The L2 approximation and shifted Grünwald methods
comparison for RFDE, in: Materiay XV Jubileuszowego Sympozjum
,,Podstawowe Problemy Energoelektroniki, Elektromechaniki i Mecha-
troniki”. PPEEm 2012, M. Szczygieł, ed., Archiwum Konferencji
PTETiS, 32, Komitet Organizacyjny Sympozjum PPEE i Seminarium
BSE, 2012, 129–132.
[23] I. Podlubny: Fractional Differential Equations. Academic Press, San
Diego (1999)
[24] D. Połap, M. Wozniak: Flexible Neural Network Architecture for Hand-
written Signatures Recognition. International Journal of Electronics and
Telecommunications, Vol. 62, No. 2, DOI: 10.1515/eletel-2016-0027,
ISSN 2300-1933, De Gruyter Open Ltd 2016, pp. 197-202.
[25] D. Połap, M. Wozniak, C. Napoli, E. Tramontana: Real-Time Cloud-
based Game Management System via Cuckoo Search Algorithm. Inter-
national Journal of Electronics and Telecommunications, Vol. 61, No.
4, DOI: 10.1515/eletel-2015-0043, ISSN 2300-1933, De Gruyter Open
Ltd 2015, pp. 333-338.
[26] D. Połap, M. Wozniak, C. Napoli, E. Tramontana: Is swarm intelli-
gence able to create mazes? International Journal of Electronics and
Telecommunications, Vol. 61, No. 4, DOI: 10.1515/eletel-2015-0039,
ISSN 2300-1933, De Gruyter Open Ltd 2015, pp. 305-310.
[27] C. Napoli, G. Pappalardo, E. Tramontana: A Mathematical Model For
File Fragment Diffusion And A Neural Predictor To Manage Priority
Queues Over Bittorrent. International Journal of Applied Mathematics
and Computer Science. Vol. 26, No. 1, Univ Zielona Gora Press 2016,
pp. 147-160
[28] C. Napoli, G. Pappalardo, E. Tramontana: Improving files availability
for BitTorrent using a diffusion model. Proceedings of IEEE Interna-
tional WETICE Conference June, 2014, pp. 191-196
[29] C. Napoli, E. Tramontana: Massively parallel WRNN reconstructors
for spectrum recovery in astronomical photometrical surveys. Neural
Networks. Vol. 83, Elsevier 2016, pp. 42-50
[30] G. Sudha Priya, P. Prakash, J.J. Nieto, Z. Kayar: Higher-order numerical
scheme for the fractional heat equation with Dirichlet and Neumann
boundary conditions. Numer. Heat Transfer B 63 (2013), 540–559.
[31] J. Ren, Z.Z. Sun, X. Zhao: Compact difference scheme for the fractional
sub-diffusion equation with Neumann boundary conditions. J. Comput.
Phys. 232 (2013), 456–467.
[32] J.P. Roop: Computational aspects of fem approximation of fractional
advection dispersion equations on bounded domains in R2 . J. Comput.
Appl. Math. 193 (2006), 243–268.
[33] J. Sabatier, O.P. Agrawal, J.A. Tenreiro Machado: Advances in Frac-
tional Calculus. Theoretical Developments and Applications in Physics
and Engineering. Springer, Dordrecht 2007.
[34] Sz. Rabsztyn, D. Słota, R. Wituła: Gamma and Beta Functions, vol. 1
(in Polish). Wyd. Pol. Śl., Gliwice (2012)
[35] Ch. Tadjeran, M.M. Meerschaert, H.P. Scheffler: A second-order
accurate numerical approximation for the fractional diffusion equation.
J. Comput. Phys. 213 (2006), 205–213.
[36] X. Zhao, Q. Xu: Efficient numerical schemes for fractional sub-
diffusion equation with the spatially variable coefficient. Appl. Math.
Modelling 38 (2014), 3848–3859.
[37] Y. Zheng, Ch. Li, Z. Zhao: A note on the finite element method for the
space-fractional advection diffusion equation. Comput. Math. Appl. 59
(2010), 1718–1726.
45