Solving the Nonlinear Heat Equilibrium Problems Using the Local Multiquadric Radial Basis Function Collocation Method

: In this article, the nonlinear heat equilibrium problems are solved by the local multiquadric (MQ) radial basis function (RBF) collocation method. The system of nonlinear algebraic equations is solved by iteration based on the residual norm-based algorithm, in which the direction of evolution is determined by a linear equation. In addition, the role of the collocation point and source point is clearly deﬁned such that in our proposed method the ﬁeld value of any interested point can be expressed. Six numerical examples are shown to check the performance of the proposed method. As the number of supporting points (m p ) increases, the accuracy of numerical solution increases. Among all examples, m p = 50 can perform well. In addition, the selection of shape parameter, c, a ﬀ ects the accuracy. However, as c < 2 the maximum relative absolute error percentage is less than 1%.


Introduction
The nonlinear heat equilibrium problems can be written as: where k(T) is the temperature-dependent thermal conductivity, and T is the temperature field. This equation is a nonlinear one, and the nonlinearity occurs since the thermal conductivity is temperature-dependent. There exist many applications of using such materials, such as the nanofluid [1], single-walled carbon nanotubes [2], and porous materials [3].
To solve the abovementioned nonlinear partial differential equation, one may use the Kirchhoff transform to convert the nonlinear partial differential equation (PDE) into a linear one. If the Kirchhoff transform is not adopted, the governing equation remains nonlinear and a nonlinear solver is required. If the problem is a well-posed one, usually the Newton-Raphson method [4] convergences very fast. While the system is not a well-posed one, there are many nonlinear iteration methods such as the Landweber iteration method [5], the fictitious time integration method [6,7], the Jacobian inverse free method [8], the exponentially convergent scalar homotopy method [9], the double iteration process [10,11], and so on.
The governing equation, as well as its associate boundary value, needs to be discretized first. The discretization method may influence the accuracy of the solution as well as the computational efficiency. Generally speaking, most discretization methods fall into two groups: the mesh-dependent method and meshless method. For the mesh-dependent methods, the finite element method [12], the boundary element method [13], and the finite difference method [14] are typical examples. However, it has been reported that a lot of computation effort is consumed in mesh generation for these mesh-dependent methods. Based on the abovementioned reason, the meshless methods then were developed. For the meshless methods, the global methods were first developed. The meshless methods represent the solution as the superposition of the shape functions at many central points (or source points). The radial basis functions (RBF) [15], the method of fundamental solutions (MFS) [16], the Trefftz collocation methods (TCM) [17] are typical examples. The global meshless methods have some drawbacks such as dense matrix and ill-conditioned coefficient matrix which makes them not so popular.
To overcome the drawbacks for the global meshless methods, the local meshless methods then were developed, such as the local RBF method [18][19][20][21][22], the generalized finite difference method [23][24][25], and so on. The local meshless methods usually have a sparse coefficient matrix for a linear problem and the coefficient matrix is more well-posed than that resulting from the global meshless methods.
In this paper, the nonlinear heat equilibrium problems are solved by the local multiquadric radial basis functions (MQ RBF). It is difficult to compare the performance of various numerical techniques since most numerical algorithms have their merits as well as drawbacks. However, the current approach can have some merits over other techniques as stated in the following. While it is compared with the mesh-dependent methods, the current approach which is a meshless method such that it does not require mesh generation and thus can save computation efforts in mesh generation. In comparison with the boundary type methods such as the boundary element method or method of fundamental solution, the current approach, although it does not have the merit of reducing dimensionality, it does not, however, need to perform the inversion of the Kirchhoff transform, which sometimes becomes difficult when the integration of k(T) is not easy to be found. In comparison with the global radial basis function collocation method, the current approach does not encounter a full matrix form for the Jacobian matrix, and this merit makes it possible for us to increase the number of collocation points. Even this localized MQ RBF method is compared with other localized RBF such as localized Gaussian RBF, the current approach is more accurate since multiquadric RBF is recognized as the most accurate RBF. While it is compared with other localized methods such as the generalized finite difference method, the current approach still has better performance in accuracy under the same number of supporting points.
Aside from the introduction section, the rest of the contents of this paper are arranged as the following. In Section two, the mathematical backgrounds, such as the mathematical modeling of the nonlinear heat equilibrium problems, the local MQ RBF, and the residual norm-based algorithm [26], are introduced. Six examples are then given to show the validity of current approach in Section 3. In Section 4, the conclusions are drawn based on the results and findings.

