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

Authors: Arshad Hussain; Marjan Uddin; Sirajul Haq; Hameed Ullah Jan
Journal: Journal of Mathematics
 (2021)

---

## 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.


---

## Body

## 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),ϒr=1+c2r2, 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 c is called the shape parameter which shapes the stability and accuracy of the method. The selection of c 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: Section2 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 coordinates(1)∂u∂t=∂2u∂ρ2+1ρ∂u∂ρ+1ρ2∂2u∂θ2+a∂2u∂z2,ρ0≤ρ≤ρ1,0≤θ≤2π,z0≤z≤z1,t>0,where u=uρ,θ,z,t, a=0 or 1. Here, we study the radially symmetric case(2)∂u∂t=∂2u∂ρ2+1ρ∂u∂ρ+a∂2u∂z2,ρ0≤ρ≤ρ1,z0≤z≤z1,t>0,where u=uρ,z,t. The above equation can be written as(3)∂u∂t=ℒu,where ℒ=∂2/∂ρ2+1/ρ∂/∂ρ+a∂2/∂z2 is the spatial differential operator. The initial condition (IC) is(4)u=u0,t=0,where u0=u0ρ or u0=u0ρ,z according to a=0 or 1, respectively. The boundary conditions (BCs) along ρ are(5)∂u∂ρ=0,ρ=ρ0=0,u=h1,ρ=ρ1,or(6)u=h2,ρ=ρ0,u=h3,ρ=ρ1,and along z,(7)u=g1ρ,t,z=z0,u=g2ρ,t,z=z1,where h1, h2, and h3 are the functions of t, or they are functions of z and t, when a=0 or 1, respectively.


### 2.1. Space Discretization Using RBF

We chooseN distinct nodes, xi, in the region Ω⊂ℝd, d is the space dimension, for approximating a function ux,t by Ux,t as(8)Ux,t=∑i=1Nbitϒx−xi=∑i=1Nbitϒir,where ϒir=ϒx−xi is some RBF, r=x−xi, and bit are the unknowns i=1,2,…,N. Using the collocation points xj, we obtain(9)Uxj,t=∑i=1Nbitϒirj,j=1,2,…,N.Applyingℒ on the above equation, we have(10)ℒUxj,t=∑i=1Nbitℒϒirj,j=1,2,…,NI,where NI is the number of nodes at which the PDE is enforced. Similarly, let NB is the number of nodes at which BCs are enforced; then,(11)ℬUxj,t=∑i=1Nbitℬϒirj,j=NI+1,NI+2,…,N,where ℬ is the boundary operator. Writing equation (9) in matrix notation, we get(12)U=Ab,such that(13)U=U1,U2,…,UNT,Uj=Uxj,t,forj=1,2,…,N,A=ϒ1r1ϒ2r1⋯ϒNr1ϒ1r2ϒ2r2⋯ϒNr2⋮⋮⋱⋮ϒ1rNϒ2rN⋯ϒNrN,b=b1b2⋮bN.Micchelli [32] proved the invertibility of the system matrix A. It follows from equations (10)–(12) that(14)ℒℬU=ℒAIℬABA−1U,where AI and AB are the submatrices of A evaluated at the interior and boundary nodes, respectively. Let ℒAIℬABA−1=ℳ; then, the above equation becomes(15)ℒℬU=ℳU,where ℳ 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(16)x=ρcosθ,y=ρsinθ,z=z,ρ=x2+y2,θ=tan−1yx.The distance between two pointsx=x,y,z and xi=xi,yi,zi is given by(17)r=x−xi=x−xi2+y−yi2+z−zi2,substituting the values of x and y, we obtain(18)r=ρ2+ρi2−2ρρicosθ−θi+z−zi2.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,ϒr, with respect to the space variables ρ and z, are given by(19)∂ϒ∂ρ=∂r∂ρdϒdr,∂2ϒ∂ρ2=∂2r∂ρ2dϒdr+∂r∂ρ2d2ϒdr2,∂2ϒ∂z2=∂2r∂z2dϒdr+∂r∂z2d2ϒdr2,where(20)∂r∂ρ=ρ−ρicosθ−θir,∂2r∂ρ2=ρ2−ρ−ρicosθ−θi2r3,∂r∂z=z−zir,∂2r∂z2=1−∂r/∂z2r,and ϒ=ϒr=1+cr2 is the multiquadric (MQ) RBF. Applying RBF approximation to equations (3), (5)–(7), we have the following system of ODEs:(21)dUdt=ℳU,with the corresponding ICs(22)U0=u0.


### 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 variablet 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 δt in time t 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.


## 2.1. Space Discretization Using RBF

