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

Authors: Baiyu Wang; Anping Liao
Journal: Chinese Journal of Mathematics
 (2013)

---

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


---

## Body

## 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 Section2, 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.Letr be the distance between a fixed point x*∈Rd and any x∈Rd; that is, r=∥x-x*∥2.Defineϕ:R+→R, and ϕ(∥x-x*∥2) is a univariate function.Letx1,x2,…xN be a given set of distinct points. According to the definition of the function ϕ, we denote ϕj(x)=ϕ(∥x-xj∥2). For a given data {xj,uj},j=1,2,…,N, the interpolating function u(x) can be written as follows:
(1)u(x)=∑j=1N‍λjϕ(∥x-xj∥2),x∈Rd,
where λj is the unknown RBF coefficient and N is the number of data points.In order to computeλj, assume that we want to interpolate the values u(xj); that is,
(2)u(xj)=uj.By substitutingxi(i=1,2,…,N) for x in (1), we obtain the following linear system. 						(3)ΦΛ=U,
where
(4)Φ=[ϕ1(x1)ϕ2(x1)⋯ϕN(x1)ϕ1(x2)ϕ2(x2)⋯ϕN(x2)⋮⋮⋮ϕ1(xN)ϕ2(xN)⋯ϕN(xN)],Λ=[λ1λ2⋮λN],U=[u1u2⋮uN],
so we can get the unknown coefficient Λ from the linear system ΦΛ=U. 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 Table1.Table 1
Some well-known functions that generate RBFs.
Definition
Gaussian (GA)
ϕ
(
r
)
=
exp
(
-
c
r
2
)
Hardy multiquadrics (MQ)
ϕ
(
r
)
=
r
2
+
c
2
Inverse multiquadrics (IMQ)
ϕ
(
r
)
=
(
r
2
+
c
2
)
-
1
Inverse quadric (IQ)
ϕ
(
r
)
=
(
r
2
+
c
2
)
-
1


## 3. The Meshless Method for the Inverse Problem

In this section, we consider a class of inverse problems as follows.Findu=u(x,t) and f(t) which satisfy
(5)ut(x,t)=uxx(x,t)+f(t),(x,t)∈QT=(0,l)×(0,T),u(x,0)=φ(x),x∈(0,l),u(0,t)=h0(t),u(l,t)=hl(t),t∈(0,T),
and the additional specification
(6)u(x*,t)=E(t),x*∈(0,l),t∈(0,T),
where φ(x), h0(t), hl(t), and E(t)≠0 are known functions and x* is a fixed prescribed point in the admissible set.If the functionsφ(x), h0(t), hl(t) and E(t) 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(u,f) solves the inverse problem (5)-(6), applying the trace type functional formulation, we have
(7)E′(t)=uxx|x=x*+f(t).
From this, we get
(8)f(t)=E′(t)-uxx|x=x*.Substituting (8) into (5), we have the following initial boundary value problem:
(9)ut=uxx+(E′(t)-uxx|x=x*),(x,t)∈QT,u(x,0)=φ(x),x∈(0,l),u(0,t)=h0(t),u(l,t)=hl(t),t∈(0,T).The approximate functionu^(x,t) of u(x,t) can be represented as
(10)u^(x,t)=∑j=1N‍λj(t)ϕj(x),
where ϕj(x) is the radial basis functions described in Section 2.Substituting (10) into (9), the system (9) can be rewritten as
(11)∑j=1N‍λj′(t)ϕj(x)=∑j=1N‍λj(t)ϕj′′(x)+(E′(t)-∑j=1N‍λj(t)ϕj′′(x*)),(x,t)∈QT,∑j=1N‍λj(0)ϕj(x)=φ(x),x∈(0,l),∑j=1N‍λj(t)ϕj(0)=h0(t),∑j=1N‍λj(t)ϕj(l)=hl(t),t∈(0,T).For the solution of the system (11), we apply a simple one-step forward difference formula to the time and obtain
(12)∑j=1N‍λj(tm+1)-λj(tm)Δtϕj(x)=∑j=1n‍λj(tm)ϕj′′(x)+(E(tm+1)-E(tm)Δt-∑j=1N‍λj(tm)ϕj′′(x*)),(x,t)∈QT,∑j=1N‍λj(0)ϕj(x)=φ(x),x∈(0,l),∑j=1N‍λj(tm)ϕj(0)=h0(tm),∑j=1N‍λj(tm)ϕj(l)=hl(tm),t∈(0,T).From the iterative method, we can obtainλj(tm), j=1,2,…,N, m=1,2,…,M, and combining (10) and (8), we get the numerical solution as follows:
(13)u^(x,tm+1)=∑j=1N‍λj(tm+1)ϕj(x)=∑j=1N‍λj(tm)ϕj(x)+Δt∑j=1N‍λj(tm)ϕj′′(x)+E(tm+1)-E(tm)-Δt∑j=1N‍λj(tm)ϕj′′(x*)
and the approximate solution as follows:
(14)f^(tm+1)=E(tm+1)-E(tm)Δt-∑j=1N‍λj(tm+1)ϕj′′(x*).


