A Meshless Method to Determine a Source Term in Heat Equation with Radial Basis Functions

Abstract

This paper considers a numerical method based on the radial basis functions for the inverse problem of heat equation; the inverse problem is determining an unknown source term subject to the overdetermination along with the usual initial boundary conditions, and the unknown source term is only time-dependent. The radial basis functions method is a meshless method with high accuracy for the inverse problem. Some numerical experiments using this method are presented and discussed.

1. Introduction

The inverse problem of parabolic equations appears naturally in a wide variety of physical and engineering settings, such as elasticity, hydrology, material sciences, heat transfer, medical imaging, transport problems, and control theory.

In this paper, we consider an inverse problem for determining an unknown source term in a heat equation. There are many numerical methods which have been developed to solve the inverse problem of parabolic equations. One approach is the method of output least squares, which assumes that the unknown source function is a specific functional form depending on some parameters and then seeks to determine optimal parameter values which minimized an error function based on the overdetermined data. Other effective and broadly applicable techniques such as finite difference method, finite elements method, boundary elements method, and finite volume method are widely used in inverse problems; however, the finite elements method approximates the solution by using low-order piecewise polynomials, and in the finite difference method, the derivatives of the solution are approximated by difference quotients. These methods depend on a suitable generation of meshes, which is difficult for problems with very complicated and irregular geometries. The method we used in this paper is the meshless method based on the radial basis functions.

Meshless methods are a class of numerical methods for solving partial differential equations. In these methods, mesh generation on the spatial domain of the problem is not needed, and this property is the main advantage of these techniques over the mesh-dependent methods such as finite difference methods and finite element methods. In addition, the radial basis functions method is a technique for interpolation on scattered data; this method uses the distributed nodal points to approximate the unknown function; the distribution of nodes could be selected regularly or randomly in the analyzed domain, and there is no demand for the geometry of the domain. The meshless method based on the radial basis functions provides a global interpolation formula not only for the solution but also for the derivatives of the solution and possesses a high order of accuracy. A number of investigators have considered the meshless method based on the radial basis functions: Buhmann [1], Wendland [2, 3], Dehghan and Tatari [4, 5], Wu et al. [6, 7], and so forth.

In Section 2, we will give an outline of the radial basis functions. In Section 3, we solve the inverse problem using the meshless method based on the radial basis functions. Numerical experiments will be given in Section 4 to illustrate the application of the method.

2. An Outline of the Radial Basis Functions (RBF)

The radial basis functions method is a meshless method, it has been introduced for interpolation of scattered data, and it is very suitable for the numerical solution of partial differential equations.

Let 𝑟 be the distance between a fixed point 𝑥^∗∈𝑅^𝑑 and any 𝑥∈𝑅^𝑑; that is, 𝑟=‖𝑥−𝑥^∗‖₂.

Define 𝜙∶𝑅⁺→𝑅, and 𝜙(‖𝑥−𝑥^∗‖₂) is a univariate function.

Let 𝑥₁,𝑥₂,…𝑥_𝑁 be a given set of distinct points. According to the definition of the function 𝜙, we denote 𝜙_𝑗(𝑥)=𝜙(‖𝑥−𝑥_𝑗‖₂). For a given data {𝑥_𝑗,𝑢_𝑗},𝑗=1,2,…,𝑁, the interpolating function 𝑢(𝑥) can be written as follows: 𝑢(𝑥)=^𝑁_𝑗=1𝜆_𝑗𝜙𝑥−𝑥_𝑗2,𝑥∈𝑅^𝑑, (1) where 𝜆_𝑗 is the unknown RBF coefficient and 𝑁 is the number of data points.

In order to compute 𝜆_𝑗, assume that we want to interpolate the values 𝑢(𝑥_𝑗); that is, 𝑢𝑥_𝑗=𝑢_𝑗. (2)

By substituting 𝑥_𝑖(𝑖=1,2,…,𝑁) for 𝑥 in (1), we obtain the following linear system. ΦΛ=𝑈, (3) where 𝜙Φ=¹𝑥₁𝜙₂𝑥₁⋯𝜙_𝑁𝑥₁𝜙₁𝑥₂𝜙₂𝑥₂⋯𝜙_𝑁𝑥₂𝜙⋮⋮⋮₁𝑥_𝑁𝜙₂𝑥_𝑁⋯𝜙_𝑁𝑥_𝑁,𝜆Λ=¹𝜆₂⋮𝜆_𝑁,𝑢𝑈=¹𝑢₂⋮𝑢_𝑁, (4) so we can get the unknown coefficient Λ from the linear system ΦΛ=𝑈. The matrix Φ has been shown to be positive definite and is invertible for distinct set of scattered set of points. But in most cases, the matrix Φ is very ill-conditioned, so we usually use the more stable method instead of the direct method to solve the coefficient Λ.