We chooseN distinct nodes, xi, in the region Ω⊂ℝd, d is the space dimension, for approximating a function ux,t by Ux,t as(8)Ux,t=∑i=1Nbitϒx−xi=∑i=1Nbitϒir,where ϒir=ϒx−xi is some RBF, r=x−xi, and bit are the unknowns i=1,2,…,N. Using the collocation points xj, we obtain(9)Uxj,t=∑i=1Nbitϒirj,j=1,2,…,N.Applyingℒ on the above equation, we have(10)ℒUxj,t=∑i=1Nbitℒϒirj,j=1,2,…,NI,where NI is the number of nodes at which the PDE is enforced. Similarly, let NB is the number of nodes at which BCs are enforced; then,(11)ℬUxj,t=∑i=1Nbitℬϒirj,j=NI+1,NI+2,…,N,where ℬ is the boundary operator. Writing equation (9) in matrix notation, we get(12)U=Ab,such that(13)U=U1,U2,…,UNT,Uj=Uxj,t,forj=1,2,…,N,A=ϒ1r1ϒ2r1⋯ϒNr1ϒ1r2ϒ2r2⋯ϒNr2⋮⋮⋱⋮ϒ1rNϒ2rN⋯ϒNrN,b=b1b2⋮bN.Micchelli [32] proved the invertibility of the system matrix A. It follows from equations (10)–(12) that(14)ℒℬU=ℒAIℬABA−1U,where AI and AB are the submatrices of A evaluated at the interior and boundary nodes, respectively. Let ℒAIℬABA−1=ℳ; then, the above equation becomes(15)ℒℬU=ℳU,where ℳ 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(16)x=ρcosθ,y=ρsinθ,z=z,ρ=x2+y2,θ=tan−1yx.The distance between two pointsx=x,y,z and xi=xi,yi,zi is given by(17)r=x−xi=x−xi2+y−yi2+z−zi2,substituting the values of x and y, we obtain(18)r=ρ2+ρi2−2ρρicosθ−θi+z−zi2.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,ϒr, with respect to the space variables ρ and z, are given by(19)∂ϒ∂ρ=∂r∂ρdϒdr,∂2ϒ∂ρ2=∂2r∂ρ2dϒdr+∂r∂ρ2d2ϒdr2,∂2ϒ∂z2=∂2r∂z2dϒdr+∂r∂z2d2ϒdr2,where(20)∂r∂ρ=ρ−ρicosθ−θir,∂2r∂ρ2=ρ2−ρ−ρicosθ−θi2r3,∂r∂z=z−zir,∂2r∂z2=1−∂r/∂z2r,and ϒ=ϒr=1+cr2 is the multiquadric (MQ) RBF. Applying RBF approximation to equations (3), (5)–(7), we have the following system of ODEs:(21)dUdt=ℳU,with the corresponding ICs(22)U0=u0.


## 2.1.1. RBF in Polar Form

The cylindrical and Cartesian coordinates are related by(16)x=ρcosθ,y=ρsinθ,z=z,ρ=x2+y2,θ=tan−1yx.The distance between two pointsx=x,y,z and xi=xi,yi,zi is given by(17)r=x−xi=x−xi2+y−yi2+z−zi2,substituting the values of x and y, we obtain(18)r=ρ2+ρi2−2ρρicosθ−θi+z−zi2.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,ϒr, with respect to the space variables ρ and z, are given by(19)∂ϒ∂ρ=∂r∂ρdϒdr,∂2ϒ∂ρ2=∂2r∂ρ2dϒdr+∂r∂ρ2d2ϒdr2,∂2ϒ∂z2=∂2r∂z2dϒdr+∂r∂z2d2ϒdr2,where(20)∂r∂ρ=ρ−ρicosθ−θir,∂2r∂ρ2=ρ2−ρ−ρicosθ−θi2r3,∂r∂z=z−zir,∂2r∂z2=1−∂r/∂z2r,and ϒ=ϒr=1+cr2 is the multiquadric (MQ) RBF. Applying RBF approximation to equations (3), (5)–(7), we have the following system of ODEs:(21)dUdt=ℳU,with the corresponding ICs(22)U0=u0.


