# Eighth-Order Compact Finite Difference Scheme for 1D Heat Conduction Equation

Authors: Asma Yosaf; Shafiq Ur Rehman; Fayyaz Ahmad; Malik Zaka Ullah; Ali Saleh Alshomrani
Journal: Advances in Numerical Analysis
 (2016)

---

## Abstract

The purpose of this paper is to develop a high-order compact finite difference method for solving one-dimensional (1D) heat conduction equation with Dirichlet and Neumann boundary conditions, respectively. A parameter is used for the direct implementation of Dirichlet and Neumann boundary conditions. The introduced parameter adjusts the position of the neighboring nodes very next to the boundary. In the case of Dirichlet boundary condition, we developed eighth-order compact finite difference method for the entire domain and fourth-order accurate proposal is presented for the Neumann boundary conditions. In the case of Dirichlet boundary conditions, the introduced parameter behaves like a free parameter and could take any value from its defined domain but for the Neumann boundary condition we obtained a particular value of the parameter. In both proposed compact finite difference methods, the order of accuracy is the same for all nodes. The time discretization is performed by using Crank-Nicholson finite difference method. The unconditional convergence of the proposed methods is presented. Finally, a set of 1D heat conduction equations is solved to show the validity and accuracy of our proposed methods.


---

## Body

## 1. Introduction

Heat conduction problems with suitable boundary conditions exist in many areas of engineering applications [1–7]. Historically, highly accurate compact finite difference schemes are developed in the work by Lele [8]. But these higher-order compact finite difference schemes only offer good accuracy at the interior nodes or for periodic boundary conditions. Usually, compact finite difference schemes have first- or second-order accuracy [9, 10]. The low-order of accuracy near boundary grid points affects the whole numerical results and it reduces the accuracy of overall numerical solution [11]. Some authors offer one-side finite difference approximations for the Dirichlet boundary condition [8, 12] but they cannot offer unconditional stability for the whole finite difference scheme. Recently, Dai et al. [11, 13–16] proposed a new idea to achieve higher-order accuracy with unconditional stability. Actually, authors introduced a new parameter θ that adjusts the location of nodes near the boundaries in symmetric way.


## 2. Higher-Order Compact Finite Difference Method

The 1D heat conduction equation can be written as(1a)
∂
T
x
,
t
∂
t
=
k
∂
2
T
x
,
t
∂
x
2
+
s
x
,
t
,
0
<
x
<
l
,
0
<
t
≤
t
f
(1b)
T
x
,
0
=
f
x
,
0
≤
x
≤
l
.Dirichlet boundary conditions are as follows:(2)
T
0
,
t
=
T
0
t
,
T
l
,
t
=
T
l
t
,
0
≤
t
≤
t
f
.Neumann boundary conditions are as follows:(3)
T
x
0
,
t
=
T
0
t
,
T
x
l
,
t
=
T
l
t
,
0
≤
t
≤
t
f
.Han and Dai [17] have proposed a compact finite difference method for the spatial discretization of (1a) that has eighth-order accuracy at interior nodes and sixth-order accuracy for boundary nodes. Actually, they use eighth-order accurate approximation for the second-order derivative developed in [8] and given as (4)
β
3
T
i
-
2
′
′
+
α
3
T
i
-
1
′
′
+
T
i
′
′
+
α
3
T
i
+
1
′
′
+
β
3
T
i
+
2
′
′
=
b
3
T
i
+
2
-
2
T
i
+
T
i
-
2
4
h
2
+
a
3
T
i
+
1
-
2
T
i
+
T
i
-
1
h
2
,where double dash always means derivative with respect to spatial variable and α
3
=
344
/
1179, β
3
=
38
α
3
-
9
/
214, a
3
=
696
-
1191
α
3
/
428, and b
3
=
2454
α
3
-
294
/
535. The implicit compact finite difference scheme (4) without the contribution of boundary nodes in matrix form is given by(5)
A
T
x
x
=
B
T
.Here,(6)
A
=
β
3
α
3
1
α
3
β
3
β
3
α
3
1
α
3
β
3
β
3
α
3
1
α
3
β
3
⋱
⋱
⋱
⋱
β
3
α
3
1
α
3
β
3
β
3
α
3
1
α
3
β
3
B
=
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
⋱
⋱
⋱
⋱
⋱
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
.The governing equation of implicit compact finite difference approximation of second-order derivative is(7)
β
T
i
-
2
′
′
+
α
T
i
-
1
′
′
+
T
i
′
′
+
α
T
i
+
1
′
′
+
β
T
i
+
2
′
′
=
c
T
i
+
3
-
2
T
i
+
T
i
-
3
9
h
2
+
b
T
i
+
2
-
2
T
i
+
T
i
-
2
4
h
2
+
a
T
i
+
1
-
2
T
i
+
T
i
-
1
h
2
.Note that (4) can be obtained from (7).We can observe that matrixA is strictly diagonally dominant and matrix B is diagonally dominant. In [17], authors constructed finite difference scheme at boundary nodes in such a way that they conserve the diagonal dominance of matrixes A and B to attain higher order of accuracy. The achieved order of accuracy at boundary nodes in [17] is six. The essence of article [17] is hidden in the calculation of a parameter θ that they use to obtain higher-order approximation near boundary nodes. The parameter θ makes the point distribution unequal near boundary but offer higher order of accuracy and diagonal dominance for matrices A and B. The diagonally dominance of A and B is the key to prove stability of whole compact finite difference scheme when they use Crank-Nicholson method for time integration.We construct a new finite difference scheme for boundary nodes to achieve eighth order of accuracy that exactly matches the order of accuracy at the interior nodes. The inclusion of two parametersθ and α’s will be introduced in our newly developed finite difference scheme for boundary nodes. The diagonally dominance will be conserved and the whole eighth-order scheme becomes stable.The stencil of eighth-order implicit finite difference scheme to approximate the second-order derivative for the interior nodes is(8)
i
-
2
,
i
-
1
,
i
,
i
+
1
,
i
+
2
for
i
=
3,4
,
…
,
M
-
2
.It means that if we are at location i, then we need two grid nodes to the left of i and two grid points to the right of it. It is noticeable that mutual distance between nodes is equal to h
=
l
/
2
θ
+
M
-
2 and M points divide the interval [
0
,
l
] with x
0
=
0, x
M
+
2
=
l, x
1
-
x
0
=
θ
h, x
M
+
1
-
x
M
=
θ
h, and x
i
=
(
θ
+
i
-
1
)
h for i
=
2,3
,
…
,
M
-
1. Figure 1 shows the location of interior and boundary nodes.Figure 1
Grid points and location of boundary nodes.In the case of Dirichlet boundary conditions, we construct the eighth-order accurate implicit finite difference approximations of second-order derivatives for the nodesx
1, x
2, x
M, and x
M
+
1, whereas, in the case of Neumann boundary conditions, we establish stable fourth-order compact finite difference scheme for the entire grid.