Mathematical Modeling
The governing equation of the nonlinear heat equilibrium problem is stated in Equation (1). To solve this partial differential equation, the boundary conditions are required for points on the boundary: where x is the position vector, Γ D is the boundary where the Dirichlet type boundary condition is given, Γ N is the boundary where the Neumann type boundary condition is given, Γ m is the boundary where the nonlinear mixed type boundary condition is given, n denoted the outnormal direction, σ is the Stefan-Boltzmann constant, T f is the temperature of the surrounding medium, ε is the temperature-dependent emissivity between the surface Γ m and radiating medium at temperature T s , h is a constant, and T, q, and f are prescribed functions. For the standard boundary value problem, the union of Γ D , Γ N , and Γ m is the entire boundary and the intersection of any two of Γ D , Γ D , Γ N , and Γ m is the empty set. This means for each boundary point, only one type of boundary condition is given and for the entire boundary the boundary condition should be prescribed. The standard boundary value problem consists a system constructed by Equations (1)-(4).

Local MQ RBF
The local MQ RBF can be described as follows. First, for a collocation point x the physical quantity at this point (for example, the temperature field T) can be expressed as the linear combination of the MQ RBF shape functions for m p source points (or central points) in the supporting neighborhood of x (which is called "star" in some literature). The expression is written as where s j is the j-th source points in the supporting neighborhood, a j is the density of the shape function of the j-th source point, j is the shape function of the s j and can be written as in which c is the shape parameter and d is the maximum distance for any two points belonging to the m p source points.
In this study, we distribute source points on the boundary and domain uniformly. While the collocation point is determined, we select m p points among these source points by assigning the distance between collocation point and the candidate source point is ranked in the m p -th shortest ones. Notice that, the collocation points can be not in the set of source points. In previous published literature, the roles of the collocation point and source point are not clearly explained, and most researchers select the collocation points from the set of source points. Consequently, it is not clear how to calculate the filed quantity when the point considered is not in the set of source points. We should make a comment here that theoretically the set of source points is not necessarily constructed by the boundary points and internal domain points. For example, the local MFS [26] defines that for each collocation point the source points distribute on a circle surrounding the collocation point.
When the m p source points for the collocation point are selected, the local MQ RBF wants to express the source density a j first. Now assume that for these m p source points, the field quantities at these points can be expressed by Equation (5) as well. That means we have where T is a column vector consisting the temperature values for these m p points, i.e., T = T(s 1 ) · · · T(s m p ) T ; a is a column vector consisting the densities of the shape functions, i.e., a = a 1 · · · a m p T ; and the H is a m p by m p matrix written as in which Φ j s i , s j is the i-th row vector inside H. From Equation (5), the source density can be obtained by After the densities of shape functions are determined, the value of field quantity can be calculated by inserting Equation (9) into Equation (5): where h = Φ 1 (x, s 1 ) · · · Φ m p x, s m p is a row vector with the shape functions of m p source points. Considering a two-dimensional domain and the x-y Cartesian coordinate system is adopted, the derivatives of T can be expressed from Equation (10) such as T x (x) (the derivative of T with respect to x) is expressed as where ∂Φ mp x,s mp ∂x Therefore, for the internal points inside the domain the governing equation (Equation (1)) is required to be satisfied and for the boundary points one of the boundary conditions (Equations (2)-(4)) should be satisfied. For each point (internal point or boundary point), one equation is constructed. Collecting all equations, the following system of algebraic equations is constructed: in which F represents the vector consisting the residuals for all equations, T is the temperature field distribution vector. Notice that each equation in the system of Equation (12) is constructed at a specific collocation point which relates to m p supporting source points only such that the Jacobian matrix B = ∂F ∂T is a sparse matrix that is beneficial for a large-scale problem with many collocation points. This algebraic equation system is nonlinear, and solving it requires a nonlinear solver. In this study, the residual norm-based algorithm (RNBA) is adopted and it will be briefly introduced in the next subsection.

