This work is licensed under the Creative Commons Attribution 4.0 International License.
Introduction
This manuscript is going to dwell on the following 1D heat equation
{{\partial u} \over {\partial t}} = {\alpha ^2}{{{\partial ^2}u} \over {\partial {x^2}}},\;\;\;x \in (0,L),\;\;\;t > 0
with initial condition
u(x,0) = f(x),\;x \in [0,L]
and boundary conditions
u(0,t) = u(L,t) = 0,\;\;\;t > 0
where α stands for the rod thermal diffusivity and the f (x) stands for a previously defined function. Many analytical and numerical studies in the literature [1,2] clearly state that this Initial and Boundary Value Problem (IBVP) can be considered as one of the most outstanding Partial Differential Equations (PDEs). The equation can mostly be encountered in various fields of applied sciences and mathematics for describing the temperature change in other words, the heat distribution and throughout a given solution domain for a given temporal domain.
This work is going to take into consideration the flow of the heat in 1D insulated thoroughly except with both ends of the investigated rod. The found solutions of this paper are going to be presented in terms of functions given along the rod in spatial and temporal direction. One can obviously see that this equation has got a vital importance in various areas of scientific world such as chemistry, physics and mechanics etc. The heat equation can be seen as a prototype for parabolic PDE in several fields of science. Due to its wide utilization and vital place, both the numerical and exact solutions of these types of equations have taken attention and essentially have become important in order to investigate various natural phenomena. Eq. (1) has at the same time been handled for economical problems in order to simulate several models with different choices [3]. The heat equation can also be described as the flow of heat given in a rod having diffusion α2uxx through the rod having the parameter α representing the thermal diffusivity in this rod and l stands for the length of the considered rod [4]. Kaskar [5] has stated that electrodynamics, fluid flow, elasticity, electrostatics and other natural phenomena rely mostly on these equations. This prototype problem has been considered both analytically and numerically for many years by several researchers. Despite this, the heat equation is still an outstanding one because of the fact that several natural phenomena can be stated as PDEs given together with appropriate conditions given at the initial time and the boundaries of the solution domain. Çağlar et al. [3] have introduced a novel numerical algorithm developed mainly with the collocation technique by utilizing cubic polynomial classical B-splines to produce approximate solutions simulating the physical treatment of the initial and boundary value problem for the heat conduction equation in one-dimension space. Dhawan and Kumar [2] have studied temperature change by virtue of the cubic B-spline finite elements method (FEM) and this gave satisfactory results. Suarez-Carreno and Rosales-Romero [6] have dealt with the heat equation as a role model for complex equations encountered in the field of finance as well as engineering. Unlike other researchers, they have presented a family of fully discretized explicit and implicit finite difference schemes producing accurate and stable approximate solutions, especially for sufficiently large time steps. Goh et al. [7] have introduced a fully dicretized stable numerical scheme utilizing the Crank-Nicolson type approach with cubic usual B-spline basis functions for spatial derivative for deriving approximate results of the heat and advection-diffusion equations governed in one-dimension. Lozada-Cruz et al. [8] have proved several results on the well-posedness of the 1D heat equation and associated stationary problem and then presented the application of the Crank-Nicolson for the numerical solution of the problem involving the Robin boundary condition. Hooshmandasl et al. [9] have derived a simple and fast scheme producing accurate results in terms of the mean error and maximum error norms constructed on the first kind Chebyshev wavelet collocation method by virtue of operational integration matrices for finding numerical approximate solutions of the classical heat equation given with Dirichlet boundary conditions. Han and Dai [10] have derived two compact schemes with fifth and sixth order accuracy based on Crank-Nicolson type approach for calculating much more accurate and unconditionally stable approximate solutions of two different initial-boundary value problems consisting of the heat equation with only homogeneous Dirichlet and only Neumann boundary conditions. Kutluay et al. [11] have developed an effective fully discretized scheme for producing approximate solutions of the heat equation by virtue of the usual collocation technique but utilizing the roots of both Legendre and Chebychev polynomials and third degree Hermite B-spline basis at the internal collocation mesh points. Sun and Zhang [12] have presented a hybrid scheme using the fourth order boundary value technique in the discretization of the temporal variable and the fourth order compact difference method in the discretization of the spatial variable to solve heat equations with two different constant coefficients. Patel and Pandya [13] have introduced a solution algorithm discretized with the cubic traditional spline polynomial function approach for the spatial derivative similar to the explicit and Crank-Nicolson type approximations and the usual forward difference approach for the time derivative to solve the heat equation given in one-dimension satisfying the prescribed initial and homogenous Dirichlet boundary conditions. Tarmizi et al. [14] have considered the initial and Dirichlet boundary value problem for the heat equation having an analytical solution, and proposed two computational schemes using spectral method and the Crank-Nicolson type approximation technique separately for the solutions of the problem. Yosaf et al. [15] have provided an implicit compact finite difference method since it gives a way which is more accurate in order to approximate the spatial derivative when compared with the explicit finite difference method.
This current article aims to develop a new robust computational scheme finding approximate solutions to the one-dimensional heat problem with initial and Dirichlet type boundary conditions given by Eqs.(1)–(3). To accomplish the goal, septic Hermite basis functions are going to be utilized. If one investigates its historical development, it can be seen that in 1946, first of all, Schoenberg initiated the basic theory of classical B-splines [2]. If one uses the iterative relation formulated in [16, 17] and is widely known by the names of the authors B-splines of any degree can easily be found. Recent years have witnessed a boom both in the theory of spline function and their applications for obtaining numerical solutions of the differential equations approximately such as those in [3]. Furthermore, various interesting problems have been studied using a wide class of approximate numerical solution techniques based on mainly FEM methods by different B-spline basis, especially with the odd-numbered degrees [18,19,20].
Many differential equations have been solved approximately by various methods, schemes and techniques using both the usual, trigonometric and Hermite spline basis functions. Although these methods have some drawbacks, they are worth utilizing in the application of numerical methods. One can also point out their advantages over other methods such as (i) one obtains a diagonal matrix and its storage and manipulation in digital computers are easy; (ii) the computer speed and digital storage area are in favor of programming.
Up to now, one-dimensional heat problem with initial and Dirichlet boundary conditions defined by Eqs.(1)–(3) and also some numerical techniques used in the literature to find the approximate solution have been mentioned. The works to be done in the remaining sections of this current article can be listed as follows: In the section of the description of the method, one dimensional heat problem governed by Eqs.(1)–(3) is fully discretized in both spatial and temporal variables. Such that the temporal variable is first discretized using Crank-Nicolson type approximation and then the spatial variable is discretized using the collocation finite element technique with septic Hermite polynomial base functions. In the application, the shifted roots of both Chebyshev and Legendre polynomials at the inner collocation mesh points are used. Additionally, the initial values required for starting the resulting recursive system of algebraic equations are computed using the initial condition of the governing equation. In the section of the Stability Analysis of the Scheme, the traditional von-Neumann technique is used and it is shown that the linear scheme is stable for every modest values of mesh sizes. In the section of the Illustrative Example and Its Numerical Results, the newly obtained numerical results using the current scheme are given in tables and compared with some ones available in the literature for the same values of parameters. The last section put forwards a brief conclusion about the present work and possible future studies.
Description of the method
Before proceeding with the description of the method to be presented, as in the abstraction of other existing numerical schemes, the solution region [0, l] × [0, T] of the problem at mesh points (xj, tn) for sufficiently small h = Δx and k = Δt must be discretized into equal width finite elements as
[0,l] = \{ {x_j}:{x_j} = jh,j = 1(1)M + 1,Mh = l\}
and
[0,T] = \{ {t_n}:{t_n} = nk,n = 0(1)N,Nk = T\} .
In this section, a useful and effective pointwise computational scheme will be given resulting in a fully discretized linear algebraic equation system for spatial and temporal integrations to be used in calculating the approximate solutions of the one-dimensional initial and Dirichlet boundary conditioned heat problem governed by Eqs. (1)–(3).
Spatial discretization
Let uM (x, t) denote an approximate solution to the analytical solution u(x, t). The approximate solution to be derived will be obtained by a combination of Crank-Nicolson type discretization for the time variable and septic Hermite B-spline collocation type finite element discretization utilizing the roots of the sixth degree of shifted Legendre polynomial at internal collocation points for the spatial variable.
To achieve the above mentioned goals, the numerical solution uM (x, t) will be searched in the format
{u_M}(x,t) = \mathop {\mathop {\sum }\limits_{j = 1} }\limits^{M + 4} {a_{j + 6k - 6}}\left( t \right){H_j}({x_i}) \approx u(x,t)
in which a’s stand for time dependent coefficients to be determined, k stands for the element number, i (i = 1(1)6) stands for the inner-collocation points and finally Hj (x) (j = 1(1)M + 1) are septic Hermite B-spline basis polynomial functions defined over the interval [xj−1, xj+1] as follows [21].
\matrix{ \hfill {{H_{6j + 1}}(x)\;\; = 35{{{{(x - {x_{j - 1}})}^4}} \over {{h^4}}} - 84{{{{(x - {x_{j - 1}})}^5}} \over {{h^5}}} + 70{{{{(x - {x_{j - 1}})}^6}} \over {{h^6}}} - 20{{{{(x - {x_{j - 1}})}^7}} \over {{h^7}}},\quad \quad {x_{j - 1}} \le x \le {x_j}} \cr \hfill {{H_{6j - 5}}(x)\;\;\; = 35{{{{({x_{j + 1}} - x)}^4}} \over {{h^4}}} - 84{{{{({x_{j + 1}} - x)}^5}} \over {{h^5}}} + 70{{{{({x_{j + 1}} - x)}^6}} \over {{h^6}}} - 20{{{{({x_{j + 1}} - x)}^7}} \over {{h^7}}},\quad \quad {x_j} \le x \le {x_{j + 1}}} \cr \hfill {{H_{6j - 4}}(x) = - 15{{{{(x - {x_{j - 1}})}^4}} \over {{h^3}}} + 39{{{{(x - {x_{j - 1}})}^5}} \over {{h^4}}} - 34{{{{(x - {x_{j - 1}})}^6}} \over {{h^5}}} + 10{{{{(x - {x_{j - 1}})}^7}} \over {{h^6}}},\quad \quad {x_{j - 1}} \le x \le {x_j}} \cr \hfill {{H_{6j + 2}}(x) = - 15{{{{({x_{j + 1}} - x)}^4}} \over {{h^3}}} + 39{{{{({x_{j + 1}} - x)}^5}} \over {{h^4}}} - 34{{{{({x_{j + 1}} - x)}^6}} \over {{h^5}}} + 10{{{{({x_{j + 1}} - x)}^7}} \over {{h^6}}},\quad \quad {x_j} \le x \le {x_{j + 1}}} \cr \hfill {{H_{6j - 3}}(x)\;\;\;\;\;\;\; = {5 \over 2}{{{{(x - {x_{j - 1}})}^4}} \over {{h^2}}} - 7{{{{(x - {x_{j - 1}})}^5}} \over {{h^3}}} + {{13} \over 2}{{{{(x - {x_{j - 1}})}^6}} \over {{h^4}}} - 2{{{{(x - {x_{j - 1}})}^7}} \over {{h^5}}},\quad \quad {x_{j - 1}} \le x \le {x_j}} \cr \hfill {{H_{6j}}(x)\;\;\;\;\;\;\; = {5 \over 2}{{{{({x_{j + 1}} - x)}^4}} \over {{h^2}}} - 7{{{{({x_{j + 1}} - x)}^5}} \over {{h^3}}} + {{13} \over 2}{{{{({x_{j + 1}} - x)}^6}} \over {{h^4}}} - 2{{{{({x_{j + 1}} - x)}^7}} \over {{h^5}}},\quad \quad {x_j} \le x \le {x_{j + 1}}} \cr \hfill {{H_{6j - 2}}(x)\;\;\;\;\; = - {1 \over 6}{{{{(x - {x_{j - 1}})}^4}} \over h} + {1 \over 2}{{{{(x - {x_{j - 1}})}^5}} \over {{h^2}}} - {1 \over 2}{{{{(x - {x_{j - 1}})}^6}} \over {{h^3}}} + {1 \over 6}{{{{(x - {x_{j - 1}})}^7}} \over {{h^4}}},\quad \quad {x_{j - 1}} \le x \le {x_j}} \cr \hfill {{H_{6j - 1}}(x)\;\;\;\; = - {1 \over 6}{{{{({x_{j + 1}} - x)}^4}} \over h} + {1 \over 2}{{{{({x_{j + 1}} - x)}^5}} \over {{h^2}}} - {1 \over 2}{{{{({x_{j + 1}} - x)}^6}} \over {{h^3}}} + {1 \over 6}{{{{({x_{j + 1}} - x)}^7}} \over {{h^4}}},\quad \quad {x_{j - 1}} \le x \le {x_{j + 1}}.} \cr }
Throughout this paper, the 1D heat conduction equation which has been widely given by Eq. (1) by the appropriate initial (2) and boundary conditions (3) is going to be handled. For this purpose, the finite element collocation method is selected as a tool. While applying this method, the interval [a, b] of the problem is discretized in M equal width finite elements by means of the mesh points xj (j = 1(1)M + 1) in such a way that xl = x1 < x2 < ⋯ < xM < xM+1 = xr and h = xj+1 − xj. To be more precise, a non-uniform partition of the solution domain might also be preferred. However, due to the fact that the non-uniform selection would have increased the digital computer memory storage requirement and at the same time simulation duration, the uniform partition is commonly preferred. In place of analytical solution u(x, t), a numerical approximation uM (x, t) will be used by means of the septic Hermite basis functions as
{u_M}(x,t) = \mathop {\mathop {\sum }\limits_{j = 1} }\limits^{M + 4} {a_{j + 6k - 6}}\left( t \right){H_j}({x_i}) \approx u(x,t).
The following roots
\xi _i^L
and
\xi _i^C\left( {i = 1\left( 1 \right)6} \right)
of the sixth-degree shifted Legendre and Chebyshev polinomials [22]
\xi _1^L = 0.033765242898424,\xi _2^L = 0.169395306766868,\xi _3^L = 0.380690406958402\xi _4^L = 0.619309593041598,\xi _5^L = 0.830604693233132,\xi _6^L = 0.966234757101576
and
\xi _1^C = {1 \over 2} - \sqrt {{{\sqrt 3 } \over {16}} + {1 \over 8}} = 0.017037086855466,\xi _2^C = {1 \over 2}(1 - {{\sqrt 2 } \over 2}) = 0.146446609406726\xi _3^C = {1 \over 2} - \sqrt {{1 \over 8} - {{\sqrt 3 } \over {16}}} = 0.37059047744874,\xi _4^C = {1 \over 2} + \sqrt {{1 \over 8} - {{\sqrt 3 } \over {16}}} = 0.62940952255126\xi _5^C = {1 \over 2}(1 + {{\sqrt 2 } \over 2}) = 0.853553390593274,\xi _6^C = {1 \over 2} + \sqrt {{{\sqrt 3 } \over {16}} + {1 \over 8}} = 0.982962913144534
are used, respectively.
In the construction process for the algorithm of the present method, firstly the Crank-Nicolson type finite difference formula for the semi-discretizing of the governed equation and then for fully discretizing the equation, the collocation technique with septic Hermite spline basis functions on spatial direction will be used. This choice has the advantages of several vital characteristics such as easy to handle algorithms which produce and low level of storage requirement. Besides those advantages, both of the linear and non-linear systems resulted from the usage of splines are usually not ill-conditioned and thus permit the required coefficients to be found out in an easy manner. Furthermore, the newly obtained approximate solutions generally don’t tend in numerical instability. Last but not least, the matrix systems resulting from the usage of splines are usually sparsly loaded band matrixes and allow easy implementation in digital computers [18].
Temporal discretization
When the aforementioned 1D heat conduction equation given by Eq. (1) as follows is discretized
{u_t} - {\alpha ^2}{u_{xx}} = 0.
one can use Crank-Nicolson type formula. Firstly, one discretizes Eq. (1) as
{{{u^{n + 1}} - {u^n}} \over k} - {\alpha ^2}\left[ {{{{{({u_{xx}})}^{n + 1}} + {{({u_{xx}})}^n}} \over 2}} \right] = 0.
If one separates the variables such that the next time level ones are placed on the left hand side and the ones for previous time level are placed on the right hand side, then
{{{u^{n + 1}}} \over k} - {\alpha ^2}{{{{({u_{xx}})}^{n + 1}}} \over 2} = {{{u^n}} \over k} + {\alpha ^2}{{{{({u_{xx}})}^n}} \over 2}
is obtained.
The newly obtained system (8) is obviously recursive in nature. Thus, the element coefficients vector
{{\bf{a}}^n} = (a_1^n, \cdots ,a_{6M + 1}^n,a_{6M + 2}^n)
can be found in a recursive manner, where tn = nk, n = 0(1)N till the desired time T. If one utilizes the conditions presented at the boundary of the solution domain given by Eq.
(3) and eliminates only the coefficients
a_1^n
and
a_{6M + 1}^n
in Eq. (8) as follows: Using the boundary condition given at the left of the solution domain gives
\matrix{ {u({x_l},t) = a_1^n{H_1}({x_l}) + a_2^n{H_2}({x_l}) + a_3^n{H_3}({x_l}) + a_4^n{H_4}({x_l})} \hfill \cr {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + a_5^n{H_5}({x_l}) + a_6^n{H_6}({x_l}) + a_7^n{H_7}({x_l}) + a_8^n{H_8}({x_l}) = 0.} \hfill \cr }
Since Hm1 = 0 for m = 2(1)8 and H11 ≠ 0, the condition
a_1^n = 0
is obtained. In a similar way, using the right boundary condition results in
\matrix{ {u({x_r},t) = a_{6M - 5}^n{H_1}({x_r}) + a_{6M - 4}^n{H_2}({x_r}) + a_{6M - 3}^n{H_3}({x_r}) + a_{6M - 2}^n{H_4}({x_r})} \hfill \cr {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + a_{6M - 1}^n{H_5}({x_r}) + a_{6M}^n{H_6}({x_r}) + a_{6M + 1}^n{H_7}({x_r}) + a_{6M + 2}^n{H_8}({x_r}) = 0.} \hfill \cr }
Since Hm6 = 0 for m = 1(1)8 but m ≠ 7 and H76 ≠ 0, the condition
a_{6M + 1}^n = 0
is found. Finally, the following solvable square algebraic system is obtained of the form
{\bf{L}}{{\bf{a}}^{n + 1}} = {\bf{R}}{{\bf{a}}^n}.
In this iterative system, L and R are square 6M × 6M diagonal band matrices and, an+1 and an are 6M column matrices.
Now, the unknown coefficients a can be obtained by solving Eq.(9). To be able to start this iterative calculation, the initial vector a0 should be found. This initial vector is constructed with the usage of the initial condition given together with the governing equation. The iterative procedure is made in a successive manner for tn = nk (n = 0(1)N) until the desired time T. To obtain the solution of this system, an algorithm written in MatlabR2015a is run on a digital computer.
Therefore, the required initial coefficients to initiate the scheme given by Eq. (9) are calculated by Eq. (10) and next the iterative procedure to generate numerical solutions is made in a repetitive manner till the final time T.
Illustrative application and its numerical results
The present section will present the application of the above proposed numerical scheme based on the collocation FEM using septic Hermite basis functions to produce some numerical results of the test problem governed by Eqs. (1)–(3). For the test problem, all computations have been done for f (x) = sin(πx), (x, t) = [0, 1] × [0, T] and α = 1.
Due to the fact that the test problem has an analytical solution of the form [10, 15]
u(x,t) = {{\sin (\pi x)} \over {{e^{{\alpha ^2}{\pi ^2}t}}}}
the error norms L2 and L∞ formulated as, respectively, will be used to test the accuracy and validity of this method
{L_2} = \sqrt {h\sum\limits_{j = 1}^M {{\left| {{u_j} - {{({u_M})}_j}} \right|}^2}} ,\quad \quad {L_\infty } = \mathop {\max }\limits_{1 \le j \le M} \left| {{u_j} - {{({u_M})}_j}} \right|.
Moreover the rate of convergence is calculated by the following formula [23] and found to be near to its expected value
Rate\;of\;Convergence \simeq {{\log (E_j^ * /E_{j + 1}^ * )} \over {\log ({h_j}/{h_{j + 1}})}}
where j corresponds to the spatial partition and
E_j^ *
denotes the error norm.
This manuscript carries out all of the numerical computations using both Septic Hermite Collocation Method with Legendre roots (SHCM-L) and Septic Hermite Collocation Method with Chebyshev roots (SHCM-C). Throughout this work, all calculations have been made with MATLAB R2015a on Intel (R) Core(TM) i5 3330SU CPU @2.70Ghz computer having 6 GB of RAM. To control the effectiveness and accuracy of the present method, the numerical results obtained by the scheme (9) using various space step sizes, time step sizes and desired values of final time are displayed in tables and compared with some of those existing in the literature.
Table 1 illustrates a clear comparison of the newly computed error norm values L2 of this scheme with those given in [11,15] for various values of k = 1/100, 5/1000, 25/10000 and h = 1/1000 at T = 1. It is obvious from this table that when the values of Δt decrease, at the same time the error norm L2 values decrease. Furthermore, it is seen that the obtained values using both SHCM-C and SHCM-L are clearly better than those in [15].
A comparison of the computed error norm L2 by the proposed scheme with those given in [11, 15] for h = 1/1000 and k = 0.01, 0.005, 0.0025 at T = 1.
Table 2 displays a comparison of the error norm values L2 calculated by the current scheme with those in [10, 11] for several values of h = 1/5, 1/10, 1/20 and k = 1/106 at T = 1. It can be easily observed from the table that when the number of elements is increased, the values of the error norm L2 decrease. Similarly, this table clearly shows that the results found out by utilizing the present scheme are in general better than those of compared ones among others as well as in [10].
A comparison of the computed error norm L2 of the proposed scheme with those given in [10, 11] for M = 5, 10, 20 and k = 1/106 at T = 1.
A comparison of the values of error norm L2 calculated from numerical solutions obtained by the current scheme with those given in [11, 15] for different values h = 1/10, 1/20, 1/40 and k = 1/106 at T = 0.1 is listed in Table 3. This table shows that the obtained values of norm L2 are sufficiently small and also much better than ones given in [15].
A comparison of the computed error norm L2 of the proposed scheme with those given in [11, 15] for M = 10, 20, 40 and k = 1/106 at T = 0.1.
Table 4 compares the calculated values of error norm L∞ with those given in [7, 11, 12] for various values of h = k at T = 1. It can be obviously understood from this table that the computed values of L∞ norm are in very good harmony with those in [11] and generally better than those in [7, 12].
A comparison of the computed error norm L∞ of the proposed scheme with those given in [7, 11, 12] for various values of h = k at T = 1.
Finally, Table 5 presents the calculated values of L2 and L∞ with those given in [9, 11] for M = 16 and k = 0.01 at various T. When the table is investigated, it is obvious that both of the newly computed norms are inherently small enough. They match with those in [11] and much better than those in [9]. It should be emphasized from the tables that the numerical results also approach to their analytical ones as the step size in spatial dimension is decreased.
A comparison of the computed error norms L2 and L∞ of the proposed scheme with those given in [9, 11] for h = 1/16, k = 0.01 at various T.
This work has successfully applied the collocation finite element method with septic Hermite spline base functions resulting in an implicit linear algebraic system to find out the numerical solutions for the 1D heat conduction equation. The obtained approximate numerical solutions show that the current scheme derives very good results which are in very good harmony with the analytical solutions and better than those given in comparison tables. As a conclusion, the newly presented numerical scheme, with easy implementation, results in efficient and accurate solutions. For a prospective future study, the current method can be successfully implemented to obtain accurate numerical solutions of various nonlinear PDEs as well as linear ones that play a vital role in describing natural phenomena which are widely encountered in our daily life.