## 3. Compact Implicit Finite Difference Scheme for Boundary Nodes in the Case of Dirichlet Boundary Conditions

In the case of Dirichlet boundary condition, we construct the following compact implicit finite difference scheme at pointx
1
=
θ
h:(9)
α
1
T
x
x
x
0
+
T
x
x
x
1
=
1
h
2
a
0
T
x
0
+
a
1
T
x
1
+
a
2
T
x
2
+
⋯
+
a
8
T
x
8
,where θ, α
1 (≠0), and a
i for i
=
0,1
,
2
,
…
,
8 are unknowns and we are interested to find them in such a way that we can get(i)
eighth-order accurate approximation of second-order derivative atx
1,
(ii)
diagonal dominance, that is,a
1
<
0 and a
1
≥
∑
j
=
1
,
j
≠
1
8
a
j.By expanding (9) around the node x
1
=
θ
h and equating the coefficients of the same-order derivatives, we get the following system of equations:(10)
-
a
0
-
a
1
-
a
2
-
a
3
-
a
4
-
a
5
-
a
6
-
a
7
-
a
8
=
0
a
0
θ
-
a
2
-
2
a
3
-
3
a
4
-
4
a
5
-
5
a
6
-
6
a
7
-
7
a
8
=
0
-
1
2
a
0
θ
2
-
a
2
2
-
2
a
3
-
9
2
a
4
-
8
a
5
-
25
2
a
6
-
18
a
7
-
49
a
8
2
+
α
1
+
1
=
0
1
6
a
0
θ
3
-
α
1
θ
-
a
2
6
-
4
3
a
3
-
9
2
a
4
-
32
a
5
3
-
125
a
6
6
-
36
a
7
-
343
a
8
6
=
0
-
1
24
a
0
θ
4
+
1
2
α
1
θ
2
-
a
2
24
-
2
3
a
3
-
27
a
4
8
-
32
a
5
3
-
625
a
6
24
-
54
a
7
-
2401
a
8
24
=
0
a
0
θ
5
120
-
1
6
α
1
θ
3
-
a
2
120
-
4
a
3
15
-
81
a
4
40
-
128
a
5
15
-
625
a
6
24
-
324
a
7
5
-
16807
a
8
120
=
0
-
a
0
θ
6
720
+
1
24
α
1
θ
4
-
a
2
720
-
4
a
3
45
-
81
a
4
80
-
256
a
5
45
-
3125
a
6
144
-
324
a
7
5
-
117649
a
8
720
=
0
a
0
θ
7
5040
-
α
1
θ
5
120
-
a
2
5040
-
8
a
3
315
-
243
a
4
560
-
1024
a
5
315
-
15625
a
6
1008
-
1944
a
7
35
-
117649
a
8
720
=
0
-
a
0
θ
8
40320
+
α
1
θ
6
720
-
a
2
40320
-
2
a
3
315
-
729
a
4
4480
-
512
a
5
315
-
78125
a
6
8064
-
1458
a
7
35
-
823543
a
8
5760
=
0
a
0
θ
9
362880
-
α
1
θ
7
5040
-
a
2
362880
-
4
a
3
2835
-
243
a
4
4480
-
2048
a
5
2835
-
390625
a
6
72576
-
972
a
7
35
-
5764801
a
8
51840
=
0
.By solving (10), we get the solution as follows:(11)
a
0
=
36
6534
θ
6
+
140336
θ
5
+
1186143
θ
4
+
4998000
θ
3
+
10913735
θ
2
+
11457768
θ
+
4357806
2
θ
8
+
63
θ
7
+
826
θ
6
+
5831
θ
5
+
23919
θ
4
+
57232
θ
3
+
75477
θ
2
+
46998
θ
+
8820
θ
+
6
θ
+
5
θ
+
4
θ
+
3
θ
+
2
θ
+
1
a
1
=
26264
θ
7
+
594463
θ
6
+
5162388
θ
5
+
20813366
θ
4
+
33289228
θ
3
-
15584793
θ
2
-
103119912
θ
-
78440508
5040
θ
7
+
123480
θ
6
+
1217160
θ
5
+
6174000
θ
4
+
17057880
θ
3
+
24819480
θ
2
+
16465680
θ
+
3175200
a
2
=
-
1
40
θ
1784
θ
7
+
40047
θ
6
+
345570
θ
5
+
1402175
θ
4
+
2454868
θ
3
+
368812
θ
2
-
3371760
θ
-
2116800
θ
+
1
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
a
3
=
3516
θ
7
+
76839
θ
6
+
641010
θ
5
+
2496750
θ
4
+
4223382
θ
3
+
1178139
θ
2
-
3322620
θ
-
1587600
θ
40
θ
+
80
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
a
4
=
-
θ
7592
θ
7
+
159775
θ
6
+
1269918
θ
5
+
4643015
θ
4
+
7205332
θ
3
+
1657932
θ
2
-
4782960
θ
-
2116800
72
θ
+
216
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
a
5
=
1
8
656
θ
7
+
13251
θ
6
+
100272
θ
5
+
345950
θ
4
+
502120
θ
3
+
102667
θ
2
-
309960
θ
-
132300
θ
θ
+
4
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
a
6
=
-
3
θ
536
θ
7
+
10367
θ
6
+
74842
θ
5
+
246119
θ
4
+
341492
θ
3
+
64588
θ
2
-
202608
θ
-
84672
40
θ
+
200
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
a
7
=
θ
4076
θ
7
+
75313
θ
6
+
520770
θ
5
+
1650050
θ
4
+
2218462
θ
3
+
397053
θ
2
-
1283940
θ
-
529200
360
θ
+
2160
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
a
8
=
-
θ
392
θ
7
+
6903
θ
6
+
45990
θ
5
+
141575
θ
4
+
186004
θ
3
+
31948
θ
2
-
105840
θ
-
43200
560
θ
8
+
17640
θ
7
+
231280
θ
6
+
1632680
θ
5
+
6697320
θ
4
+
16024960
θ
3
+
21133560
θ
2
+
13159440
θ
+
2469600
α
1
=
9
363
θ
-
140
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
.We define the error equation as(12)
E
1
=
α
1
T
x
x
x
0
+
T
x
x
x
1
-
1
h
2
a
0
T
x
0
+
a
1
T
x
1
+
a
2
T
x
2
+
⋯
+
a
8
T
x
8
.After substituting the solution (11) into (12), we get(13)
E
1
=
q
1
q
2
h
8
+
O
h
9
,where (14)
q
1
=
θ
862488
θ
8
+
20838048
θ
7
+
200482128
θ
6
+
972269760
θ
5
+
2446576440
θ
4
+
2782223136
θ
3
+
359486688
θ
2
-
1448979840
θ
-
558835200
f
10
q
2
=
53222400
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
.Similarly, at the node x
2, we have the following construction of the implicit compact finite difference scheme:(15)
α
1
T
′
′
x
0
+
α
2
T
′
′
x
1
+
T
′
′
x
2
=
1
h
2
a
0
T
x
0
+
a
1
T
x
1
+
a
2
T
x
2
+
⋯
+
a
8
T
x
8
,where θ, α
1 (≠0), α
2 (≠0), and a
i for i
=
0,1
,
2
,
…
,
8 are unknowns and we are interested to find them in such a way that we can get(i)
eighth-order accurate approximation of second-order derivative atx
2,
(ii)
diagonal dominance, that is,a
2
<
0 and a
2
≥
∑
j
=
1
,
j
≠
2
8
‍
a
j.By expanding (15) around the node x
2
=
(
θ
+
1
)
h and equating the coefficients of the same-order derivative, we get the system of the following:(16)
-
a
0
-
a
1
-
a
2
-
a
3
-
a
4
-
a
5
-
a
6
-
a
7
-
a
8
=
0
a
0
θ
+
a
0
+
a
1
-
a
3
-
2
a
4
-
3
a
5
-
4
a
6
-
5
a
7
-
6
a
8
=
0
-
1
2
a
0
θ
2
-
a
0
θ
-
a
0
2
-
a
1
2
-
a
3
2
-
2
a
4
-
9
2
a
5
-
8
a
6
-
25
2
a
7
-
18
a
8
+
α
1
+
α
2
+
1
=
0
1
6
a
0
θ
3
+
1
2
a
0
θ
2
+
1
2
a
0
θ
-
α
1
θ
+
a
0
6
+
a
1
6
-
a
3
6
-
4
3
a
4
-
9
2
a
5
-
32
a
6
3
-
125
a
7
6
-
36
a
8
-
α
1
-
α
2
=
0
-
1
24
a
0
θ
4
-
1
6
a
0
θ
3
-
1
4
a
0
θ
2
+
1
2
α
1
θ
2
-
1
6
a
0
θ
+
α
1
θ
-
a
0
24
-
a
1
24
-
a
3
24
-
2
3
a
4
-
27
a
5
8
-
32
a
6
3
-
625
a
7
24
-
54
a
8
+
α
1
2
+
α
2
2
=
0
a
0
θ
5
120
+
1
24
a
0
θ
4
+
1
12
a
0
θ
3
-
1
6
α
1
θ
3
+
1
12
a
0
θ
2
-
1
2
α
1
θ
2
+
1
24
a
0
θ
-
1
2
α
1
θ
+
a
0
120
+
a
1
120
-
a
3
120
-
4
a
4
15
-
81
a
5
40
-
128
a
6
15
-
625
a
7
24
-
324
a
8
5
-
α
1
6
-
α
2
6
=
0
-
a
0
θ
6
720
-
a
0
θ
5
120
-
1
48
a
0
θ
4
+
1
24
α
1
θ
4
-
1
36
a
0
θ
3
+
1
6
α
1
θ
3
-
1
48
a
0
θ
2
+
1
4
α
1
θ
2
-
a
0
θ
120
+
1
6
α
1
θ
-
a
0
720
-
a
1
720
-
a
3
720
-
4
a
4
45
-
81
a
5
80
-
256
a
6
45
-
3125
a
7
144
-
324
a
8
5
+
α
1
24
+
α
2
24
=
0
a
0
θ
7
5040
+
a
0
θ
6
720
+
a
0
θ
5
240
-
α
1
θ
5
120
+
a
0
θ
4
144
-
1
24
α
1
θ
4
+
a
0
θ
3
144
-
1
12
α
1
θ
3
+
a
0
θ
2
240
-
1
12
α
1
θ
2
+
a
0
θ
720
-
1
24
α
1
θ
+
a
0
5040
+
a
1
5040
-
a
3
5040
-
8
a
4
315
-
243
a
5
560
-
1024
a
6
315
-
15625
a
7
1008
-
1944
a
8
35
-
α
1
120
-
α
2
120
=
0
-
a
0
θ
8
40320
-
a
0
θ
7
5040
-
a
0
θ
6
1440
+
α
1
θ
6
720
-
a
0
θ
5
720
+
α
1
θ
5
120
-
a
0
θ
4
576
+
1
48
α
1
θ
4
-
a
0
θ
3
720
+
1
36
α
1
θ
3
-
a
0
θ
2
1440
+
1
48
α
1
θ
2
-
a
0
θ
5040
+
α
1
θ
120
-
a
0
40320
-
a
1
40320
-
a
3
40320
-
2
a
4
315
-
729
a
5
4480
-
512
a
6
315
-
78125
a
7
8064
-
1458
a
8
35
+
α
1
720
+
α
2
720
=
0
-
α
1
θ
7
5040
+
a
0
θ
9
362880
+
a
0
θ
8
40320
-
α
1
θ
6
720
+
a
0
θ
7
10080
-
α
1
θ
5
240
+
a
0
θ
6
4320
-
α
1
θ
4
144
-
α
2
5040
-
α
1
5040
-
972
a
8
35
-
2048
a
6
2835
-
243
a
5
4480
-
4
a
4
2835
-
a
3
362880
+
a
1
362880
+
a
0
362880
-
α
1
θ
3
144
+
a
0
θ
2
10080
-
α
1
θ
2
240
+
a
0
θ
40320
-
α
1
θ
720
+
a
0
θ
5
2880
+
a
0
θ
4
2880
+
a
0
θ
3
4320
-
390625
a
7
72576
=
0
.By solving (16), we get the following solution: (17)
a
0
=
36
6534
α
2
θ
7
+
140336
α
2
θ
6
-
522
θ
7
+
1186143
α
2
θ
5
-
11494
θ
6
+
4998000
α
2
θ
4
-
100695
θ
5
+
10913735
α
2
θ
3
-
448035
θ
4
+
11457768
α
2
θ
2
-
1069705
θ
3
+
4357806
α
2
θ
-
1325247
θ
2
-
745740
θ
-
131886
θ
2
θ
8
+
63
θ
7
+
826
θ
6
+
5831
θ
5
+
23919
θ
4
+
57232
θ
3
+
75477
θ
2
+
46998
θ
+
8820
θ
+
1
θ
+
2
θ
+
3
θ
+
4
θ
+
5
θ
+
6
-
1
a
1
=
26264
α
2
θ
8
+
594463
α
2
θ
7
+
3528
θ
8
+
5162388
α
2
θ
6
+
90351
θ
7
+
20813366
α
2
θ
5
+
947583
θ
6
+
33289228
α
2
θ
4
+
5259870
θ
5
-
15584793
α
2
θ
3
+
16674966
θ
4
-
103119912
α
2
θ
2
+
30375639
θ
3
-
78440508
α
2
θ
+
30308103
θ
2
+
14582268
θ
+
2373948
2520
θ
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
-
1
a
2
=
-
16056
α
2
θ
8
+
360423
α
2
θ
7
+
280
θ
8
+
3110130
α
2
θ
6
+
11055
θ
7
+
12619575
α
2
θ
5
+
168485
θ
6
+
22093812
α
2
θ
4
+
1312969
θ
5
+
3319308
α
2
θ
3
+
5706195
θ
4
-
30345840
α
2
θ
2
+
13929248
θ
3
-
19051200
α
2
θ
+
18006864
θ
2
+
10528632
θ
+
1900080
360
θ
+
360
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
-
1
a
3
=
3516
α
2
θ
8
+
76839
α
2
θ
7
-
216
θ
8
+
641010
α
2
θ
6
-
4419
θ
7
+
2496750
α
2
θ
5
-
31935
θ
6
+
4223382
α
2
θ
4
-
79182
θ
5
+
1178139
α
2
θ
3
+
122958
θ
4
-
3322620
α
2
θ
2
+
1008645
θ
3
-
1587600
α
2
θ
+
1830321
θ
2
+
1230444
θ
+
234360
40
θ
+
80
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
-
1
a
4
=
-
7592
α
2
θ
8
+
159775
α
2
θ
7
-
684
θ
8
+
1269918
α
2
θ
6
-
14895
θ
7
+
4643015
α
2
θ
5
-
124587
θ
6
+
7205332
α
2
θ
4
-
494721
θ
5
+
1657932
α
2
θ
3
-
898227
θ
4
-
4782960
α
2
θ
2
-
414864
θ
3
-
2116800
α
2
θ
+
680382
θ
2
+
740124
θ
+
158760
72
θ
+
216
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
-
1
a
5
=
5904
α
2
θ
8
+
119259
α
2
θ
7
-
536
θ
8
+
902448
α
2
θ
6
-
11361
θ
7
+
3113550
α
2
θ
5
-
92665
θ
6
+
4519080
α
2
θ
4
-
363530
θ
5
+
924003
α
2
θ
3
-
688698
θ
4
-
2789640
α
2
θ
2
-
494185
θ
3
-
1190700
α
2
θ
+
132495
θ
2
+
286452
θ
+
66780
72
θ
+
288
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
-
1
a
6
=
-
1608
α
2
θ
8
+
31101
α
2
θ
7
-
144
θ
8
+
224526
α
2
θ
6
-
2943
θ
7
+
738357
α
2
θ
5
-
23169
θ
6
+
1024476
α
2
θ
4
-
88425
θ
5
+
193764
α
2
θ
3
-
166923
θ
4
-
607824
α
2
θ
2
-
130872
θ
3
-
254016
α
2
θ
+
4716
θ
2
+
48816
θ
+
12096
40
θ
+
200
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
-
1
a
7
=
4076
α
2
θ
8
+
75313
α
2
θ
7
-
360
θ
8
+
520770
α
2
θ
6
-
7065
θ
7
+
1650050
α
2
θ
5
-
53685
θ
6
+
2218462
α
2
θ
4
-
199746
θ
5
+
397053
α
2
θ
3
-
373590
θ
4
-
1283940
α
2
θ
2
-
303057
θ
3
-
529200
α
2
θ
-
17541
θ
2
+
87372
θ
+
22680
360
θ
+
2160
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
-
1
a
8
=
-
3528
α
2
θ
8
+
62127
α
2
θ
7
-
308
θ
8
+
413910
α
2
θ
6
-
5787
θ
7
+
1274175
α
2
θ
5
-
42595
θ
6
+
1674036
α
2
θ
4
-
155141
θ
5
+
287532
α
2
θ
3
-
287511
θ
4
-
952560
α
2
θ
2
-
237160
θ
3
-
388800
α
2
θ
-
25662
θ
2
+
57492
θ
+
15480
5040
θ
8
+
158760
θ
7
+
2081520
θ
6
+
14694120
θ
5
+
60275880
θ
4
+
144224640
θ
3
+
190202040
θ
2
+
118434960
θ
+
22226400
-
1
α
1
=
9
363
α
2
θ
-
140
α
2
-
29
θ
-
9
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
-
1
.The error equation at the nodex
2 is given by(18)
E
2
=
α
1
T
x
x
x
0
+
α
2
T
x
x
x
1
+
T
x
x
x
2
-
1
h
2
a
0
T
x
0
+
a
1
T
x
1
+
a
2
T
x
2
+
⋯
+
a
8
T
x
8
.By substituting (17) into (18), we get(19)
E
2
=
q
1
q
2
h
8
+
O
h
9
,where (20)
q
1
=
7762392
α
2
θ
9
+
187542432
α
2
θ
8
-
620136
θ
9
+
1804339152
α
2
θ
7
-
15798024
θ
8
+
8750427840
α
2
θ
6
-
164081808
θ
7
+
22019187960
α
2
θ
5
-
894701808
θ
6
+
25040008224
α
2
θ
4
-
2742600168
θ
5
+
3235380192
α
2
θ
3
-
4665991176
θ
4
-
13040818560
α
2
θ
2
-
3959356896
θ
3
-
5029516800
α
2
θ
-
1062129024
θ
2
+
359764416
θ
+
131725440
f
10
q
2
=
479001600
2
θ
7
+
49
θ
6
+
483
θ
5
+
2450
θ
4
+
6769
θ
3
+
9849
θ
2
+
6534
θ
+
1260
.In a similar way to that for the nodes x
1 and x
2, the implicit compact finite difference approximation of second-order derivatives at nodes x
M and x
M
-
1 is constructed.