Residual Norm Based Algorithm
The details of the derivation of residual norm-based algorithm can be found in [26]. It begins with the following scalar equation, which is equivalent to Equation (12): Now, let us construct a space-time manifold h as where T 0 is the initial guess vector, Q(t) is a function of time t and satisfies that Q(t) > 0, Q(0) = 1, and it is a monotonically increasing function of t with Q(∞) = ∞. If the trajectory of T can be kept on the manifold, as the time t increases T approaches to the solution of Equation (13). In order to keep the trajectory of T on the manifold, the consistency equation is required to be satisfied: in which ∇ T denotes the gradient operator with respect to vector T. Since Equation (13) is a scalar equation, it is not possible to determine the evolution of vector T (that is dT dt ) uniquely. Assuming the evolution of vector T is in the direction of u (i.e., u is the direction of evolution for vector T), we have: where λ is a proportional constant. Skipping the lengthy derivation, the following iteration equation of T can be obtained: where r is the relaxation parameter which is designed for making the convergence of T more stable, and v = Bu.
In the process of derivation, a parameter, a 0 , is defined as: where θ is the angle between the residual vector F and v. It is proved that in order to ensure the trajectory of T on the manifold, the value of a 0 needs to satisfy 1 ≤ a 0 ≤ 4. In addition, the optimal relaxation parameter is r = 1 − a 0 2 . The selection of u becomes the important key to guarantee the trajectory of T on the manifold. Theoretically speaking, the best selection of u makes Bu − F = 0 (or u = B −1 F if the inverse of the Jacobian matrix exists). In such a manner, Equation (15) reduces to the well-known Newton-Raphson method. For an ill-posed system, if it is difficult to find the inverse of the Jacobian matrix then other direction is needed.
For the standard BVP of the nonlinear heat equilibrium problem, it is expected that the mathematical model is a well-posed one. However, if the discretization method makes the Jacobian matrix a dense one and the truncation error in the mathematical operation may make the numerical process more unstable. The local MQ RBF collocation method makes the Jacobian matrix a sparse one, which evidently can help reduce the risk of encountering an unstable process. In addition, while a large-scale problem is solved, the localized method, such as the local MQ RBF collocation method, will have advantage over the global method due to the memory limitation of the machine. In such a case, to seek the solution of Bu − F = 0 is much easier than seeking the inverse of the Jacobian matrix. In our programming with MATLAB, the solution of Bu − F = 0 is determined by the left division operator (\), which basically is based on lower-upper (LU) decomposition method. The left division operator in MATLAB can support the operations on a sparse matrix.
The iteration will be terminated when one of the following two conditions has been achieved: (1) the root mean square error (RMSE) for Equation (13) is less than a prescribed value, ε, i.e., RMSE < ε.
(2) the number of iteration steps exceeds the maximum number I max .

