Numerical Solution of Heat Equation in Polar Cylindrical Coordinates by the Meshless Method of Lines

We propose a numerical solution to the heat equation in polar cylindrical coordinates by using the meshless method of lines approach.-e 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. -is work is an attempt to explore the versatility of RBFs in nonrectangular coordinates as well.-e 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.


Introduction
e 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. ey 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. ey 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. ese methods do not require mesh generation, and thus, they are more flexible and save computational time. e method of lines (MOL) is a strong numerical technique that solves a partial differential equation (PDE) in two phases.
e spatial derivatives are approximated by using some algebraic approximations such as finite difference, finite element, or radial basis functions (RBFs). is 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 spacetime approach [10] where the time integration procedure preferably avoided. Many authors used the MOL approach for solving various engineering problems (see for example [11][12][13][14][15][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][20][21][22][23][24][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][27][28][29]. Aminataei and Mazarei have used RBF on the polar coordinates to solve Poisson's equation [30]. In the MQ method, the parameter c is called the shape parameter which shapes the stability and accuracy of the method. e selection of c requires much attention as there is a tradeoff between the stability and accuracy of the approximation method [19]. e 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. e 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.

Formulation of the Method
We consider the heat equation in polar cylindrical coordinates where u � u(ρ, θ, z, t), a � 0 or 1. Here, we study the radially symmetric case where u � u(ρ, z, t). e above equation can be written as where L � (z 2 /zρ 2 ) + (1/ρ)(z/zρ) + a(z 2 /zz 2 ) is the spatial differential operator. e initial condition (IC) is where u 0 � u 0 (ρ) or u 0 � u 0 (ρ, z) according to a � 0 or 1, respectively. e boundary conditions (BCs) along ρ are or and along z, where h 1 , h 2 , and h 3 are the functions of t, or they are functions of z and t, when a � 0 or 1, respectively.

Space Discretization Using
where are the unknowns (i � 1, 2, . . . , N). Using the collocation points x j , we obtain Applying L on the above equation, we have where N I is the number of nodes at which the PDE is enforced. Similarly, let N B is the number of nodes at which BCs are enforced; then, where B is the boundary operator. Writing equation (9) in matrix notation, we get such that 2 Journal of Mathematics Micchelli [32] proved the invertibility of the system matrix A. It follows from equations (10)-(12) that where A I and A B are the submatrices of A evaluated at the interior and boundary nodes, respectively. Let where M is called the differentiation matrix [33] which will be used for eigenvalue stability.

RBF in Polar
Form. e cylindrical and Cartesian coordinates are related by e distance between two points x � (x, y, z) and substituting the values of x and y, we obtain For a more detailed discussion about the RBF in non-Cartesian coordinates, see [22].

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 t alone. e transformed coupled system of reduced ODEs is to be solved by any time integration method (e.g., Runge-Kutta). e 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 M in (21), the eigenvalues scaled by step size δt in time t lie inside the region of stability of Runge-Kutta.

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.

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 a � 0 and the following IC [8]: e exact solution is [8] u In numerical computations, τ 0 � 0.1 is used. We solved the problem over the domain 0 ≤ ρ ≤ 3 at t � 0.8. 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 L 1 error, at different number of nodes when shape parameter c � 8.86. Table 1 also provides that the results are comparable with Li et al. [8]. Figure 1 shows the solution and absolute error at time t � 0.8, number of nodes N � 1000, and shape parameter c � 3.25. Convergence and eigenvalue stability are shown in Figure 2, while Figure 3 shows the shape parameter with L 1 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.
Example 2. Next, we solve equation (3) when a � 0 with the following IC [3]: where J 0 is the Bessel function of the first kind, and β is the zero of J 0 . e exact solution is [3] u(ρ, t) � J 0 (βρ)e − β 2 t .
We obtained the numerical solution over the region 0 ≤ ρ ≤ 1 at different time, t. 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). e 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 M are stable. Shape parameter and maximum absolute error and condition number are shown in Figure 6. Optimum value of the shape parameter c 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.
Example 3. Now, we consider equation (3) when a � 0 with the following IC [3]: e exact solution of this problem is [3] u(ρ, t) � J 0 (ρ)e − t . (28) We used the domain 2 ≤ ρ ≤ 3 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]. e 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. e shape parameter c plays a significant role in the accuracy and stability. Figure 9(a) shows that the optimum value of c is in the vicinity of 0.68. e exact solution of this problem is given by [5]

Journal of Mathematics
where Z � π(z − 1)/2 and ξ � 1 + π 2 /4. e IC and BCs were derived from the exact solution (29). Table 6 provides the condition number, spectral radius, maximum absolute error, and the L 1 error at different numbers of nodes when t � 1.6.
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.  Journal of Mathematics

Conclusion
e 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. e shape parameter plays a vital role in the stability and accuracy of the numerical solution. e parameter c has to be selected in such a way that the accuracy and stability are maintained simultaneously. e proposed method is meshless in nature and thus can effectively be applied to problems with irregular domain. e present approach is an efficient addition to the existing methods for the heat equation in polar cylindrical coordinates.

Data Availability
e data used to support this study are included within the article.