## 4. Implicit Compact Finite Difference Scheme for Boundary Nodes in the Case of Neumann Boundary Conditions

For the construction of fourth-order implicit compact finite difference approximation for second-order derivatives at nodex
i, i
=
2,3
,
…
,
M, we consider the following model: (21)
α
1
T
x
x
x
i
-
1
+
T
x
x
x
i
+
α
1
T
x
x
x
i
+
1
=
1
h
2
a
1
T
x
i
-
1
+
a
2
T
x
i
+
a
3
T
x
i
+
1
.The stencil in (21) is(22)
i
-
1
,
i
,
i
+
1and its length is three. It means that we have to construct a fourth-order approximation of second-order derivative at nodes x
1 and x
M. Due to symmetry, we will only present the construction for the node x
1 in the case of Neumann boundary conditions. By expanding (21) around the node x
i and comparing the coefficients of the same-order derivative of T at x
i, we find the following values:(23)
α
1
=
1
10
,
a
1
=
6
5
,
a
2
=
-
12
5
,
a
3
=
6
5
, and (21) becomes(24)
1
10
T
i
-
1
′
′
+
T
i
′
′
+
1
10
T
i
+
1
′
′
=
6
5
T
i
-
1
-
2
T
i
+
T
i
+
1
h
2
1
10
1
1
10
T
i
-
1
′
′
T
i
′
′
T
i
+
1
′
′
=
6
5
h
2
1
-
2
1
T
i
-
1
T
i
T
i
+
1
,where i
=
2,3
,
…
,
M. Equation (24) can be written as(25)
1
10
1
1
10
0
⋯
0
0
1
10
1
1
10
⋯
0
⋱
⋱
⋱
0
⋯
0
1
10
1
1
10
T
1
′
′
T
2
′
′
⋮
T
M
′
′
=
6
5
h
2
1
-
2
1
0
⋯
0
0
1
-
2
1
⋯
0
⋱
⋱
⋱
0
⋯
0
1
-
2
1
T
1
T
2
⋮
T
M
.The error equation for (24) is given by(26)
E
3
=
1
200
f
i
6
h
4
+
O
h
5
.Next, we consider the following model for the construction of implicit compact finite difference approximation of the second-order derivative at node x
1:(27)
T
x
x
x
1
+
β
1
T
x
x
x
2
=
1
h
2
b
T
x
2
-
T
x
1
+
c
T
x
3
-
T
x
1
+
d
T
′
x
1
h
.After expanding (27) around x
1 and simplifying it, we get the values of the parameters:(28)
θ
=
-
1
+
1
15
ϕ
b
=
30
-
1
+
1
15
ϕ
3
+
90
-
1
+
1
15
ϕ
2
-
72
+
16
ϕ
5
β
=
15
2
-
1
+
1
/
15
ϕ
3
+
45
-
1
+
1
/
15
ϕ
2
2
-
10
+
2
5
ϕ
c
=
-
15
-
1
+
1
15
ϕ
3
-
135
-
1
+
1
/
15
ϕ
2
4
+
6
-
1
10
ϕ
d
=
-
45
-
1
+
1
/
15
ϕ
2
2
+
60
-
3
ϕ
,where ϕ
=
225
+
30
30. The error equation for the stencil (27) is given by(29)
E
4
=
30
ϕ
27000
+
1
40
-
ϕ
600
f
i
6
h
4
+
O
h
5
.