Example 1
In this example, the interested domain is a circle with a radius of 1 and the center of the circle is located at (x,y) = (5,5). The temperature dependent heat conduction coefficient is K(T) = T. The designated solution is T(x,y) = 2xy. The shape parameter c = 1.5. On the whole boundary, the case of using Dirichlet boundary conditions is first tested. Totally, 180 boundary points and 705 uniformly distributed interior points are prescribed. The initial guess of temperature distribution is a constant distribution with T = 7.
The convergence criterion is set as RMSE ≤ 10 −8 and the maximum iteration number is I max = 100. The points inside the supporting neighborhood (m p ) are 30, 40, and 50 points. It can be seen from Table 1 that for different values of m p , the RMSE reaches the convergence criterion within 28 steps. However, the maximum absolute relative error percentage decreases as m p increases. It should be remarked here that the computation time increases as m p increases since the computation of the inverse of H matrix shown in Equation (9) requires more computation time while m p is larger. The changes in RMSE versus the number of iteration steps are illustrated in Figure 1. As shown in this figure, one can find that the RMSE continues decreasing as the number of iteration steps increases. We wonder how good the current approach can achieve while the convergence criterion is not prescribed. We then fix the value of m p = 30 and the maximum iteration number is 100, and check the RMSE versus the number of iteration step as shown in Figure 2. It then can be seen while RMSE reaches about 10 −12 , RMSE begins to fluctuate as the number of iteration step increases. The minimum RMSE within 100 steps is 6.7493 × 10 −12 and the maximum absolute relative error percentage is 0.2706% which is the same as we set RMSE ≤ 10 −8 as the convergence criterion.
uniformly distributed interior points are prescribed. The initial guess of temperature distribution is a constant distribution with T = 7.
The convergence criterion is set as RMSE 8 10 − ≤ and the maximum iteration number is Imax = 100. The points inside the supporting neighborhood (mp) are 30, 40, and 50 points. It can be seen from Table 1 that for different values of mp, the RMSE reaches the convergence criterion within 28 steps. However, the maximum absolute relative error percentage decreases as mp increases. It should be remarked here that the computation time increases as mp increases since the computation of the inverse of H matrix shown in Equation (9) requires more computation time while mp is larger. The changes in RMSE versus the number of iteration steps are illustrated in Figure 1. As shown in this figure, one can find that the RMSE continues decreasing as the number of iteration steps increases. We wonder how good the current approach can achieve while the convergence criterion is not prescribed. We then fix the value of mp = 30 and the maximum iteration number is 100, and check the RMSE versus the number of iteration step as shown in Figure 2. It then can be seen while RMSE reaches about 10 −12 , RMSE begins to fluctuate as the number of iteration step increases. The minimum RMSE within 100 steps is 6.7493 12 10 − × and the maximum absolute relative error percentage is 0.2706% which is the same as we set RMSE 8 10 − ≤ as the convergence criterion.  Now, the mixed type boundary condition is considered. For the upper circular boundary, the Dirichlet boundary conditions are prescribed. Furthermore, for the lower circular boundary, the Neumann boundary conditions are prescribed. Again, the points inside the supporting neighborhood (mp) are 30, 40, and 50 points. The convergence criterion is set as RMSE 8 10 − ≤ and the maximum iteration number is Imax = 100. It can be seen from Table 2 that for different values of mp, Now, the mixed type boundary condition is considered. For the upper circular boundary, the Dirichlet boundary conditions are prescribed. Furthermore, for the lower circular boundary, the Neumann boundary conditions are prescribed. Again, the points inside the supporting neighborhood (m p ) are 30, 40, and 50 points. The convergence criterion is set as RMSE ≤ 10 −8 and the maximum iteration number is I max = 100. It can be seen from Table 2 that for different values of m p , the RMSE reaches the convergence criterion within 29 steps. The maximum absolute error percentage reaches the minimum value for m p = 40. In addition, from Tables 1 and 2, one can find out that the computation result for Dirichlet type boundary value problem is more precise than that for the mixed type boundary value problem. The contour plots for the physical quantities inside the domain for various m p are depicted in Figure 3 and the comparison between numerical solutions and exact solution shows that the current approach can obtain the satisfactory solution.     Before we close this example, the influence of the shape parameter c is investigated using the Dirichlet boundary conditions with the convergence criterion is set as RMSE < 10 −12 . The value of m p is fixed at 30, and we examine the minimum RMSE before the I max = 100 is reached and the maximum absolute relative error percentage versus the shape parameter (c begins from 0.5 to 10, and the difference between two adjacent c values is 0.1). The results are illustrated in Figure 4. It can be found that the shape parameter does not influence RMSE, apparently; however, the maximum absolute relative error percentage increases as the shape parameter increases. When c is less than 2, the maximum absolute relative error percentage is less than 1%. In the following examples, the value of the shape parameter is fixed as c = 1.5.