The well-known radial basis functions are listed in Table 1.

Table 1 

Some well-known functions that generate RBFs.

3. The Meshless Method for the Inverse Problem

In this section, we consider a class of inverse problems as follows.

Find 𝑢=𝑢(𝑥,𝑡) and 𝑓(𝑡) which satisfy 𝑢_𝑡(𝑥,𝑡)=𝑢_𝑥𝑥()(𝑡),()𝑥,𝑡+𝑓𝑥,𝑡∈𝑄_𝑇=()×(),𝑢()(𝑥)(),𝑢()0,𝑙0,𝑇𝑥,0=𝜑,𝑥∈0,𝑙0,𝑡=ℎ₀(𝑡),𝑢(𝑙,𝑡)=ℎ_𝑙(𝑡),𝑡∈(0,𝑇), (5) and the additional specification 𝑢𝑥^∗,𝑡=𝐸(𝑡),𝑥^∗∈(0,𝑙),𝑡∈(0,𝑇), (6) where 𝜑(𝑥), ℎ₀(𝑡), ℎ_𝑙(𝑡), and 𝐸(𝑡)≠0 are known functions and 𝑥^∗ is a fixed prescribed point in the admissible set.

If the functions 𝜑(𝑥), ℎ₀(𝑡), ℎ_𝑙(𝑡) and 𝐸(𝑡) satisfy consistence, then the inverse problem (5)-(6) has a unique solution [8].

This kind of inverse problems of identifying unknown source control function has been studied by several authors [8, 9]. The approaches used in these papers are the finite difference method based on the trace type functional formulation and the boundary element method. In this paper, we use the meshless method based on radial basis functions combined with the trace type functional formulation to solve the inverse problem as follows.

If the functions pair (𝑢,𝑓) solves the inverse problem (5)-(6), applying the trace type functional formulation, we have 𝐸^′(𝑡)=𝑢_{𝑥𝑥𝑥=𝑥∗}+𝑓(𝑡). (7) From this, we get 𝑓(𝑡)=𝐸^′(𝑡)−𝑢_{𝑥𝑥𝑥=𝑥∗}. (8)

Substituting (8) into (5), we have the following initial boundary value problem: 𝑢_𝑡=𝑢_𝑥𝑥+𝐸^′(𝑡)−𝑢_𝑥𝑥|_𝑥=𝑥∗,(𝑥,𝑡)∈𝑄_𝑇,𝑢()(𝑥)(),𝑢()𝑥,0=𝜑,𝑥∈0,𝑙0,𝑡=ℎ₀(𝑡),𝑢(𝑙,𝑡)=ℎ_𝑙(𝑡),𝑡∈(0,𝑇). (9)

The approximate function 𝑢(𝑥,𝑡) of 𝑢(𝑥,𝑡) can be represented as 𝑢(𝑥,𝑡)=^𝑁_𝑗=1𝜆_𝑗(𝑡)𝜙_𝑗(𝑥), (10) where 𝜙_𝑗(𝑥) is the radial basis functions described in Section 2.

Substituting (10) into (9), the system (9) can be rewritten as ^𝑁_𝑗=1𝜆^′_𝑗(𝑡)𝜙_𝑗(𝑥)=^𝑁_𝑗=1𝜆_𝑗(𝑡)𝜙_𝑗^′′(𝑥)+𝐸^′(𝑡)−^𝑁_𝑗=1𝜆_𝑗(𝑡)𝜙_𝑗^′′𝑥^∗,(𝑥,𝑡)∈𝑄_𝑇,_𝑁𝑗=1𝜆_𝑗(0)𝜙_𝑗(𝑥)=𝜑(𝑥),𝑥∈(0,𝑙_𝑁𝑗=1),𝜆_𝑗(𝑡)𝜙_𝑗(0)=ℎ₀(𝑡),_𝑁𝑗=1𝜆_𝑗(𝑡)𝜙_𝑗(𝑙)=ℎ_𝑙(𝑡),𝑡∈(0,𝑇). (11)

For the solution of the system (11), we apply a simple one-step forward difference formula to the time and obtain ^𝑁_𝑗=1𝜆_𝑗𝑡_𝑚+1−𝜆_𝑗𝑡_𝑚𝜙Δ𝑡^𝑗(𝑥)=^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗^′′(𝑥)+𝐸𝑡_𝑚+1−𝐸𝑡_𝑚−Δ𝑡^𝑁_𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗^′′𝑥^∗,(𝑥,𝑡)∈𝑄_𝑇,_𝑁𝑗=1𝜆_𝑗(0)𝜙_𝑗(𝑥)=𝜑(𝑥),𝑥∈(0,𝑙_𝑁𝑗=1),𝜆_𝑗𝑡_𝑚𝜙_𝑗(0)=ℎ₀𝑡_𝑚,_𝑁𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗(𝑙)=ℎ_𝑙𝑡_𝑚,𝑡∈(0,𝑇). (12)

