Linear and Nonlinear Matrix Equations: Solutions, Convergence, Perturbation, and OptimizationView this Special Issue
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. 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.
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  used the Crank–Nicolson method to solve the cylindrical heat equation. Mitchell and Pearce  applied the explicit finite difference method to obtain numerical solution of the cylindrical heat conduction equation. Iyengar and Mittal  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  to solve the heat equation. Iyengar and Manohar  developed high-order difference methods for the solution of heat equation in polar cylindrical form. Manohar et al.  derived two-level, three-point difference schemes to solve the heat equation with linear variable coefficient. Al-Odat and Al-Nimr  studied thermal stability of composite superconducting tape by considering the heat conduction equation. Li et al.  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  or space-time approach  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 . 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 . The MMOL has advantage over the traditional mesh-based methods . 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  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  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 .
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 . 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  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 :The exact solution is 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,  at different number of nodes when shape parameter . Table 1 also provides that the results are comparable with Li et al. . 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.
Example 2. Next, we solve equation (3) when with the following IC :where is the Bessel function of the first kind, and is the zero of .
The exact solution is 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  (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.
Example 3. Now, we consider equation (3) when with the following IC :The exact solution of this problem is 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 . The table provides that our result is a bit poor than , 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.
Example 4. As a last example, we consider equation (3) when over the region , .
The exact solution of this problem is given by 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.
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.
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.
E. L. Albasiny, “On the numerical solution of a cylindrical heat-conduction problem,” Quarterly Journal of Mechanics & Applied Mathematics, vol. 13, no. 3, pp. 374–384, 1960.View at: Publisher Site | Google Scholar
A. R. Mitchell and R. P. Pearce, “Explicit difference methods for solving the cylindrical heat conduction equation,” Mathematics of Computation, vol. 17, no. 84, pp. 426–432, 1963.View at: Publisher Site | Google Scholar
S. R. K. Iyengar and R. C. Mittal, “High accuracy difference schemes for the cylindrical heat conduction equation,” IMA Journal of Applied Mathematics, vol. 22, no. 3, pp. 321–330, 1978.View at: Publisher Site | Google Scholar
E. Livne and A. Glasner, “A finite difference scheme for the heat conduction equation,” Journal of Computational Physics, vol. 58, no. 1, pp. 59–66, 1985.View at: Publisher Site | Google Scholar
S. R. K. Iyengar and R. Manohar, “High order difference methods for heat equation in polar cylindrical coordinates,” Journal of Computational Physics, vol. 77, no. 2, pp. 425–438, 1988.View at: Publisher Site | Google Scholar
R. Manohar, S. R. K. Iyengar, and U. A. Krishnaiah, “High order difference methods for linear variable coefficient parabolic equation,” Journal of Computational Physics, vol. 77, no. 2, pp. 513–523, 1988.View at: Publisher Site | Google Scholar
M. Q. Al-Odat and M. A. Al-Nimr, “Thermal stability of composite superconducting tape under the effect of a two-dimensional dual-phase-lag heat conduction model,” Heat and Mass Transfer, vol. 40, no. 3-4, pp. 211–217, 2004.View at: Publisher Site | Google Scholar
H. Li, X. Wu, and J. Zhang, “Local artificial boundary conditions for Schrödinger and heat equations by using high-order azimuth derivatives on circular artificial boundary,” Computer Physics Communications, vol. 185, no. 6, pp. 1606–1615, 2014.View at: Publisher Site | Google Scholar
M. Uddin, A. Kamran, and A. Ali, “A localized transform-based meshless method for solving time fractional wave-diffusion equation,” Engineering Analysis with Boundary Elements, vol. 92, pp. 108–113, 2018.View at: Publisher Site | Google Scholar
M. Uddin, H. Ali, and H. Ali, “Space-time kernel based numerical method for generalized black-scholes equation,” Discrete & Continuous Dynamical Systems-S, vol. 13, no. 10, pp. 2905–2915, 2020.View at: Publisher Site | Google Scholar
S. S. Hu and W. E. Schiesser, “An adaptive grid method in the numerical method of lines,” in Advances in Computer Methods for Partial Differential Equations, R. Vichnevetsky and R. Stepleman, Eds., pp. 305–311, Pergamon Press, Oxford, UK, 1981.View at: Google Scholar
W. E. Schiesser, The Numerical Method of Lines: Integration of Partial Differential Equations, Academic Press, San Diego, CA, USA, 1991.
S. Hamdi, W. H. Enright, Y. Ouellet, and W. E. Schiesser, “Method of lines solutions of the extended Boussinesq equations,” Journal of Computational and Applied Mathematics, vol. 183, no. 2, pp. 327–342, 2005.View at: Publisher Site | Google Scholar
P. Saucez, A. Vande Wouwer, W. E. Schiesser, and P. Zegeling, “Method of lines study of nonlinear dispersive waves,” Journal of Computational and Applied Mathematics, vol. 168, no. 1-2, pp. 413–423, 2004.View at: Publisher Site | Google Scholar
A. V. Wouwer, P. Saucez, W. E. Schiesser, and S. Thompson, “A Matlab implementation of upwind finite differences and adaptive grids in the method of lines,” Journal of Computational and Applied Mathematics, vol. 183, no. 2, pp. 245–258, 2005.View at: Publisher Site | Google Scholar
G. W. Griffiths and W. E. Schiesser, Traveling Wave Analysis of Partial Differential Equations: Numerical and Analytical Methods with Matlab and Maple, Academic Press, Cambridge, MA, USA, 2012.
E. J. Kansa, “Multiquadrics-A scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates,” Computers & Mathematics with Applications, vol. 19, no. 8-9, pp. 127–145, 1990.View at: Publisher Site | Google Scholar
E. J. Kansa, “Multiquadrics-A scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations,” Computers & Mathematics with Applications, vol. 19, no. 8-9, pp. 147–161, 1990.View at: Publisher Site | Google Scholar
R. Schaback, “Error estimates and condition numbers for radial basis function interpolation,” Advances in Computational Mathematics, vol. 3, no. 3, pp. 251–264, 1995.View at: Publisher Site | Google Scholar
B. Fornberg, T. A. Driscoll, G. Wright, and R. Charles, “Observations on the behavior of radial basis function approximations near boundaries,” Computers & Mathematics with Applications, vol. 43, no. 3-5, pp. 473–490, 2002.View at: Publisher Site | Google Scholar
G. E. Fasshauer and J. G. Zhang, “On choosing “optimal” shape parameters for RBF approximation,” Numerical Algorithms, vol. 45, no. 1-4, pp. 345–368, 2007.View at: Publisher Site | Google Scholar
N. Flyer and G. B. Wright, “Transport schemes on a sphere using radial basis functions,” Journal of Computational Physics, vol. 226, no. 1, pp. 1059–1084, 2007.View at: Publisher Site | Google Scholar
S. Haq, A. Hussain, and M. Uddin, “On the numerical solution of nonlinear Burgers’-type equations using meshless method of lines,” Applied Mathematics and Computation, vol. 218, no. 11, pp. 6280–6290, 2012.View at: Publisher Site | Google Scholar
A. Hussain, S. Haq, and M. Uddin, “Numerical solution of Klein-Gordon and sine-Gordon equations by meshless method of lines,” Engineering Analysis with Boundary Elements, vol. 37, no. 11, pp. 1351–1366, 2013.View at: Publisher Site | Google Scholar
S. Ul-Islam, N. Haider, and I. Aziz, “Meshless and multi-resolution collocation techniques for parabolic interface models,” Applied Mathematics and Computation, vol. 335, pp. 313–332, 2018.View at: Publisher Site | Google Scholar
N. Flyer and G. B. Wright, “A radial basis function method for the shallow water equations on a sphere,” Proceedings of the Royal Society A: Mathematical, Physical & Engineering Sciences, vol. 465, no. 2106, pp. 1949–1976, 2009.View at: Publisher Site | Google Scholar
N. Flyer and B. Fornberg, “Radial basis functions: developments and applications to planetary scale flows,” Computers & Fluids, vol. 46, no. 1, pp. 23–32, 2011.View at: Publisher Site | Google Scholar
V. Bayona, N. Flyer, G. M. Lucas, and A. J. G. Baumgaertner, “A 3-D RBF-FD solver for modeling the atmospheric global electric circuit with topography (GEC-RBFFD v1.0),” Geoscientific Model Development, vol. 8, no. 10, pp. 3007–3020, 2015.View at: Publisher Site | Google Scholar
Y. Gu, C.-M. Fan, and Z. Fu, “Localized method of fundamental solutions for three-dimensional elasticity problems: Theory,” Advances in Applied Mathematics and Mechanics, vol. 13, no. 6, 2021.View at: Google Scholar
A. Aminataei and M. M. Mazarei, “Numerical solution of Poisson’s equation using radial basis function networks on the polar coordinate,” Computers & Mathematics with Applications, vol. 56, no. 11, pp. 2887–2895, 2008.View at: Publisher Site | Google Scholar
H. Xia and Y. Gu, “Generalized finite difference method for electroelastic analysis of three-dimensional piezoelectric structures,” Applied Mathematics Letters, vol. 117, Article ID 107084, 2021.View at: Publisher Site | Google Scholar
C. A. Micchelli, “Interpolation of scattered data: distance matrices and conditionally positive definite functions,” Constructive Approximation, vol. 2, no. 1, pp. 11–22, 1986.View at: Publisher Site | Google Scholar
S. A. Sarra, “A numerical study of the accuracy and stability of symmetric and asymmetric RBF collocation methods for hyperbolic PDEs,” Numerical Methods for Partial Differential Equations, vol. 24, no. 2, pp. 670–686, 2008.View at: Publisher Site | Google Scholar
M. Uddin, “RBF-PS scheme for solving the equal width equation,” Applied Mathematics and Computation, vol. 222, pp. 619–631, 2013.View at: Publisher Site | Google Scholar