Example 2
In this example, an annular region is considered. The outer radius for this annular region is 3 and the inner radius is 1.5. The center of this annular region is (0,0). The temperature dependent conductivity k(T) = T. The exact solution is designed as T(x,y) =  Table 3. It can be seen that the maximum absolute relative error percentage for case 1 is the lowest one and that for case 3 is the highest one. Both cases 2 and 3 have part of the boundary with the Dirichlet boundary condition and Neumann boundary condition on the remaining boundary. It seems that when the Neumann boundary condition is prescribed on the inner circle (case 3) the maximum absolute relative error percentage is higher than the case when the Neumann boundary condition is prescribed on the outer circle (case 2).

Example 2
In this example, an annular region is considered. The outer radius for this annular region is 3 and the inner radius is 1.5. The center of this annular region is (0,0). The temperature dependent conductivity k(T) = T. The exact solution is designed as T(x,y) = 2 ln(x 2 + y 2 ). Totally, 180 points are arranged on the boundary, 90 on the inner circle and 90 on the outer circle. Furthermore, 1456 uniformly distributed interior points are selected. The shape parameter c = 1.5 and the convergence criterion is set as RMSE ≤ 10 −16 . The maximum number of iteration steps is 100. The value of m p is fixed at 30. The initial guess of temperature distribution is a constant distribution with T = 7.
Three cases are performed: (1) case 1: on whole boundary, the Dirichlet boundary conditions are prescribed; (2) case 2: on the inner circle, the Dirichelt boundary conditions are given and on the outer circle the Neumann boundary conditions are given; (3) case 3: on the inner circle, the Neumann boundary conditions are given and on the outer circle the Dirichlet boundary conditions are given. The computation performances for these three cases are tabulated in Table 3. It can be seen that the maximum absolute relative error percentage for case 1 is the lowest one and that for case 3 is the highest one. Both cases 2 and 3 have part of the boundary with the Dirichlet boundary condition and Neumann boundary condition on the remaining boundary. It seems that when the Neumann boundary condition is prescribed on the inner circle (case 3) the maximum absolute relative error percentage is higher than the case when the Neumann boundary condition is prescribed on the outer circle (case 2).

Example 3
In this example, a square region defined by (x, y) 0 ≤ x ≤ 1; 0 ≤ y ≤ 1 is considered. The temperature-dependent thermal conductivity is k(T) = exp(T). In the previous two examples, the temperature-dependent thermal conductivity is linear with respect to the temperature. In this example, the temperature-dependent thermal conductivity now is nonlinear with respect to the temperature. The exact solution is designed as: T(x, y) = log(xy + 5). The shape parameter c = 1.5 and the convergence criterion is set as RMSE ≤ 10 −16 . The maximum number of iteration steps is I max = 100. The value of m p is first fixed at 30. The initial guess of temperature distribution is a constant distribution with T = 1.7. Totally, 396 uniformly distributed boundary points and 9604 interior points are used. The Dirichlet boundary conditions are prescribed on the whole boundary.
The result is illustrated in Figure 5. The minimum RMSE is 4.3416 × 10 −11 and the maximum absolute relative error percentage is 3.5798%. The reason for larger absolute relative error percentage comes from the fact that the nonlinearity for the temperature-dependent thermal conductivity is higher for this case and in addition more points are considered. From Figure 5, one can say that the current approach using m p = 30 does not yield reasonable results. When the value of m p is changed to 50, the minimum RMSE increases to 4.7730 × 10 −11 ; however, the maximum absolute relative error percentage reduces to 0.5232%. It implies that when the nonlinearity increases, increasing m p may help us to obtain a more accurate result.