From the iterative method, we can obtain 𝜆_𝑗(𝑡_𝑚), 𝑗=1,2,…,𝑁, 𝑚=1,2,…,𝑀, and combining (10) and (8), we get the numerical solution as follows: 𝑢𝑥,𝑡_𝑚+1=^𝑁_𝑗=1𝜆_𝑗𝑡_𝑚+1𝜙_𝑗(𝑥)=^𝑁_𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗(𝑥)+Δ𝑡^𝑁_𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗^′′(𝑥)+𝐸𝑡_𝑚+1−𝐸𝑡_𝑚−Δ𝑡^𝑁_𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗^′′𝑥^∗ (13) and the approximate solution as follows: 𝑓𝑡_𝑚+1=𝐸𝑡_𝑚+1−𝐸𝑡_𝑚−Δ𝑡^𝑁_𝑗=1𝜆_𝑗𝑡_𝑚+1𝜙_𝑗^′′𝑥^∗. (14)

4. Numerical Experiments and Discussions

In this section, we give an example to illustrate the feasibility of the method. In the following experiments, we use the Gauss radial basis function and let 𝑐=0.1.

Example. Consider the equation 𝑢_𝑡(𝑥,𝑡)=𝑢_𝑥𝑥(𝑥,𝑡)+𝑓(𝑡), with the conditions 𝑢()(𝑥)𝑢()𝑥,0=𝜑=2+cos𝑥,0,𝑡=ℎ₀(𝑡)=(3+𝑡)𝑒^−𝑡,𝑢(𝑙,𝑡)=ℎ_𝑙(𝑡)=(2+𝑡+cos𝑙)𝑒^−𝑡,𝑢𝑥^∗,𝑡=𝐸(𝑡)=2+𝑡+cos𝑥^∗𝑒^−𝑡, (15) and we let 𝑙=2, 𝑇=2, and 𝑥^∗=1.

The exact solutions are 𝑢(𝑥,𝑡)=(2+𝑡+cos𝑥)𝑒^−𝑡,𝑓(𝑡)=−(1+𝑡)𝑒^−𝑡. (16)

At first, in order to get the approximation effect of 𝑓(𝑡) according to different nodes of 𝑡, we let Δ𝑥=0.5 and Δ𝑡=0.05,0.005,0.0005, respectively. The exact solution 𝑓(𝑡), the numerical solution 𝑓(𝑡), and the error 𝑓(𝑡)−𝑓(𝑡) are plotted in Figures 1, 2, and 3.

(a)

(b)

Figure 1 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) (a) and the error 𝑓(𝑡)−𝑓(𝑡) (b) with Δ𝑥=0.5, Δ𝑡=0.05.

(a)

(b)

Figure 2 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) (a) and the error 𝑓(𝑡)−𝑓(𝑡) (b) with Δ𝑥=0.5, Δ𝑡=0.005.

(a)

(b)

Figure 3 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) (a) and the error 𝑓(𝑡)−𝑓(𝑡) (b) with Δ𝑥=0.5, Δ𝑡=0.0005.

As seen from the Figures 1–3, the error decreases and the approximation effect becomes better as Δ𝑡 decreases.

In Figure 4, we plot the exact solution 𝑢 and the numerical solutions 𝑢 at 𝑡=0.1,1,1.5,2, with Δ𝑥=0.5 and Δ𝑡=0.005.

(a)

(b)

(c)

(d)

Figure 4 

The exact solution 𝑢(𝑥,𝑡) and numerical solution 𝑢(𝑥,𝑡) at 𝑡=0.1,1,1.5,2, with Δ𝑥=0.5, Δ𝑡=0.005.

Figure 4 shows clearly that the approximate degree of numerical solution 𝑢 is very good for different value of 𝑡.

Figures 5 and 6 give the numerical solution 𝑢 and the error 𝑢−𝑢 with Δ𝑥=0.5, Δ𝑡=0.005,0.0005, respectively.

(a)

(b)

Figure 5 

The numerical solution of 𝑢(𝑥,𝑡) (a) and the error 𝑢(𝑥,𝑡)−𝑢(𝑥,𝑡) (b) with Δ𝑥=0.5, Δ𝑡=0.005.

(a)

(b)

Figure 6 

The numerical solution of 𝑢(𝑥,𝑡) (a) and the error 𝑢(𝑥,𝑡)−𝑢(𝑥,𝑡) (b) with Δ𝑥=0.5, Δ𝑡=0.0005.

