Abstract
We propose a numerical solution to the heat equation in polar cylindrical coordinates by using the meshless method of lines approach. The space variables are discretized by multiquadric radial basis function, and time integration is performed by using the Runge-Kutta method of order 4. In radial basis functions (RBFs), much of the research are devoted to the partial differential equations in rectangular coordinates. This work is an attempt to explore the versatility of RBFs in nonrectangular coordinates as well. The results show that application of RBFs is equally good in polar cylindrical coordinates. Comparison with other cited works confirms that the present approach is accurate as well as easy to implement to problems in higher dimensions.
1. Introduction
The importance of the study of the heat equation in polar cylindrical form is evident from the fact that it is a mathematical model of many real world problems. It is a model of the heat transfer which occurs in many areas ranging from heat transfer in a disc to bioheat transfer in biological systems. Several authors have attempted to solve the heat equation. Albasiny [1] used the Crank–Nicolson method to solve the cylindrical heat equation. Mitchell and Pearce [2] applied the explicit finite difference method to obtain numerical solution of the cylindrical heat conduction equation. Iyengar and Mittal [3] derived high accuracy implicit methods for the cylindrical heat equation. They developed unconditionally stable implicit formulas to solve problem. A symmetric semiimplicit finite difference method was used by Livne and Glasner [4] to solve the heat equation. Iyengar and Manohar [5] developed high-order difference methods for the solution of heat equation in polar cylindrical form. Manohar et al. [6] derived two-level, three-point difference schemes to solve the heat equation with linear variable coefficient. Al-Odat and Al-Nimr [7] studied thermal stability of composite superconducting tape by considering the heat conduction equation. Li et al. [8] obtained the numerical solution of Schrödinger and heat equation by using the finite difference scheme. They proposed high-order artificial boundary conditions on unbounded domains. In most of the cases, a mesh-based numerical method was used to investigate the heat equation. In this study, we propose a meshless method of lines (MMOL) approach to study the heat equation in polar cylindrical coordinates. Meshless numerical methods have advantage over the traditional mesh-based methods. These methods do not require mesh generation, and thus, they are more flexible and save computational time.
The method of lines (MOL) is a strong numerical technique that solves a partial differential equation (PDE) in two phases. The spatial derivatives are approximated by using some algebraic approximations such as finite difference, finite element, or radial basis functions (RBFs). This reduces the PDE to a system of ordinary differential equations (ODEs), which can be solved by any time-stepping method contrary to the application of the Laplace transform method [9] or space-time approach [10] where the time integration procedure preferably avoided. Many authors used the MOL approach for solving various engineering problems (see for example [11–16] and the references therein).
In this research, we use multiquadric (MQ), , as RBF for space discretization and the Runge-Kutta method of order 4 (RK4) for time integration. Kansa [17, 18] used MQ for solving problems in computational fluid dynamics. Since then, several authors have solved PDEs by using RBFs, e.g., [19–25]. Very little attention has been given to the RBFs in nonrectangular coordinates. Much of the work in nonrectangular coordinates have been carried out by Flyer and her co-authors [22, 26–29]. Aminataei and Mazarei have used RBF on the polar coordinates to solve Poisson’s equation [30]. In the MQ method, the parameter is called the shape parameter which shapes the stability and accuracy of the method. The selection of requires much attention as there is a tradeoff between the stability and accuracy of the approximation method [19]. The MMOL has advantage over the traditional mesh-based methods [31]. It saves the computational time and can be applied to problems with more complicated geometries. In the nonrectangular coordinate system, RBF is best suited because of its radial nature. In this work, we use RBF to find numerical solution of the heat equation in the polar cylindrical form. The study explores the richness of RBF in polar cylindrical coordinates.
Rest of the study is organized as follows: Section 2 presents formulation of the proposed method. Section 3 gives the numerical examples and discusses the results. Section 4 concludes the study.
2. Formulation of the Method
We consider the heat equation in polar cylindrical coordinateswhere , or 1. Here, we study the radially symmetric casewhere . The above equation can be written aswhere is the spatial differential operator. The initial condition (IC) iswhere or according to or 1, respectively. The boundary conditions (BCs) along areorand along ,where , , and are the functions of , or they are functions of and , when or 1, respectively.
2.1. Space Discretization Using RBF
We choose distinct nodes, , in the region , is the space dimension, for approximating a function by aswhere is some RBF, , and are the unknowns . Using the collocation points , we obtain
Applying on the above equation, we havewhere is the number of nodes at which the PDE is enforced. Similarly, let is the number of nodes at which BCs are enforced; then,where is the boundary operator. Writing equation (9) in matrix notation, we getsuch that
Micchelli [32] proved the invertibility of the system matrix . It follows from equations (10)–(12) thatwhere and are the submatrices of evaluated at the interior and boundary nodes, respectively. Let ; then, the above equation becomeswhere is called the differentiation matrix [33] which will be used for eigenvalue stability.
2.1.1. RBF in Polar Form
The cylindrical and Cartesian coordinates are related by
The distance between two points and is given bysubstituting the values of and , we obtain
For a more detailed discussion about the RBF in non-Cartesian coordinates, see [22].
2.1.2. Derivatives of RBF
The derivatives of the RBF, , with respect to the space variables and , are given bywhereand is the multiquadric (MQ) RBF. Applying RBF approximation to equations (3), (5)–(7), we have the following system of ODEs:with the corresponding ICs
2.2. Time Integration and Stability of the Proposed Numerical Scheme
In the current procedure, the time-dependent model is reduced to the ODEs system in time variable alone. The transformed coupled system of reduced ODEs is to be solved by any time integration method (e.g., Runge-Kutta). The stability of the proposed method can be discussed numerically by incorporating the rule of thumbs [34]. Consequently, the proposed method will be stable if for the operator in (21), the eigenvalues scaled by step size in time lie inside the region of stability of Runge-Kutta.
2.3. Shape Parameter and Its Selection
Most of the available basis functions of radial type contain a scaling factor also known as shape parameter which changes the shape of the basis functions and which is greatly related to solution accuracy. For better accuracy, we always need to use the optimal value of the scaling factor. Many criteria are available to find the optimal value of the RBF shape parameter. But in the current work, the criterion developed in [21] is used for finding the good value scaling factor for better accuracy.
3. Numerical Examples
In this section, we consider numerical examples of the heat equation in polar cylindrical coordinates. In all the examples, a radially symmetric case (independent of ) is considered. To validate the numerical solution, we compare the results with the exact solution and some other numerical methods from the literature. Convergence and eigenvalue stability is also discussed to further verify the numerical computations.
Example 1. Consider equation (3) when and the following IC [8]:The exact solution is [8]In numerical computations, is used. We solved the problem over the domain at . BCs were derived from the exact solution, equation (24). Numerical solution was obtained at different numbers of nodes and different values of the shape parameter to show the convergence of the method. We examined stability of the method in terms of eigenvalues of the differentiation matrix. Table 1 provides the condition number of the system matrix, spectral radius of the differentiation matrix, and error, [8] at different number of nodes when shape parameter . Table 1 also provides that the results are comparable with Li et al. [8]. Figure 1 shows the solution and absolute error at time , number of nodes , and shape parameter . Convergence and eigenvalue stability are shown in Figure 2, while Figure 3 shows the shape parameter with error and condition number. From Table 1 and Figure 2(a), it is evident that the method converges as the number of nodes increases. In Figure 2(b), we can see that all eigenvalues of the differentiation matrix lie within the stability region of RK4. Figure 3 shows that small error can be obtained when condition number is relatively large, but at the same time, eigenvalues must remain inside the stability domain.
(a)
(b)
(a)
(b)
(a)
(b)
Example 2. Next, we solve equation (3) when with the following IC [3]:where is the Bessel function of the first kind, and is the zero of .
The exact solution is [3]We obtained the numerical solution over the region at different time, . We utilized the exact solution (26) to derive the BCs. Various number of nodes and shape parameters were used in the computations. Table 2 compares the results with [3] (formula (2.1) of the reference). The table provides that results of the current method are better than the cited work. Table 3 and Figure 4(a) indicate the convergence of the method. Figure 5 displays the solution and absolute error of the method. It can be noted from Figure 4(b) that eigenvalues of are stable. Shape parameter and maximum absolute error and condition number are shown in Figure 6. Optimum value of the shape parameter lies between 7.1 and 7.2, as shown in Figure 6(a). In order to get a small error, conditioning has to be sacrificed to some extent.
(a)
(b)
(a)
(b)
(a)
(b)
Example 3. Now, we consider equation (3) when with the following IC [3]:The exact solution of this problem is [3]We used the domain to solve the problem numerically at different times. BCs were derived from the exact solution (28). Table 4 provides a comparison of maximum absolute error with Iyengar and Mittal [3]. The table provides that our result is a bit poor than [3], but still, it is quite reasonable as compared with the exact solution. Convergence and stability of the method are given in Table 5 and Figure 7. Solution and absolute error are shown in Figure 8. Shape parameter and maximum absolute error and condition number are shown in Figure 9. Error increases by decreasing the condition number. We must keep a balance between the error and the condition number. For a small error, we need an acceptable condition number. The shape parameter plays a significant role in the accuracy and stability. Figure 9(a) shows that the optimum value of is in the vicinity of 0.68.
(a)
(b)
(a)
(b)
(a)
(b)
Example 4. As a last example, we consider equation (3) when over the region , .
The exact solution of this problem is given by [5]where and . The IC and BCs were derived from the exact solution (29). Table 6 provides the condition number, spectral radius, maximum absolute error, and the error at different numbers of nodes when . Error decreases with the increase of nodes. We have a reasonable error even for a small number of nodes. Figure 10 plots the solution and absolute error. Table 6 and Figure 11 show the convergence and stability. Figure 12 shows the shape parameter and maximum absolute error and condition number of the system matrix.
(a)
(b)
(a)
(b)
(a)
(b)
4. Conclusion
The heat equation is one of the fundamental PDE models as it captures several physical phenomena ranging from temperature of a plate to the heat transfer in tumors. We solved the equation in polar cylindrical form which is a natural coordinate system for problems with the circular domain. Furthermore, the radial characteristic of RBFs makes them more suitable in polar form. Here, we used the MQ–RBF, which has a shape parameter. The shape parameter plays a vital role in the stability and accuracy of the numerical solution. The parameter has to be selected in such a way that the accuracy and stability are maintained simultaneously. The proposed method is meshless in nature and thus can effectively be applied to problems with irregular domain. The present approach is an efficient addition to the existing methods for the heat equation in polar cylindrical coordinates.
Data Availability
The data used to support this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.