Example 3
In this example, a square region defined by ( ) In the previous two examples, the temperature-dependent thermal conductivity is linear with respect to the temperature. In this example, the temperature-dependent thermal conductivity now is nonlinear with respect to the temperature. The exact solution is designed as: . The shape parameter c = 1.5 and the convergence criterion is set as RMSE 16 10 − ≤ . The maximum number of iteration steps is Imax = 100. The value of mp is first fixed at 30. The initial guess of temperature distribution is a constant distribution with T = 1.7. Totally, 396 uniformly distributed boundary points and 9604 interior points are used. The Dirichlet boundary conditions are prescribed on the whole boundary.
The result is illustrated in Figure 5. The minimum RMSE is 4.3416 × 10 −11 and the maximum absolute relative error percentage is 3.5798%. The reason for larger absolute relative error percentage comes from the fact that the nonlinearity for the temperature-dependent thermal conductivity is higher for this case and in addition more points are considered. From Figure 5, one can say that the current approach using mp = 30 does not yield reasonable results. When the value of mp is changed to 50, the minimum RMSE increases to 4.7730 × 10 −11 ; however, the maximum absolute relative error percentage reduces to 0.5232%. It implies that when the nonlinearity increases, increasing mp may help us to obtain a more accurate result.

Example 4
In this example, we redo the problem in example 1. The only difference is that for the lower half circle the nonlinear type boundary condition as Equation (4) is provided. All setup and parameters are the same as that in example 1. The parameters used in equation (4) are: h = 10,σ = 5.667 × 10 −8 , T f = 1, T s = 1 and ε = 1.
The contour plots for numerical results are depicted in Figure 6 The minimum RMSE is 8.5895 × 10 −12 and the maximum absolute relative error percentage is 1.0617%. It can be seen that the maximum relative error percentage of the result in this example is higher than that using Dirichlet and Neumann conditions shown in Table 2. This reveals that when the Neumann boundary condition is replaced by a nonlinear one like Equation (4), the accuracy of solution decreases. When we change m p to 50, the minimum RMSE reduces to 7.7318 × 10 −12 and the maximum absolute relative error percentage reduces to 0.1180%. Once again, increasing m p may help us to obtain a more accurate result.

Example 4
In this example, we redo the problem in example 1. The only difference is that for the lower half circle the nonlinear type boundary condition as Equation (4) is provided. All setup and parameters are the same as that in example 1. The parameters used in equation (4)  The contour plots for numerical results are depicted in Figure 6 The minimum RMSE is 8.5895 × 10 −12 and the maximum absolute relative error percentage is 1.0617%. It can be seen that the maximum relative error percentage of the result in this example is higher than that using Dirichlet and Neumann conditions shown in Table 2. This reveals that when the Neumann boundary condition is replaced by a nonlinear one like Equation (4), the accuracy of solution decreases. When we change mp to 50, the minimum RMSE reduces to 7.7318 × 10 −12 and the maximum absolute relative error percentage reduces to 0.1180%. Once again, increasing mp may help us to obtain a more accurate result.

Example 5
In this example, a peanut shape domain is considered. The boundary curve of this peanut shape is described as: where r, θ are the cylindrical coordinate components. The value of m p is 50. The maximum number of iteration steps is fixed at 100, and the convergence criterion is set as RMSE ≤ 10 −16 .
The thermal conductivity k(T) = 1/T. The designated solution is T(x,y) = exp(xy + 2). The initial guess of temperature is set as T(x,y) = 7. Totally, 90 uniformly distributed boundary points and 2508 interior points are used. The Dirichlet boundary conditions are prescribed for all boundary points first. The minimum RMSE is reported as 7.548 × 10 −13 and the maximum absolute relative error percentage is 0.0697%. Next, the mixed type boundary condition with Dirichlet boundary conditions given for boundary points having y ≥ 0 and Neumann boundary conditions given for boundary points having y < 0 is considered. The minimum RMSE for this case is 7.9948967 × 10 −13 and the maximum absolute relative error percentage is 0.2187%. It can be found that while boundary condition belongs to mixed type, the accuracy of numerical solution decreases which matches results from previous examples and many previous research studies.