As seen from the Figures 5 and 6, the effect of different Δ𝑡 to the approximate solution 𝑢 is the same as 𝑓; that is, the error decreases according to decreasing Δ𝑡.

Next, in order to illustrate the effect of the noisy data on numerical solution, we introduce artificial error function into the additional specification data as follows: 𝐸_𝛾(𝑡)=𝐸(𝑡)1+𝛾𝑡,𝑞, (17) where 𝛾(𝑡,𝑞) represents the noisy parameter, and we consider two cases in this paper:(i)𝛾(𝑡,𝑞)=const=𝛾, (ii)𝛾(𝑡,𝑞) is a random function of 𝑡 uniformly distributed on (−𝑞,𝑞).

We first consider case (i), Figures 7 and 8 plot the exact solution 𝑓(𝑡), the numerical solution 𝑓(𝑡), and the error 𝑓(𝑡)−𝑓(𝑡) with Δ𝑥=0.5, Δ𝑡=0.05, and 𝛾=0.01,0.001, respectively.

(a)

(b)

Figure 7 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) (a) and the error 𝑓(𝑡)−𝑓(𝑡) (b) with Δ𝑥=0.5, Δ𝑡=0.05, and 𝛾=0.01.

(a)

(b)

Figure 8 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) (a) and the error 𝑓(𝑡)−𝑓(𝑡) (b) with Δ𝑥=0.5, Δ𝑡=0.05, and 𝛾=0.001.

In case (ii), let the parameter 𝑞=0.002 and 0.001, respectively; we plot the exact solution 𝑓(𝑡), the numerical solution 𝑓(𝑡), and the error 𝑓(𝑡)−𝑓(𝑡) with Δ𝑥=0.5, Δ𝑡=0.05, in Figures 9 and 10.

(a)

(b)

Figure 9 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) (a) and the error 𝑓(𝑡)−𝑓(𝑡) (b) with Δ𝑥=0.5, Δ𝑡=0.05, and 𝑞=0.002.

(a)

(b)

Figure 10 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) (a) and the error 𝑓(𝑡)−𝑓(𝑡) (b) with Δ𝑥=0.5, Δ𝑡=0.05, and 𝑞=0.001.

The results shown in Figures 7–10 illustrate sensitivity of the inverse problem solution with respect to the noisy data, but when the noisy parameter is const, the result is relatively stable, with the parameter 𝛾 decreasing and the error decreasing too. When the noisy parameter is a random function, the result we obtain is relatively worse, but from the figures we can see that the error can be controlled in a certain range.

In the last experiment, in order to illustrate the approximate effect, we define 𝐸𝑢₂=∑^𝑀_𝑖=1∑^𝑁_𝑗=1𝑢𝑥_𝑖,𝑡_𝑗−𝑢𝑥_𝑖,𝑡_𝑗²,𝑀𝑁𝐸𝑓₂=∑^𝑁_𝑗=1𝑓𝑡_𝑗−𝑓𝑡_𝑗²𝑁,𝐸𝑢_∞=max_{1≤𝑖≤𝑀1≤𝑗≤𝑁}𝑢𝑥_𝑖,𝑡_𝑗−𝑢𝑥_𝑖,𝑡_𝑗,𝐸𝑓_∞=max_{1≤𝑖≤𝑀1≤𝑗≤𝑁}𝑓𝑡_𝑗−𝑓𝑡_𝑗, (18) where 𝑢(𝑥_𝑖,𝑡_𝑗) and 𝑓(𝑡_𝑗) are the exact solution and 𝑢(𝑥_𝑖,𝑡_𝑗) and 𝑓(𝑡_𝑗) are the approximate numerical solution. The error results with no noisy data are listed in Table 2. In Tables 3, 4, and 5, we list the error results with different noisy data in both cases.

Table 2 

The errors for different Δ𝑡.

Table 3 

The errors for Δ𝑡=0.05 with different 𝛾.

Table 4 

The errors for Δ𝑡=0.005 with different 𝛾.

Table 5 

The errors for Δ𝑡=0.05 with different 𝑞.

5. Conclusion

In this paper, the meshless method based on the radial basis functions is used for solving the inverse problem of heat equation with the time-dependent source term. From the numerical experiment, we can see that this method is high-efficiency and stable with no noisy data. When the noisy parameter is constant, the result is worsening, but not much. The worst result is the case that the noisy parameter is a random function; however, the error can be controlled in a certain range.

Acknowledgments

The project is supported by the Hunan Provincial Innovation Foundation for Postgraduate CX2011B156 and by NNSF of China Grants nos. 11271117 and 11201136.