Efficient Finite Difference Methods for the Numerical Analysis of One-Dimensional Heat Equation ()
1. Introduction
The heat equation is a parabolic second-order linear partial differential equation that describes the distribution of heat (or temperature variation) in a given region over time. It demonstrates the fact of how heat flows through an object. For example, if a heated body is placed in a box of cold water, the body’s temperature will decrease. Eventually, the temperature in the box will equalize at a time if no other heat source is applied there.
To solve the one-dimensional heat equation, many effective methods can be used, such as separation of variables [1] [2] [3] [4] , Fourier series methods [5] [6] , domain decomposition method [7] [8] [9] , finite difference method [4] [10] [11] [12] [13] , finite element method [14] [15] , finite volume method [16] [17] and many more. In this paper, we use the finite difference method to solve and analyze the one-dimensional time-dependent heat equation. The finite difference method (FDM) is a powerful numerical technique for obtaining the approximate solution to the initial boundary value problems. It is easy to solve PDEs using this method [18] . Taylor series approach is used to find the solutions of a partial differential equation in this method [19] . The technique of this method is that all the derivatives in a partial differential equation are replaced by their finite difference approximations, and the PDEs are converted to a set of linear algebraic equations. Any iterative procedure can be used to solve this system of linear equations. Then, the solution of this system of linear equations is referred to as the solution of the PDE [1] . The finite difference method has been extensively used by researchers over the past years due to the fact that most of the calculations of this method can be carried out on the computer, and hence the solution is easy to obtain [20] [21] .
The aim of this work is to investigate a one-dimensional heat equation by using explicit Forward Time Centered Space (FTCS) and implicit Crank-Nicolson (CN) finite difference schemes. The rest of this paper is organized as follows. We present the one-dimensional heat equation with Neumann boundary conditions in the next Section. In Section 3, a brief introduction of thermal diffusivity is reported. The derivatives of the temperature distribution
are approximated by finite differences using the Taylor series in Section 4. The explicit and Crank-Nicolson schemes are developed for the one-dimensional heat equation in Section 5. The stability criterion is discussed, and the stability conditions for both schemes are verified in Section 6. We have exhibited the results and finally compared the results between the exact and the approximate solutions. The algorithm is executed by standard computer programming languages in Sections 6 and 7. We estimate error between the exact and the approximate solutions for a specific numerical problem and show it graphically in Section 7. Finally, a short conclusion is reported in Section 8.
2. Model Equations
In our model problem, we consider the one-dimensional time-dependent heat conduction equation of the form [3] [11] [12] [21] ,
(1)
where,
and
, with initial condition
(2)
and boundary conditions
(3)
(4)
where,
is the temperature of the substance, α2 is called the thermal diffusivity of the substance. Equations (1)-(4) together are called initial boundary-value problem.
3. Thermal Diffusivity (α2)
Thermal diffusivity is a significant physical quantity that describes the heat transfer rate through a medium. Different materials have different thermal diffusivity. It tells us how fast heat transfer occurs within a material from the hotter to the colder side. Thus, thermal diffusivity measures how quickly a material reacts to temperature changes.
It is defined as the ratio between the thermal conductivity and the product of the density and the specific heat capacity of the body with the pressure held constant, mathematically,
(5)
where,
K0 = thermal conductivity of the substance;
ρ = mass per unit volume of the substance;
cp = substance’s specific heat capacity.
4. Finite Difference Discretization
In order to describe the finite difference techniques, let us consider the rectangular domain
into a finite number of nodes
, where Ω is spatial domain and Γ is temporal domain.
Now, we develop the approximations of finite difference formulas according to [22] .
Let I be an open interval and
be differentiable
times. Now, we expand u in Taylor series about a point x, evaluate at
and at
, respectively as follows,
(6)
(7)
Rearranging (6), we get
(8)
which is forward difference formula for
.
Similarly, from (7), we get the backward difference formula as
(9)
By adding (6) and (7), we have
This gives us the finite difference formula for the second derivative with second order accuracy,
(10)
Now, we can define
as,
Thus, we get,
(11)
which is called a forward difference approximation for
.
In order to find an approximation to the second derivative,
; we start with the forward difference
(12)
We need to approximate the terms in the numerator. It is customary to use a backward difference approximation. This is given by letting
in the forward difference form,
(13)
applying (13) to
and evaluated at
and
we have,
(14)
and
(15)
5. Schemes for the Heat Equation
In order to develop the schemes for solving the one-dimensional heat Equation (1), both the time and space derivatives are required to replace by their finite differences. To perform this task, we need to specify both the time and spatial locations of the u values in the finite difference formulas [1] . Therefore, we need to introduce subscript i and j to designate the spatial and the time steps of the discrete solution, respectively. Hence, we can rewrite Equation (11) under these notations as
(16)
and Equation (10) as
(17)
For Crank-Nicolson scheme, the derivation of Equation (10) as
(18)
The second order partial derivative
can also be expressed by using the central difference at
and
and for any constant
(
),
(19)
(20)
where, the notation
is used to represent the numerical solution on a finite rod (substance)
and
.
5.1. Forward Time Centered Space (FTCS) Explicit Scheme
In this section, we consider the Forward in Time and Central in Space (FTCS) scheme, where we replace the time derivative in (1) by the forward difference scheme (16) and the space derivative in (1) by the central difference scheme (17). This yields,
(21)
Now we use the Courant-Friedrichs-Lewy or CFL condition in Equation (21).
The CFL condition for one dimension heat equation is
. By using this condition in Equation (21), we get
(22)
Equation (22) is called the Forward Time, Centered Space or FTCS approximation to the heat equation. The FTCS solution can be extracted from the grid of Figure 1.
Figure 1. Forward time centered space (FTCS) explicit scheme.
We exhibit the nature of the FTCS scheme (22) in Figure 2 and Figure 3 for the one-dimensional heat Equation (1) by the standard computer programming languages with boundary conditions
, and initial condition
,
.
5.2. Crank-Nicolson Scheme
To develop the Crank-Nicolson scheme [23] for our model Equation (1), we put
in Equation (20), for each
, [24] , we have,
(23)
Now, using (23) and (16) in Equation (1), we get
(24)
For the simplicity of calculation, we use the CFL condition,
, thus Equation (24) becomes,
(25)
Equation (25) is the Crank-Nicolson scheme for the one dimensional heat equation.
Setting
, i.e.,
then Equation (25) reduces to
(26)
(a)(b)
Figure 2. Temperature distribution and mesh plot for the explicit scheme in (a) and (b) are presented, respectively, for
and spatial domain
to
. (a) Temperature distribution; (b) Temperature distribution in mesh.
(a)(b)
Figure 3. Temperature distribution and mesh plot for the explicit scheme in (a) and (b) are presented, respectively, for
and spatial domain
to
. (a) Temperature distribution; (b) Temperature distribution in mesh.
The CN solution (26) can be made sense from the grid of Figure 4. Moreover, We show the temperature distribution of the one-dimensional heat Equation (1) by the standard computer programming languages using the CN scheme (26) under the boundary conditions
, and the initial conditions
,
in Figure 5 and Figure 6, respectively.
6. Stability Analysis
Stability is one of the most important properties of a system. It is a measure of the performance of a system which refers to the ability of the system to perform satisfactorily. It tells us the fact of how good the performance of the system is, whether it is stable or unstable. In this section, we discuss the Von-Neumann stability analysis for the one-dimensional heat equation for the explicit and Crank-Nicolson Schemes. The stability analysis was based on the citation [20] [24] [25] [26] .
From the definition of the error in the numerical approximation, we get
(27)
where,
is the exact solution at grid point
and
is the approximate solution.
For linear differential equations with periodic boundary condition, the spatial variation of error may be expanded in a finite Fourier series, in the interval Ω as
(28)
where the wave number
and
,
and x is the independent space variable,
is the complex exponential.
The time dependence of the error is included by assuming the amplitude of error
is a function of time. Since, the error tends to grow or decay exponentially with time, it is reasonable to assume that the amplitude varies exponentially with time; hence,
Figure 4. Crank-Nicolson implicit scheme.
(a)(b)
Figure 5. Temperature distribution and mesh plot for the Crank-Nicolson scheme in (a) and (b) are presented, respectively, for
and spatial domain
to
. (a) Temperature distribution; (b) Temperature distribution in mesh.
(a)(b)
Figure 6. Temperature distribution and mesh plot for the Crank-Nicolson scheme in (a) and (b) are presented, respectively, for
and spatial domain
to
. (a) Temperature distribution; (b) Temperature distribution in mesh.
(29)
Since the difference equation for error is linear (the behavior of each term of the series is the same as the series itself), it is enough to consider the growth of error of a typical term:
(30)
The stability characteristics can be studied using just this form for the error with no loss in generality.
We now define the amplification factor as,
(31)
Then, the necessary and sufficient condition for the error to remain bounded, maintaining numerical stability is,
(32)
6.1. Stability Criterion for FTCS Scheme
Now, we are going to show that Equation (32) holds for our FTCS numerical scheme (22). From Equation (30), we get by [20] [26] ,
(33)
(34)
(35)
(36)
Since, both the numerical solution
and the exact solution
satisfy Equation (22), therefore, the error
also follows the discretized PDE Equation (22). Thus we have,
(37)
Using Equations (33)-(36) in Equation (37), we get
Divide by
to yield after simplification,
(38)
We know that,
(39)
and
(40)
Subtracting (39) from (40) we get,
(41)
Putting
in (41), we get
squaring both sides, we get
(42)
From (38) and (42), we get
(43)
However, by Equations (33)-(36) in Equation (31), the last yield
(44)
Thus, by Equations (43), (44) plugged into (32), the condition for stability is given by
(45)
And from (45), we can write,
(46)
And
(47)
For the above conditions, Equations (46) and (47), to hold for all m (and therefore all
), we need to consider (46) and (47) separately.
Equation (46) yields
(48)
Note that the term
is always positive, i.e., always in the range [0, 1]. The worst case is when
, therefore, Equation (48) yields
(49)
Again Equation (47) yields
(50)
Since
range is [0, 1] Equation (50) yields
(51)
which always holds. Combining (49) and (51), we get
. Hence, the CFL condition
plays an important role to the stability analysis [27] .
Figure 7 and Figure 8 show that the explicit scheme is stable for
and
respectively. But, the explicit scheme is unstable when the value of r exceeds 0.5. The instability of the explicit scheme is shown in Figure 9 for
. The spatial domain of the substance is considered as
to
in Figures 7(a)-9(a) and the spatial domain
to
is used in Figures 7(b)-9(b), respectively.
We observe that our explicit numerical scheme satisfy the stability property. Hence, we can conclude that our scheme is numerically efficient.
6.2. Stability Criterion for Crank-Nicolson Scheme
By the Crank-Nicolson scheme, the discretized form of Equation (1)
Using similar arguments as in the FTCS method, the error
also follows the discretized ODE equation of the Crank-Nicolson scheme
(52)
Divide (52) by
,
(53)
Simplifying each term in Equation (53), we get
(54)
(55)
(56)
(57)
(58)
(a)(b)
Figure 7. Stability for the explicit numerical scheme for the spatial domain
to
in (a) and spatial domain
to
in (b), respectively, using
.
(a)(b)
Figure 8. Stability for the explicit numerical scheme for the spatial domain
to
in (a) and spatial domain
to
in (b), respectively, using
.
(a)(b)
Figure 9. Instability for the explicit numerical scheme for the spatial domain
to
in (a) and spatial domain
to
in (b), respectively, using
.
Plugging Equations (54)-(58) into Equation (53) yields
(59)
Using a trigonometric identity involving complex exponentials, such as
(60)
Using (60) into (59) yields
Simplifying, we get
(61)
Remembering that the cosine range is [−1, 1], the value we should explore for cosine values are
, hence we arrive at three facts.
Fact-01: If
, then from (61), we get
(62)
Since
, then
, always.
Fact-02: If
, then from (61), we get
(63)
Since
, then
, always.
Fact-03: If
, then from (61), we get
(64)
Then
, always. Therefore, we can simply say, the Crank-Nicolson scheme is unconditionally stable.
Figures 10-12 show that the Crank-Nicolson scheme is always stable. The mesh ratio
,
and
are used in Figures 10-12, respectively. The spatial domain is considered as
to
in Figures 10(a)-12(a) and the spatial domain
to
is used in Figures 10(b)-12(b), respectively.
The above figures illustrate that the Crank-Nicolson numerical scheme satisfies the stability conditions and thus we can conclude that the scheme is numerically efficient.
7. Error Analysis and Results Discussion
In this section, we focus on an error analysis between the exact solution and the approximate solution. To derive the relative error between the exact solution and the approximate solution for FTCS scheme, we compute the relative error in L1-norm defined by
(a)(b)
Figure 10. Stability for Crank-Nicolson scheme for the spatial domain
to
in (a) and spatial domain
to
in (b), respectively, using
.
(a)(b)
Figure 11. Stability for Crank-Nicolson scheme for the spatial domain
to
in (a) and spatial domain
to
in (b), respectively, using
.
(a)(b)
Figure 12. Stability for Crank-Nicolson scheme for the spatial domain
to
in (a) and spatial domain
to
in (b), respectively, using
.
where,
and
are the exact solution and the approximate solution computed by finite difference scheme, respectively.
The exact solution of heat Equation (1) with boundary conditions
, and initial condition
is [1] ,
(a)(b)
Figure 13. Comparison between the exact solution and the explicit numerical solution (a) and L1-norm of the error (b). (a) Comparison of the exact and the approximate solutions; (b) Error.
(65)
In Figure 13, we compare the exact and the approximate solution of the one-dimensional heat equation. In this case, we use FTCS explicit scheme to approximate the numerical solution. The L1-norm of the error is plotted in Figure 13(b). From the figure, it is clear that our numerical scheme converges and efficiently accurate.
8. Conclusion
In this work, two finite difference schemes for the one-dimensional heat equation are presented. In the first step, the explicit and Crank-Nicolson schemes have been developed. Later, the stability criterion is set, and the stability conditions are verified for both the schemes, which shows the validity of the numerical schemes. By using standard computer programming codes, we examine the stability criteria. We find that the explicit scheme is conditionally stable with
, andthe Crank-Nicolson scheme is unconditionally stable. Finally, we
discuss the error between the exact and the approximate solution of the one-dimensional heat equation, and validate the efficiency of our numerical schemes.