## 5. Stability

Consider a partitionP
=
x
0
,
x
1
,
x
2
,
…
,
x
M
-
1
,
x
M
,
x
M
+
1 of [
0
,
l
]. We define T
(
t
)
=
[
T
(
x
1
,
t
)
,
T
(
x
2
,
t
)
,
…
,
T
(
x
M
,
t
)
]
T and S
(
t
)
=
[
s
(
x
1
,
t
)
,
s
(
x
2
,
t
)
,
…
,
s
(
x
M
,
t
)
]
T. The implicit compact finite difference scheme for the second-order derivative can be written as (30)
A
T
x
x
t
=
-
1
h
2
B
T
t
T
x
x
t
=
-
1
h
2
A
-
1
B
T
t
.The vector form of 1D heat conduction equation (1a) can be written as(31)
T
t
t
=
k
T
x
x
t
+
S
t
T
t
t
=
-
k
h
2
A
-
1
B
T
t
+
S
t
.For time integration, we adopt Crank-Nicholson method for (31):(32)
T
m
+
1
-
T
m
Δ
t
=
-
k
h
2
A
-
1
B
T
m
+
1
-
T
m
2
+
S
m
+
1
+
S
m
2
I
+
r
A
-
1
B
T
m
+
1
=
I
-
r
A
-
1
B
T
m
+
S
m
+
1
/
2
T
m
+
1
=
I
+
r
A
-
1
B
-
1
I
-
r
A
-
1
B
T
m
+
S
m
+
1
/
2
,where U
m
=
U
(
t
0
+
m
Δ
t
), S
m
+
1
/
2
=
I
+
r
A
-
1
B
-
1
S
m
+
1
+
S
m
/
2, and r
=
k
Δ
t
/
2
h
2. For stability analysis, we assume two different solutions T
1 and T
2 of (32). Equation (32) can be written as(33)
T
1
m
+
1
=
I
+
r
A
-
1
B
-
1
I
-
r
A
-
1
B
T
1
m
+
S
m
+
1
/
2
(34)
T
2
m
+
1
=
I
+
r
A
-
1
B
-
1
I
-
r
A
-
1
B
T
2
m
+
S
m
+
1
/
2
.By subtracting (34) from (33), we get(35)
T
1
m
+
1
-
T
2
m
+
1
=
I
+
r
A
-
1
B
-
1
I
-
r
A
-
1
B
T
1
m
-
T
2
m
T
1
m
+
1
-
T
2
m
+
1
=
I
+
r
A
-
1
B
-
1
I
-
r
A
-
1
B
m
T
1
0
-
T
2
0
.If we assume that λ is the eigenvalue of A
-
1
B, then (
1
-
r
λ
)
/
(
1
+
r
λ
) is the eigenvalue of I
+
r
A
-
1
B
-
1
I
-
r
A
-
1
B. By definition of r as positive and for stability of the algorithm, we have to show that(36)
1
-
r
λ
1
+
r
λ
<
1
1
-
r
λ
1
2
+
λ
2
2
<
1
+
r
λ
1
2
+
λ
2
2
λ
1
>
0
,where λ
=
λ
1
+
i
λ
2. It means that we have to show that the real part of the eigenvalues of A
-
1
B is positive. According to Taussky theorem [18], an irreducible diagonally dominant matrix is positive definite if every main diagonal entry of the matrix is real and positive; then, every eigenvalue of the matrix has positive real part; that is, the matrix is positive definite. It can be seen that matrices (47), (48), (49), and (50) are positive definite (positive definite in the sense that the eigenvalues have positive real part) according to Taussky theorem. In fact, matrices (47), (48), (49), and (50) (i)
are irreducible,
(ii)
are diagonally dominant with at least one strict diagonal dominance,
(iii)
have real diagonal entries that are positive.Theorem 1.
Prove that if matricesA and B are positive definite, then matrix A
-
1
B is also positive definite in the sense that its eigenvalues have positive real parts.Proof.
Letλ
=
λ
1
+
i
λ
2 be the eigenvalue of A
-
1
B with corresponding eigenvector x
=
x
1
+
i
x
2. Then, we have(37)
A
-
1
B
x
=
λ
x
B
x
=
λ
A
x
(38)
B
x
1
=
λ
1
A
x
1
-
λ
2
A
x
2
B
x
2
=
λ
2
A
x
1
+
λ
1
A
x
2
x
1
T × (38) ⇒
(39)
x
1
T
B
x
1
=
λ
1
x
1
T
A
x
1
-
λ
2
x
1
T
A
x
2
(40)
x
1
T
B
x
2
=
λ
2
x
1
T
A
x
1
+
λ
1
x
1
T
A
x
2
x
2
T × (38) ⇒
(41)
x
2
T
B
x
1
=
λ
1
x
2
T
A
x
1
-
λ
2
x
2
T
A
x
2
(42)
x
2
T
B
x
2
=
λ
2
x
2
T
A
x
1
+
λ
1
x
2
T
A
x
2
λ
1 × (39) + λ
2 × (40) ⇒
(43)
λ
1
x
1
T
B
x
1
+
λ
2
x
1
T
B
x
2
=
λ
1
2
+
λ
2
2
x
1
T
A
x
1
λ
1 × (42) + λ
2 × (41) ⇒
(44)
λ
1
x
2
T
B
x
2
-
λ
2
x
2
T
B
x
1
=
λ
1
2
+
λ
2
2
x
2
T
A
x
2
.The quadratic forms x
1
T
A
x
1, x
2
T
A
x
2, x
1
T
B
x
1, and x
2
T
B
x
2 are positive, because matrixes A and B are positive definite. From (43) and (44), we have (45)
λ
1
+
λ
2
x
1
T
B
x
2
x
1
T
B
x
1
>
0
(46)
λ
1
-
λ
2
x
2
T
B
x
1
x
2
T
B
x
2
>
0
.Suppose that the eigenvalues of A
-
1
B are real; that is, λ
2
=
0; then, (45) and (46) imply that λ
1
>
0. But if the eigenvalues of A
-
1
B are not real, then we have the required positivity of λ
1 from the fact that the signs of λ
2, x
1
T
B
x
2, and x
2
T
B
x
1 are arbitrary for arbitrary x
1 and x
2. The arbitrariness of signs does not affect inequalities (45) and (46); hence, we conclude that λ
1
>
0 which completes the proof.