## 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 variablet 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 δt in time t 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 a=0 and the following IC [8]:(23)u0ρ=14πτ0exp−ρ24τ0.
The exact solution is [8](24)uρ,t=14πt+τ0exp−ρ24t+τ0.
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 L1 error, L1=1/N∑i=1Nui−Ui [8] 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 L1 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.Table 1
Convergence of the method corresponding to Example1.
NCondition numberSpectral radiusL1 errorMQ[8]10002.44071e + 0050.388561.42044e − 0041.029e − 004 (32 × 32)a20001.81550e + 0061.036114.17586e − 0052.715e − 005 (64 × 64)30006.93409e + 0061.847693.16201e − 0057.130e − 006 (128 × 128)aMesh size.Figure 1
Example1: solution and absolute error at t=0.8,N=1000, and c=3.25. (a) Exact solution Δ, numerical solution •. (b) Absolute error.
(a)(b)Figure 2
Example1: convergence and eigenvalue stability at t=0.8. (a) Convergence, c=8.86. (b) Scaled eigenvalues, N=1000, c=3.25.
(a)(b)Figure 3
Example1: shape parameter, condition number, and L1 error at t=0.8,N=1000. (a) Shape parameter and L1 error. (b) Shape parameter and condition number.
(a)(b)Example 2.
Next, we solve equation (3) when a=0 with the following IC [3]:(25)u0ρ=J0βρ,where J0 is the Bessel function of the first kind, and β is the zero of J0.
The exact solution is [3](26)uρ,t=J0βρe−β2t.
We obtained the numerical solution over the region0≤ρ≤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). 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 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.Table 2
Maximum absolute error whenN=300, c=7.14 corresponding to Example 2.
TimeMQ[3]1.52.10044e − 0074.4900e − 0063.03.47126e − 0111.5100e − 0093.754.53141e − 0132.4500e − 011Table 3
Convergence of the method corresponding to Example2.
NCondition numberSpectral radiusMaximum absolute erroraL1 error3005.61645e + 0041.048734.08133e − 0122.24683e − 0124001.33014e + 0051.790491.69287e − 0121.04455e − 0125002.66429e + 0052.329086.64945e − 0134.55461e − 013aMaximum absolute error.Figure 4
Example2: convergence and eigenvalue stability at t=3.75. (a) Convergence, c=11.85. (b) Scaled eigenvalues, N=300, c=7.14.
(a)(b)Figure 5
Example2: solution and absolute error at t=3.75,N=300, and c=7.14. (a) Exact solution Δ, numerical solution •. (b) Absolute error.
(a)(b)Figure 6
Example2: shape parameter, condition number, and maximum absolute error at t=3.75,N=300. (a) Shape parameter and maximum absolute error. (b) Shape parameter and condition number.
(a)(b)Example 3.
Now, we consider equation (3) when a=0 with the following IC [3]:(27)u0ρ=J0ρ.
The exact solution of this problem is [3](28)uρ,t=J0ρe−t.
We used the domain2≤ρ≤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]. 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 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.Table 4
Maximum absolute error whenN=300, c=0.68 corresponding to Example 3.
TimeMQ[3]1.56.31429e − 0052.9500e − 0063.01.42345e − 0056.5900e − 0073.756.72931e − 0063.1100e − 007Table 5
Convergence of the method corresponding to Example3.
NCondition numberSpectral radiusMaximum absolute errorL1 error3001.87647e + 0070.166862.21488e − 0051.01050e − 0056003.93730e + 0090.435399.46609e − 0066.47678e − 0069001.47083e + 0110.451585.22090e − 0065.78701e − 006Figure 7
Example3: convergence and eigenvalue stability at t=3.75. (a) Convergence, c=1.62. (b) Scaled eigenvalues, N=300, c=0.68.
(a)(b)Figure 8
Example3: solution and absolute error at t=3.75,N=300, and c=0.68. (a) Exact solution Δ, numerical solution •. (b) Absolute error.
(a)(b)Figure 9
Example3: shape parameter, condition number, and maximum absolute error at t=3.75,N=300. (a) Shape parameter and maximum absolute error. (b) Shape parameter and condition number.
(a)(b)Example 4.
As a last example, we consider equation (3) when a=1 over the region 0≤ρ≤1, −1≤z≤1.
The exact solution of this problem is given by [5](29)uρ,z,t=J0ρcosZ+sinZe−ξt,where Z=πz−1/2 and ξ=1+π2/4. The IC and BCs were derived from the exact solution (29). Table 6 provides the condition number, spectral radius, maximum absolute error, and the L1 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.Table 6
Convergence of the method corresponding to Example4.
NCondition numberSpectral radiusMaximum absolute errorL1 error1211.19944e + 0040.134801.43233e − 0043.76819e − 0052526.74837e + 0050.497373.31975e − 0058.09290e − 0069611.99564e + 0071.104141.24169e − 0052.52067e − 006Figure 10
Example4: solution and absolute error at t=1.6,N=121, and c=2.53. (a) Exact and numerical solution. (b) Absolute error.
(a)(b)Figure 11
Example4: convergence and eigenvalue stability at t=1.6. (a) Convergence, c=13.78. (b) Scaled eigenvalues, N=121, c=2.53.
(a)(b)Figure 12
Example4: shape parameter, condition number, and maximum absolute error at t=1.6,N=121. (a) Shape parameter and maximum absolute error. (b) Shape parameter and condition number.
(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 parameterc 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.



---

*Source: 8862139-2021-12-13.xml*