Numerical simulations of 1 D inverse heat conduction problems using overdetermined RBF-MLPG method

This paper proposes a numerical method to deal with the one-dimensional inverse heat conduction problem (IHCP). The initial temperature, a condition on an accessible part of the boundary and an additional temperature measurements in time at an arbitrary location in the domain are known, and it is required to determine the temperature and the heat flux on the remaining part of the boundary. Due to the missing boundary condition, the solution of this problem does not depend continuously on the data and therefore its numerical solution requires special care especially when noise is present in the measured data. In the proposed method, the time variable is eliminated by using finite differences approximation. The method uses a weak formulation of the problem to enjoy the stability condition. To avoid the numerical integration on the whole domain, the weak form equations are constructed on local subdomains. The approximate solution is assumed to be a linear combination of Multi Quadric (MQ) radial basis function (RBF) constructed on nodal points in the domain and on the boundary. Since the problem is known to be ill-posed, Thikhonov regularization strategy is employed to solve effectively the discrete ill-posed resultant linear system.


Introduction
Development of constructive methods for the numerical solution of mathematical problems is a main branch of mathematics.For the numerical solutions of PDEs, finite differences, finite elements, finite volume and spectral methods are traditional well-developed main numerical methods.In the both mathematics and engineering community, meshless methods have attracted much attention in recent years.Extensive developments have been made in several varieties of meshless methods and applied to many applications in science and engineering.These methods exist under different names, such as: the diffuse element method (DEM) [1], the hp-cloud method [2], Meshless Local Petrov-Galerkin (MLPG) method [3,4,5,6,7,8,9,10,11,12], the meshless local boundary integral equation (LBIE) method [13,14], the partition of unity method (PUM) [15], the meshless collocation method based on radial basis functions (RBFs) [16], the smooth particle hydrodynamics (SPH) [17], the reproducing kernel particle method (RKPM) [18], the radial point interpolation method [19] and so on.Inverse problems arise in many heat transfer situations when experimental difficulties are encountered in measuring or producing the appropriate boundary conditions.It is sometimes necessary to determine, the surface temperature of a body from a measured temperature history at a given fixed location inside the body, when the surface itself is inaccessible for measurements.Due to the missing boundary conditions and therefore maximum principle, the solution of the under consideration problem does not depend continuously on the data, and therefore it is known as an ill-posed problem in the sense as described by Hadamard [20] and therefore, its numerical solution requires special care.This ill-posedness is motivated by the fact that in general, in applications the data will be measured quantities and therefore always contaminated by errors.In this paper, we are interested on the development of a new approach called RBF-MLPG method, for numerical solutions of inverse heat conduction problems.We proposed this method in [21] for solving second order elliptic PDEs and proved that the exponential convergence rate found in MQ interpolation can be transformed into overdetermined RBF-MLPG method.The mathematical formulation of the IHCP considered in this study can be described by the following boundary value problem.The governing heat conduction equation in a slab geometry, namely has to be solved subject to the initial and boundary condition where u is the temperature, ∂ u(x,t) ∂ n is the heat flux, n is the outward normal at the boundary of the slab and u 0 and g are prescribed functions.The thermal diffusivity is assumed constant and, for simplicity, taken to be unity.In addition, temperature readings are provided at an arbitrary space location where h is a known function.When d = 0 the equations (1.1)-(1.4)reduce to the direct problem.Based on the formulation in equations (1.1)-(1.4) it is required to determine the temperature and the heat flux See [22,23,24,25] for some papers on inverse heat conduction problems.Since the problem is not well-posed, to overcome the instability problem, one possible approach is to record the temperature history at more locations than the number of unknowns in the algebraic equations to be solved.Then the system of equations becomes overdetermined and can be solved using a least squares method.It is shown in [21] that in the case of overdetermined RBF-MLPG, since more data are known from the problem, the solution accuracy increases.Moreover, a regularization technique which modifies the least squares method by adding a regularization coefficient, as described by Tikhonov and Arsenin [26] is proposed to solve the finial Overdetermined linear system of equations.Two test problems are investigated and the results are compared with the analytical solutions.The effects of the location of internal temperature measurements in relation to the distance from the surface where the temperature and the heat flux are unknown are also investigated.
For exact data the results for the surface temperature and heat flux obtained using the numerical methods show good agreement with the analytical solution.However, when noise is present in the data, instabilities can be alleviated by using Tikhonov regularization.

Finite differences approximation
Suppose the time interval [0, T ] is discretized uniformly into K subintervals; define t k = k△t, k = 0, 1, ..., K, where △t = T /K is the time step.Let u k = u k (x) := u(x,t k ) be the exact solution u restricted to time t k .Then, the finitedifference approximation of the time derivatives in the θ method is given as follows Considering (1.1) at the time instants k∆t and (k + 1)∆t, one obtains, respectively Hence and from (2.1), we have Using Crank-Nicholson scheme, θ = 1 2 , (2.2) becomes: So, assuming the field u k being known from the computation in the previous time step, the field variable u k+1 can be obtained via elliptic type PDE (2.3).

Approximation of the solution
For the RBF-MLPG formulation of the problem (1.1)-(1.4)we consider a set of N scattered RBF centers and for any radial basis kernel Φ, the numerical solution at time instant k∆t is expanded in the form of where λ k := (λ k 1 , ..., λ k N ) T are the unknown coefficients to be determined.In this paper we consider the multiquadric(MQ) kernel where r ∈ R and c > 0 is a constant.So, the unknown solution u k+1 in (2.3) is approximated by a linear combination of MQs in the form of (2.4)