## 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 letc=0.1.Example. Consider the equationut(x,t)=uxx(x,t)+f(t), with the conditions
(15)u(x,0)=φ(x)=2+cosx,u(0,t)=h0(t)=(3+t)e-t,u(l,t)=hl(t)=(2+t+cosl)e-t,u(x*,t)=E(t)=(2+t+cosx*)e-t,
and we let l=2, T=2, and x*=1.The exact solutions are(16)u(x,t)=(2+t+cosx)e-t,f(t)=-(1+t)e-t.At first, in order to get the approximation effect off^(t) according to different nodes of t, we let Δx=0.5 and Δt=0.05,0.005,0.0005, respectively. The exact solution f(t), the numerical solution f^(t), and the error f(t)-f^(t) are plotted in Figures 1, 2, and 3.The exact solutionf(t) and the numerical solution of f^(t) (a) and the error f(t)-f^(t) (b) with Δx=0.5, Δt=0.05.
(a)
(b)The exact solutionf(t) and the numerical solution of f^(t) (a) and the error f(t)-f^(t) (b) with Δx=0.5, Δt=0.005.
(a)
(b)The exact solutionf(t) and the numerical solution of f^(t) (a) and the error f(t)-f^(t) (b) with Δx=0.5, Δt=0.0005.
(a)
(b)As seen from the Figures1–3, the error decreases and the approximation effect becomes better as Δt decreases.In Figure4, we plot the exact solution u and the numerical solutions u^ at t=0.1,1,1.5,2, with Δx=0.5 and Δt=0.005.The exact solutionu(x,t) and numerical solution u^(x,t) at t=0.1,1,1.5,2, with Δx=0.5, Δt=0.005.
(a)
(b)
(c)
(d)Figure4 shows clearly that the approximate degree of numerical solution u^ is very good for different value of t.Figures5 and 6 give the numerical solution u^ and the error u^-u with Δx=0.5, Δt=0.005,0.0005, respectively.The numerical solution ofu^(x,t) (a) and the error u(x,t)-u^(x,t) (b) with Δx=0.5, Δt=0.005.
(a)
(b)The numerical solution ofu^(x,t) (a) and the error u(x,t)-u^(x,t) (b) with Δx=0.5, Δt=0.0005.
(a)
(b)As seen from the Figures5 and 6, the effect of different Δt to the approximate solution u^ is the same as f^; that is, the error decreases according to decreasing Δt.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:(17)Eγ(t)=E(t)(1+γ(t,q)),
where γ(t,q) represents the noisy parameter, and we consider two cases in this paper:(i)
γ(t,q)=const=γ,(ii)
γ(t,q) is a random function of t uniformly distributed on (-q,q).We first consider case (i), Figures7 and 8 plot the exact solution f(t), the numerical solution f^(t), and the error f(t)-f^(t) with Δx=0.5, Δt=0.05, and γ=0.01,0.001, respectively.The exact solutionf(t) and the numerical solution of f^(t) (a) and the error f(t)-f^(t) (b) with Δx=0.5, Δt=0.05, and γ=0.01.
(a)
(b)The exact solutionf(t) and the numerical solution of f^(t) (a) and the error f(t)-f^(t) (b) with Δx=0.5, Δt=0.05, and γ=0.001.
(a)
(b)In case (ii), let the parameterq=0.002 and 0.001, respectively; we plot the exact solution f(t), the numerical solution f^(t), and the error f(t)-f^(t) with Δx=0.5, Δt=0.05, in Figures 9 and 10.The exact solutionf(t) and the numerical solution of f^(t) (a) and the error f(t)-f^(t) (b) with Δx=0.5, Δt=0.05, and q=0.002.
(a)
(b)The exact solutionf(t) and the numerical solution of f^(t) (a) and the error f(t)-f^(t) (b) with Δx=0.5, Δt=0.05, and q=0.001.
(a)
(b)The results shown in Figures7–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(18)Eu2=∑i=1M∑j=1N(u(xi,tj)-u^(xi,tj))2MN,Ef2=∑j=1N(f(tj)-f^(tj))2N,Eu∞=max1≤i≤M1≤j≤N∥u(xi,tj)-u^(xi,tj)∥,Ef∞=max1≤i≤M1≤j≤N∥f(tj)-f^(tj)∥,
where u(xi,tj) and f(tj) are the exact solution and u^(xi,tj) and f^(tj) 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Δt.
Δ
t
E
u
2
E
u
∞
E
f
2
E
f
∞
0.05
4.7798
×
1
0
-
4
1.3944
×
1
0
-
3
5.6777
×
1
0
-
3
1.4795
×
1
0
-
2
0.01
1.8430
×
1
0
-
4
5.3488
×
1
0
-
4
1.0676
×
1
0
-
3
2.9989
×
1
0
-
3
0.005
1.4814
×
1
0
-
4
4.3282
×
1
0
-
4
4.8902
×
1
0
-
4
1.6424
×
1
0
-
3
0.001
1.1931
×
1
0
-
4
3.5192
×
1
0
-
4
3.9699
×
1
0
-
5
5.6000
×
1
0
-
4
0.0005
1.1571
×
1
0
-
4
3.4185
×
1
0
-
4
4.6209
×
1
0
-
5
4.2486
×
1
0
-
4
0.0001
1.1283
×
1
0
-
4
3.3380
×
1
0
-
4
8.6191
×
1
0
-
5
3.1679
×
1
0
-
4Table 3
The errors forΔt=0.05 with different γ.
γ
E
u
2
E
u
∞
E
f
2
E
f
∞
0.01
8.3499
×
1
0
-
3
1.9258
×
1
0
-
2
9.7874
×
1
0
-
3
5.1589
×
1
0
-
2
0.001
9.4176
×
1
0
-
4
1.9258
×
1
0
-
3
4.4862
×
1
0
-
3
1.3189
×
1
0
-
2Table 4
The errors forΔt=0.005 with different γ.
γ
E
u
2
E
u
∞
E
f
2
E
f
∞
0.01
8.3369
×
1
0
-
3
1.9258
×
1
0
-
2
1.3919
×
1
0
-
2
4.3541
×
1
0
-
2
0.001
8.3550
×
1
0
-
4
1.9258
×
1
0
-
3
1.0025
×
1
0
-
3
5.2231
×
1
0
-
3Table 5
The errors forΔt=0.05 with different q.
q
E
u
2
E
u
∞
E
f
2
E
f
∞
0.002
2.2191
×
1
0
-
3
5.8045
×
1
0
-
3
1.8645
×
1
0
-
2
1.1786
×
1
0
-
1
0.001
8.7714
×
1
0
-
4
3.5146
×
1
0
-
3
1.0989
×
1
0
-
2
5.9674
×
1
0
-
2


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



---

*Source: 761272-2013-09-12.xml*