## 6. Numerical Simulations

To observe the validity and accuracy of our proposed finite difference schemes, we solve the following two test problems, one with Dirichlet boundary conditions and the other one with Neumann boundary conditions:(47)
A
=
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
2
5
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
β
3
α
3
1
α
3
β
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
2
5
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
(48)
B
=
4.094
-
2.303
0.556
-
0.160
0.034
-
0.003
-
0.0004
0.0001
0
0
0
0
-
0.141
2.008
-
1.371
0.088
0.030
-
0.019
0.005
-
0.0005
0
0
0
0
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
0
0
0
0
0
0
0
0
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
0
0
0
0
0
0
0
0
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
0
0
0
0
0
0
0
0
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
0
0
0
0
0
0
0
0
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
0
0
0
0
0
0
0
0
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
0
0
0
0
0
0
0
0
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
0
0
0
0
0
0
0
0
b
3
4
a
3
-
2
a
3
-
b
3
2
a
3
b
3
4
0
0
0
0
-
0.0005
0.005
-
0.019
0.030
0.088
-
1.371
2.008
-
0.141
0
0
0
0
0.0001
-
0.0004
-
0.003
0.034
-
0.160
0.556
-
2.303
4.094
(49)
A
=
1.0
0.366
0
0
0
0
0
0
0
0
0.1
1.0
0.1
0
0
0
0
0
0
0
0
0.1
1.0
0.1
0
0
0
0
0
0
0
0
0.1
1.0
0.1
0
0
0
0
0
0
0
0
0.1
1.0
0.1
0
0
0
0
0
0
0
0
0.1
1.0
0.1
0
0
0
0
0
0
0
0
0.1
1.0
0.1
0
0
0
0
0
0
0
0
0.1
1.0
0.1
0
0
0
0
0
0
0
0
0.1
1.0
0.1
0
0
0
0
0
0
0
0
0.366
1.0
(50)
B
=
1.233
-
1.0342
-
0.19872
0
0
0
0
0
0
0
-
1.2
2.4
-
1.2
0
0
0
0
0
0
0
0
-
1.2
2.4
-
1.2
0
0
0
0
0
0
0
0
-
1.2
2.4
-
1.2
0
0
0
0
0
0
0
0
-
1.2
2.4
-
1.2
0
0
0
0
0
0
0
0
-
1.2
2.4
-
1.2
0
0
0
0
0
0
0
0
-
1.2
2.4
-
1.2
0
0
0
0
0
0
0
0
-
1.2
2.4
-
1.2
0
0
0
0
0
0
0
0
-
1.2
2.4
-
1.2
0
0
0
0
0
0
0
-
0.19872
-
1.0342
1.233
(51)
Dirichlet  problem
=
T
t
x
,
t
=
T
x
x
x
,
t
,
0
<
x
<
1
,
t
>
0
,
T
x
,
0
=
sin
⁡
π
x
,
0
≤
x
≤
1
,
T
0
,
t
=
T
1
,
t
=
0
,
t
>
0
.The analytical solution to the above Dirichlet problem is T
(
x
,
t
)
=
exp
⁡
-
π
2
t
sin
⁡
(
π
x
). And(52)
Neumann  problem
=
T
t
x
,
t
=
T
x
x
x
,
t
,
0
<
x
<
1
,
t
>
0
,
T
x
,
0
=
cos
⁡
π
x
,
0
≤
x
≤
1
,
T
x
0
,
t
=
T
x
1
,
t
=
0
,
t
>
0
.The analytical solution to the above Neumann problem is T
(
x
,
t
)
=
exp
⁡
-
π
2
t
cos
⁡
(
π
x
).In Table1, we computed the spatial rate of convergence of our proposed implicit compact finite difference scheme for the Dirichlet problem. In all cases, Table 1 shows that the spatial rate of convergence is at least eight. In Table 2, we measure the temporal rate of convergence of Crank-Nicholson method which is at least two and it is according to the theoretical rate of convergence. Similarly, in Tables 3 and 4, we computed the spatial and temporal rates, respectively, for Neumann problem and found them according to theoretical values. In the case of Dirichlet problem, the variation of θ does not affect error in temporal dimension but does affect the error in spatial dimension and that is why Table 2 shows the same values of error with respect to different values of θ and Δ
t. We also observe that as we decrease either h, Δ
t, or both, we get reduction in the error E
(
M
,
Δ
). Decreasing of h means that we increase the number of grid points M in the spatial domain and this is also valid in temporal dimension. The temporal domain for the Dirichlet and Neumann problems is chosen to [
0,1
] as the solution decays rapidly when the time passes 1. The maximum error in the numerical solution occurs around 0.1 in temporal dimension; this is the reason why we just integrate the Dirichlet problem in the temporal domain [
0,0.1
] (see Table 1).Table 1
Dirichlet problem: Max.l
2-norm errors and convergence rate when t
=
0.1 and Δ
t
=
1.0
e
-
6.
M
E
(
M
,
Δ
)
Rate
E
(
M
,
Δ
)
Rate
θ
=
0.1
θ
=
0.2
10
1.6534
e
-
04
—
2.4968
e
-
04
—
20
4.3906
e
-
07
8.5568
7.8152
e
-
07
8.3196
40
8.6154
e
-
10
8.9933
8.0111
e
-
10
9.9301
θ
=
0.3
θ
=
0.4
10
3.2357
e
-
04
—
3.6276
e
-
04
—
20
1.0791
e
-
06
8.2281
1.2556
e
-
06
8.1746
40
9.2171
e
-
10
10.193
1.1166
e
-
09
10.135
θ
=
0.5
θ
=
0.6
10
3.5225
e
-
04
—
2.8316
e
-
04
—
20
1.2494
e
-
06
8.1393
1.0109
e
-
06
8.1299
40
1.2167
e
-
09
10.000
1.1107
e
-
09
9.8299
θ
=
0.7
θ
=
0.8
10
1.5954
e
-
04
—
1.9640
e
-
04
—
20
5.4419
e
-
07
8.1956
9.1503
e
-
07
7.7457
40
7.4804
e
-
10
9.5068
8.6751
e
-
10
10.043
θ
=
0.9
θ
=
0.95
10
4.5197
e
-
04
—
6.1645
e
-
04
—
20
2.0820
e
-
06
7.7621
2.8818
e
-
06
7.7409
40
2.4668
e
-
09
9.7211
3.6858
e
-
09
9.6108Table 2
Dirichlet problem: Max.l
2-norm errors and convergence rate when M
=
1000 and t
=
1.
Δ
t
E
(
M
,
Δ
)
Rate
E
(
M
,
Δ
)
Rate
θ
=
0.1
θ
=
0.2
1.0
e
-
2
4.2273
e
-
04
—
4.2273
e
-
04
—
5.0
e
-
3
1.0560
e
-
04
2.0011
1.0560
e
-
04
2.0011
2.5
e
-
3
2.6395
e
-
05
2.0003
2.6395
e
-
05
2.0003
θ
=
0.3
θ
=
0.4
1.0
e
-
2
4.2273
e
-
04
—
4.2273
e
-
04
—
5.0
e
-
3
1.0560
e
-
04
2.0011
1.0560
e
-
04
2.0011
2.5
e
-
3
2.6395
e
-
05
2.0003
2.6395
e
-
05
2.0003
θ
=
0.5
θ
=
0.6
1.0
e
-
2
4.2273
e
-
04
—
4.2273
e
-
04
—
5.0
e
-
3
1.0560
e
-
04
2.0011
1.0560
e
-
04
2.0011
2.5
e
-
3
2.6395
e
-
05
2.0003
2.6395
e
-
05
2.0003
θ
=
0.7
θ
=
0.8
1.0
e
-
2
4.2273
e
-
04
—
4.2273
e
-
04
—
5.0
e
-
3
1.0560
e
-
04
2.0011
1.0560
e
-
04
2.0011
2.5
e
-
3
2.6395
e
-
05
2.0003
2.6395
e
-
05
2.0003
θ
=
0.9
θ
=
0.95
1.0
e
-
2
4.2273
e
-
04
—
4.2273
e
-
04
—
5.0
e
-
3
1.0560
e
-
04
2.0011
1.0560
e
-
04
2.0011
2.5
e
-
3
2.6395
e
-
05
2.0003
2.6395
e
-
05
2.0003Table 3
Neumann problem: Max.l
2-norm errors and convergence rate when t
=
1, Δ
t
=
1.0
e
-
5, and θ
=
3.1541
e
-
01.
M
E
(
M
,
Δ
)
Rate
40
3.5993
e
-
08
—
80
2.2134
e
-
09
4.0234
160
5.5122
e
-
11
5.3275Table 4
Neumann problem: Max.l
2-norm errors and convergence rate when t
=
1, M
=
1000, and θ
=
3.1541
e
-
01.
Δ
t
E
(
M
,
Δ
)
Rate
1.0
e
-
2
2.1144
e
-
04
—
5.0
e
-
3
5.2819
e
-
05
2.0011
2.5
e
-
3
1.3202
e
-
05
2.0003


## 7. Conclusions

The implicit compact finite difference methods provide a more accurate way to approximate the spatial derivatives compared to explicit finite difference methods. The construction of compact finite difference operators for the interior nodes provides diagonal dominance and positivity of the diagonal entries which is in fact a very nice property which finally appears in the form of positive deftness of the compact operator. We have observed that the positive deftness helps us to prove the stability of numerical algorithm to solve 1D heat conduction equations. However, the diagonal dominance and positivity of diagonal entries for the interior nodes are not enough, because we also have to deal with the boundary conditions and usually one-sided compact finite difference schemes do not respect the nice property of positive definiteness with high order of convergence rate. In this project, the designing of the high-order accurate boundary conditions is established in such a way that we can maintain the positive definiteness of compact operator for the numerical scheme.



---

*Source: 8376061-2016-05-16.xml*