Local weak equations
To find the unknown coefficient in the approximation (2.4), it should be tested to satisfy in the governing equations (1.1-1.4).Our testing method in this paper contains a combination of strong and local weak testing.So, we consider a set of M test nodes in the domain and on the boundary of the problem and: 1-For x i = 1 using Boundary condition (1.3)we have: 2-For x i = d using Condition (1.4)we have: 3-For x i in the interior of domain, we consider a subdomain Ω s i around x i and the weak form of Eq.( 2.3) for some test function v on subdomain Ω s i will be written as follows: ∫ Using integration by parts, one has: and using the test function: x / ∈ Ω i s , the following local weak equation will be obtained: and considering the approximation (2.4) we have: Combining equations (2.5),(2.6)and (2.9) finally we obtain an overdetermined system of M equations with N unknowns λ k+1 .Note that λ 0 will be obtained from initial condition.

Evaluation of integrals
For evaluating ∫ √ c 2 + x 2 dx, by using the triangular change of the variable )) International Scientific Publications and Consulting Services and, substituting Z j i = x i − x j we have: Due to small size of the radius of the subdomains, r, the right hand domain integral is approximated as follows: ∫

Regularization
After discretization, the ill conditioning of the original problem (1.1)-(1.4)do not vanish and it will be moved to the final system of linear equations (2.5)-(2.9).Several regularization methods have been developed for solving this kind of ill-conditioning problems [27].In this paper we use the Tikhonov regularization [26] to solve the final linear system (2.5)-(2.9)which we name this system as Aλ = b.The Tikhonov regularized solution of the final linear system, say λ α * , is defined to be the solution to the following least squares problem: where ∥ • ∥ denotes the usual Euclidean norm and α is the regularization parameter.The interpretation of Tikhonov regularization indicates that it keeps the residual ∥ Aλ − b ∥ 2 small and stablizes by preventing λ α from becoming large through the penalty term α 2 ∥ λ ∥ 2 .we use the generalized cross-validation (GCV) method [28] to find the regularization parameter α.The GCV method chooses the optimal regularization parameters by minimizing the following GCV function: 2  (2.11) where K = (A T A + α 2 I) −1 A T is a matrix which produces the regularized solution λ α = K b.In our computation, we use the Matlab code developed by Hansen [29] for solving the discrete ill-conditioned system of equations (2.5)-(2.9).

Numerical demonstrations
In this section, first we are interested on numerically analyzing the convergence of the proposed method and its sensitivity with respect to the parameters by using one test problem, where no noise is inserted in the data.Then,the noise effect on the approximate numerical solution is studied by using two test problems.The root mean square (RMS) error of the approximate solution presented in this paper is defined as: where Z is the set of all evaluation points and z i ∈ Z.Although from Fig. 1 more number of test nodes give better approximation, but the increase is not remarkable and we use 2n × n matrix in the rest of this paper.Fig. 2 presents the RMS error versus the shape parameter c obtained with m = 2n which suggests a relatively large domain for choosing c,the interval (0.5, 1.8).It also should be noted that in MQ RBF, ill conditioning happens as c increases.So, hereafter, we fix the shape parameter at c = 1.Fig. 3 shows the numerical and exact solution at x = 0 while Fig. 4 is presented to show the numerical approximation for the first derivative at x = 0 and its exact value versus time variable t.To see that how the numerical solutions depend on the location of internal temperature measurement, we present Fig. 5 and Fig. 6.Comparing these figures, it is observed that as d approaches to the right boundary, x = 1, the error increases.This increase is more significant for d > 0.5 as shown in Fig. 6.All parameters in this subsection is chosen as∆t = 0.1, h = 0.05 (the distance between RBF centers), n = 21,and d = 0.4.

Example 3.2. measured data
In this test problem, we consider the case of error measured data which can be simulated by adding random noisy data , ε i , to the exact data at each node on x = d and x = b, as follows: 1), l = 1, 2. Two test problems with the data taken from the exact solutions u(x,t) = x 3 + t 2 and u(x,t) = 1 4 e t cos(x) with different noise level is considered and the results are presented in Fig. 7 and Fig. 8.The parameters are chosen as follows: c = 1, ∆t = 0.2, d = 0.4, h = 0.05, h 0 = 0.025, m = 2n, r = 0.5h 0 , where h is the distance between RBF centers and h 0 is the distance between collocation points, r is the radius of subdomains, m is the number of equations and n is the number of unknowns.A numerical approach was proposed to effectively solve the one-dimensional inverse heat conduction problem.The major problem was determination of temperature on an inaccessible part of the boundary.The time variable was discretized by using finite differences approximation.Then, weak form of the discretized equations were constructed on local subdomains.MQ RBF approximation was used to discretize the spatial variations of the field variable.Tikhonov regularization was used to overcome the problem of ill-conditioning.Test problems verified the accuracy of the proposed approach.

Example 3 . 1 . 1 .
Convergence and stability analysis Consider Problem (1.1)-(1.4)with f (x,t) = 0 and other data chosen in such a way that the exact solution is u(x,t) = exp(x + t).The RMS error of the numerical solution at time instant t = 1, versus the number of test nodes, m, is presented in Fig. From this figure, the error with square 21 × 21 matrix is 5.4 × 10 −4 and for the 70 × 21 matrix is 4.2 × 10 −4 .

Figure 1 :Figure 2 :
Figure 1: RMS error versus number of test nodes

Figure 3 :
Figure 3: Numerical and exact solution at x = 0 versus time

Figure 4 :Figure 5 :Figure 6 :
Figure 4: Numerical approximation for the first derivative at x = 0 and its exact value versus time

Figure 7 : 2 0
Figure 7: The exact and numerical solution obtained with different noise level for test problem 2 with the exact solution u(x,t) = x 3 + t 2

Figure 8 :
Figure 8: The exact and numerical solution obtained with different noise level for test problem 2 with the exact solution u(x,t) = 1 4 e t cos(x)