Example 6
In this example, we consider a boundary value problem without the exact solution. The interested domain is a circle with a radius of 1 and the center of the circle is located at (x,y) = (5,5). The temperature dependent heat conduction coefficient is K(T) = T. The shape parameter c = 1.5. On the whole boundary, the Dirichlet boundary conditions as T(x, y) = xy are given. Totally, 180 boundary points and 705 uniformly distributed interior points are adopted. The initial guess of temperature distribution is a constant distribution with T = 7. The value for m P is 50.
Using the current approach, the temperature distribution at all selected points (180 boundary points and 705 interior points) can be found. However, the exact solution is unknown now and the only way to check the performance of the current approach is to examine the coincidences of boundary values for boundary points. In order to perform that, one extra collocation point between two adjacent boundary points is added. After the temperature distribution for all points are obtained, then we calculate the temperature values for these extra boundary collocation points by selecting 50 nearest supporting points. The computed values and the designated boundary values are compared as shown in Figure 7. For this case, the minimum RMSE is 9.7741 × 10 −11 and the maximum relative error percentage for all extra boundary collocation points is 5.3472 × 10 −4 %. It can be seen that the current approach method still can obtain reasonable results while the exact solution cannot be found. The contour plot of temperature distribution is depicted in Figure 8.

Remarks
Before we close this section, some remarks are made here. When the proposed approach is used for calculating the thermal equilibrium problem in which the thermal conductivity coefficient is not expressed as an analytical expression but is given as a numerical array, the following steps can be adopted. One can first find the mathematical representation for the thermal conductivity coefficient simply by the curve fitting technique, then the required representation of the thermal conductivity coefficient, the derivative of thermal conductivity coefficient with respect to temperature and the second order derivative of thermal conductivity coefficient with respect to temperature, which are necessary for the current approach, can be obtained. On the other hand, one also can obtain the thermal conductivity coefficient from the numerical array, and the derivatives by numerical differentiation.
In this article, only the energy balance is considered. However, in the real world the problem may need to consider the balance equations for energy, mass, and momentum. It means we might have nonlinear coupled equations to be solved simultaneously. When such a case is encountered, the current approach has no difficulty to solve. One only needs to express all unknown physical fields by the localized radial basis function and substitute the expressions into the governing equation and boundary condition and construct a system of nonlinear algebraic equations. The nonlinear solver used in the current approach can then find the numerical solution from the initial guess.

Remarks
Before we close this section, some remarks are made here. When the proposed approach is used for calculating the thermal equilibrium problem in which the thermal conductivity coefficient is not expressed as an analytical expression but is given as a numerical array, the following steps can be adopted. One can first find the mathematical representation for the thermal conductivity coefficient simply by the curve fitting technique, then the required representation of the thermal conductivity coefficient, the derivative of thermal conductivity coefficient with respect to temperature and the

Remarks
Before we close this section, some remarks are made here. When the proposed approach is used for calculating the thermal equilibrium problem in which the thermal conductivity coefficient is not expressed as an analytical expression but is given as a numerical array, the following steps can be adopted. One can first find the mathematical representation for the thermal conductivity coefficient simply by the curve fitting technique, then the required representation of the thermal conductivity coefficient, the derivative of thermal conductivity coefficient with respect to temperature and the

Conclusions
In this paper, the standard boundary value problems of the nonlinear heat equilibrium equation are investigated. The roles of collocation points and central points are clearly defined such that seeking values for extra collocation points becomes possible. Six numerical examples are provided to show the validity of the current approach. It is found that when the nonlinearity of the problem becomes more apparent the solution become less accurate. Increasing the number of supporting points (m p ) can help to obtain a more accurate result. In all cases, m p = 50 can achieve good results. The boundary value problem without the exact solution is also examined; it is found that the boundary values for extra boundary collocation points can match the designated values very well.