Moving Least Squares Method for a One-Dimensional Parabolic Inverse Problem

Abstract

This paper investigates the numerical solution of a class of one-dimensional inverse parabolic problems using the moving least squares approximation; the inverse problem is the determination of an unknown source term depending on time. The collocation method is used for solving the equation; some numerical experiments are presented and discussed to illustrate the stability and high efficiency of the method.

1. Introduction

The inverse problem of parabolic equations appears in a wide variety of physical and engineering settings, such as elasticity, plasticity, hydrology, material sciences, heat transfer, combustion systems, medical imaging, transport problems, and control theory. Many researchers solved this problem using different methods [1–10].

For many decades, many numerical methods have been developed to solve partial differential equations, such as finite difference method, finite elements method, and boundary elements method. However, these above methods are mesh dependent methods, and, for some problems such as large deformation, the mesh generation is very difficult, especially with the irregular geometries. As a result, the meshless method appeared.

Compared with these mesh dependent numerical methods, in the meshless methods, mesh generation on the spatial domain of the problem is not needed; this property is the main advantage of these techniques over the mesh dependent methods, and it is a technique 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 need for the geometry of the domain.

Many meshless methods have been discussed, such as smooth particle hydrodynamics method [11], the finite point method [12], diffuse element method [13], element-free Galerkin method [14], hp-clouds method [15], and meshless local Petrov-Galerkin method [16]. Moving least squares is very important for constructing the shape functions in meshless method, and it had been studied by a number of investigators, Lancaster and Salkauskas [17], Levin [18], Mirzaei and Dehghan [19], Armentano [20, 21], and so forth.

In this paper, a collocation meshless method based on moving least squares approximation is used for numerical solution of a class of inverse problems in parabolic equation. The rest of the paper is organized as follows. In Section 2, we will give an outline of the moving least squares. In Section 3, we solve the inverse problem using the meshless method based on the moving least squares. Numerical experiments will be given in Section 4 in order to illustrate the feasibility of the method.

2. An Outline of the Moving Least Squares (MLS)

The MLS method as an approximation and interpolating method has been introduced by Lancaster and Salkauskas.

Given data values 𝐮={𝑢_𝑗}^𝑁_𝑗=1 at nodes 𝑥_𝑗, 𝑁 is the number of scattered nodes, the MLS method produces a function 𝑢^ℎ that approximates data 𝐮 in a weighted square sense.

It can be defined as 𝑢^ℎ(𝐱)=^𝑚_𝑖=1𝑝_𝑖(𝐱)𝑎_𝑖(𝐱)=𝑃^𝑇(𝐱)𝐚(𝐱),∀𝐱∈Ω, (1) where 𝑚 is the number of terms in the basis, 𝑝_𝑖(𝐱) are the monomial basis functions, 𝑎_𝑖(𝐱) are the coefficients of the basis functions, and 𝑃^𝑇(𝐱)=𝑝₁(𝐱),𝑝₂(𝐱),…,𝑝_𝑚(𝐱,𝑎(𝐱))=[𝑎₁(𝐱),𝑎₂(𝐱),…,𝑎_𝑚(𝐱)]^𝑇. (2)

The unknown coefficient 𝐚(𝐱) is determined by minimizing the functional 𝐽, which is defined as 𝐽=^𝑛_𝑗=1𝜔𝐱−𝐱_𝑗𝑢^ℎ(𝐱)−𝑢_𝑗²=^𝑛_𝑗=1𝜔𝐱−𝐱_𝑗^𝑚_𝑖=1𝑝_𝑖(𝐱_𝑗)𝑎_𝑖(𝐱)−𝑢_𝑗²=^𝑛_𝑗=1𝜔𝐱−𝐱_𝑗𝑃^𝑇(𝐱_𝑗)𝐚(𝐱)−𝑢_𝑗², (3) where 𝑛 is the number of nodes in the support domain of the point 𝐱, 𝜔(𝐱−𝐱_𝑗) is the weight function, and 𝐱_𝑗 are the nodes in the influence domain of 𝐱.

Equation (3) can be rewritten in the vector form 𝐽=(𝑃𝐚−𝐮)^𝑇𝑊(𝑃𝐚−𝐮), (4) where 𝑝𝑃=¹𝐱₁𝑝₂𝐱₁⋯𝑝_𝑚𝐱₁𝑝₁𝐱₂𝑝₂𝐱₂⋯𝑝_𝑚𝐱₂𝑝⋮⋮⋮₁𝐱_𝑛𝑝₂𝐱_𝑛⋯𝑝_𝑚𝐱_𝑛,𝐮=𝑢₁,𝑢₂,…,𝑢_𝑛^𝑇,𝜔𝑊=𝐱−𝐱₁0⋯00𝜔𝐱−𝐱₂⋯0⋮⋮⋮00⋯𝜔𝐱−𝐱_𝑛. (5)

According to the above conditions (4)-(5), taking the derivative 𝐚(𝐱) to zero, we have 𝜕𝐽(𝐱)𝐚(𝐱)(𝐱)𝜕𝐚=𝐴−𝐵𝐮=0, (6) which is equivalent to 𝐴(𝐱)𝐚(𝐱)=𝐵(𝐱)𝐮, (7) where 𝐴(𝐱)=𝑃^𝑇𝐵(𝐱)𝑊𝑃,=𝑃^𝑇𝑊. (8) Then, we get 𝐚(𝐱)=𝐴⁻¹(𝐱)𝐵(𝐱)𝑢, (9) and, by substituting (9) into (1), we have 𝑢^ℎ(𝐱)=𝑃^𝑇(𝐱)𝐚(𝐱)=Φ^𝑇(𝐱)𝐮=^𝑁_𝑗=1𝜙_𝑗(𝐱)𝑢_𝑗, (10) where Φ^𝑇(𝐱)=𝑃^𝑇(𝐱)𝐴⁻¹(𝐱)𝐵(𝐱) (11) and 𝜙_𝑗(𝐱) is called the shape function.

3. The Inverse Problem and Its Numerical Solution

In this section, the collocation meshless method based on moving least squares approximation is implemented for the numerical solution of the inverse problem of parabolic equations, and the problem is as follows.

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

Assume that these known functions satisfy the compatibility conditions; the unique solvability of this inverse problem is presented in [5]. However, this inverse problem is ill-posed; that is to say, the solution is unstable, and any small change in the additional specification may result in a dramatic change in the numerical solution. To obtain a stable numerical solution, some regularization strategies should be applied.

According to (13), consider the following transformation [22]: 𝐸^′(𝑡)=𝑢_𝑥𝑥|_𝑥=𝑥∗+𝑓(𝑡), (14) and using (14), we get 𝑓(𝑡)=𝐸^′(𝑡)−𝑢_𝑥𝑥|_𝑥=𝑥∗. (15)

By substituting (14) into (12), (12) is transformed to the following equation: 𝑢_𝑡=𝑢_𝑥𝑥+𝐸^′(𝑡)−𝑢_𝑥𝑥|_𝑥=𝑥∗,(𝑥,𝑡)∈𝑄_𝑇,𝑢()(𝑥)(),𝑢()𝑥,0=𝜑,𝑥∈0,𝑙0,𝑡=ℎ₀(𝑡),𝑢(𝑙,𝑡)=ℎ_𝑙(𝑡),𝑡∈(0,𝑇). (16) So the inverse problem is transformed to a direct problem.

The approximate function 𝑢(𝑥,𝑡) of 𝑢(𝑥,𝑡) can be represented as 𝑢(𝑥,𝑡)=^𝑛_𝑗=1𝜆_𝑗(𝑡)𝜙_𝑗(𝑥), (17) where 𝜙_𝑗(𝑥) is the shape function described in Section 2.

By substituting (17) into (16), the system (16) can be rewritten as ^𝑛_𝑗=1𝜆^′_𝑗(𝑡)𝜙_𝑗(𝑥)=^𝑛_𝑗=1𝜆_𝑗(𝑡)𝜙^′′_𝑗(𝑥)+𝐸^′(𝑡)−^𝑛_𝑗=1𝜆_𝑗(𝑡)𝜙^′′_𝑗𝑥^∗,(𝑥,𝑡)∈𝑄_𝑇,_𝑛𝑗=1𝜆_𝑗(0)𝜙_𝑗(𝑥)=𝜑(𝑥),𝑥∈(0,𝑙_𝑛𝑗=1),𝜆_𝑗(𝑡)𝜙_𝑗(0)=ℎ₀(𝑡),^𝑛_𝑗=1𝜆_𝑗(𝑡)𝜙_𝑗(𝑙)=ℎ_𝑙(𝑡),𝑡∈(0,𝑇). (18)

It is well known that the step lengths in the standard numerical method have regularization effect; hence, we can choose appropriate step lengths to compute the derivative of 𝑡. We apply a simple one-step forward difference formula, and we get ^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚+1−𝜆_𝑗𝑡_𝑚𝜙Δ𝑡^𝑗(𝑥)=^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙^′′_𝑗(𝑥)+𝐸𝑡_𝑚+1−𝐸𝑡_𝑚−Δ𝑡^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙^′′_𝑗𝑥^∗,(𝑥,𝑡)∈𝑄_𝑇,_𝑛𝑗=1𝜆_𝑗(0)𝜙_𝑗(𝑥)=𝜑(𝑥),𝑥∈(0,𝑙_𝑛𝑗=1),𝜆_𝑗𝑡_𝑚𝜙_𝑗(0)=ℎ₀𝑡_𝑚,^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗(𝑙)=ℎ_𝑙𝑡_𝑚,𝑡∈(0,𝑇). (19) Then, from the first formula of (19), we obtain ^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚+1𝜙_𝑗(𝑥)=^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗(𝑥)+Δ𝑡^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙^′′_𝑗(𝑥)+𝐸𝑡_𝑚+1−𝐸𝑡_𝑚−Δ𝑡^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙^′′_𝑗𝑥^∗;(𝑥,𝑡)∈𝑄_𝑇, (20) by substituting each 𝑥_𝑘 for 𝑥 in (20), we get ^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚+1𝜙_𝑗𝑥_𝑘=^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗𝑥_𝑘+Δ𝑡^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙^′′_𝑗𝑥_𝑘+𝐸𝑡_𝑚+1−𝐸𝑡_𝑚−Δ𝑡^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙^′′_𝑗𝑥^∗.(𝑥,𝑡)∈𝑄_𝑇. (21)

Let 𝑢𝑥_𝑘,𝑡_𝑚=^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙_𝑗𝑥_𝑘; (22) then (21) can be rewritten as follows: 𝑢𝑥_𝑘,𝑡_𝑚+1=^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚+1𝜙_𝑗𝑥_𝑘=𝑢𝑥_𝑘,𝑡_𝑚+Δ𝑡^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙^′′_𝑗𝑥_𝑘+𝐸𝑡_𝑚+1−𝐸𝑡_𝑚−Δ𝑡^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚𝜙^′′_𝑗𝑥^∗; (23) using the iterative method, by (19) and (22), we can obtain 𝜆_𝑗(𝑡_𝑚),  𝑗=1,2,…,𝑛,  𝑚=1,2,…,𝑀, and the numerical solution 𝑢(𝑥_𝑘,𝑡_𝑚); then we get the numerical solution 𝑓𝑡_𝑚+1=𝐸𝑡_𝑚+1−𝐸𝑡_𝑚−Δ𝑡^𝑛_𝑗=1𝜆_𝑗𝑡_𝑚+1𝜙^′′_𝑗𝑥^∗. (24)

4. Numerical Experiments and Discussions

In this section, we give two examples to illustrate the correctness of the theoretical result and the feasibility of the method.

Example 1. Consider (12)-(13), with the conditions 𝑢()(𝑥)𝑢()𝑥,0=𝜑=2+sin𝑥,0,𝑡=ℎ₀(𝑡)=(2+𝑡)𝑒^−𝑡,𝑢(𝑙,𝑡)=ℎ_𝑙(𝑡)=(2+𝑡+sin𝑙)𝑒^−𝑡,𝑢𝑥^∗,𝑡=𝐸(𝑡)=2+𝑡+sin𝑥^∗𝑒^−𝑡, (25) and we let 𝑙=2, 𝑇=2, and 𝑥^∗=1.

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

At first, in order to get the approximation effect of 𝑓(𝑡) according to different nodes of 𝑥 and 𝑡, we let Δ𝑥=0.5,  Δ𝑡=0.01,0.001 and Δ𝑥=0.1,  Δ𝑡=0.01,0.001, respectively; the exact solution 𝑓(𝑡), the numerical solution 𝑓(𝑡), and the error 𝑓(𝑡)−𝑓(𝑡) are plotted in Figures 1–4.

(a)

(b)

Figure 1 

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

(a)

(b)

Figure 2 

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

(a)

(b)

Figure 3 

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

(a)

(b)

Figure 4 

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

As seen from Figures 1 and 2 if the Δ𝑥 is fixed, the error decreases as Δ𝑡 decreases, but around the initial point, the error changes a little. When we decrease Δ𝑥, from Figures 3 and 4, it is clear that the error around the initial point decreases obviously, and with the decrease of Δ𝑡, the total error decreases too.

In Figure 5, we plot the exact solution 𝑢 and the numerical solutions 𝑢 at 𝑡=0.1,1,1.5,2, respectively, with Δ𝑥=0.1 and Δ𝑡=0.001.

(a)

(b)

(c)

(d)

Figure 5 

The exact solution 𝑢(𝑥,𝑡) and numerical solution 𝑢(𝑥,𝑡) at 𝑡=0.1 (a), 1 (b), 1.5 (c), and 2 (d).

Figure 5 shows clearly that the approximate solution 𝑢 is well for different value of 𝑡.

Figures 6 and 7 give the numerical solution 𝑢 and the error 𝑢−𝑢 with Δ𝑥=0.1,  Δ𝑡=0.01,0.001, respectively.

(a)

(b)

Figure 6 

The numerical solution of 𝑢(𝑥,𝑡) (a), the error 𝑢(𝑥,𝑡)−𝑢(𝑥,𝑡) (b) with Δ𝑥=0.1,  Δ𝑡=0.01.

(a)

(b)

Figure 7 

The numerical solution of 𝑢(𝑥,𝑡) (a), the error 𝑢(𝑥,𝑡)−𝑢(𝑥,𝑡) (b) with Δ𝑥=0.1,  Δ𝑡=0.001.

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

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

Let 𝛾=0.01,0.001, respectively; in case (i), we plot the exact solution 𝑓(𝑡), the numerical solution 𝑓(𝑡), and the error 𝑓(𝑡)−𝑓(𝑡) with Δ𝑥=0.1,  Δ𝑡=0.001, in Figures 8 and 9.

(a)

(b)

Figure 8 

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

(a)

(b)

Figure 9 

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

From Figures 8 and 9, we can see that when the noisy data exist in additional specification, the effect of numerical solution is worse than when there is no noisy data. However, the effect of approximation becomes better with the noisy parameter 𝛾 decrease; this suggests that the method in this paper is stable.

In case (ii), we let the parameter 𝑞=0.01 and 0.001, respectively, and plot the exact solution 𝑓(𝑡), the numerical solution 𝑓(𝑡), and the error 𝑓(𝑡)−𝑓(𝑡) with Δ𝑥=0.1,  Δ𝑡=0.001, and the results are shown in Figures 10 and 11.

(a)

(b)

Figure 10 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) (a), the error 𝑓(𝑡)−𝑓(𝑡) (b) with 𝑞=0.01.

(a)

(b)

Figure 11 

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

Figures 10 and 11 show that when the noisy parameter is random, the result of numerical solution is also bad, but it is the same to the case that 𝛾 is const.; with the parameter 𝑞 decrease, the effect of approximation becomes better. Moreover, comparing with the fact that 𝛾 is const, when the noisy parameter is random, although the result fluctuation is very big, they can be controlled in a certain range.

At last, in order to illustrate the approximate effect, we define the root mean square error of the functions 𝑢 and 𝑓 as follows: 𝐸𝑢=∑^𝑀_𝑖=1∑^𝑁_𝑗=1𝑢(𝑥_𝑖,𝑡_𝑗)−𝑢(𝑥_𝑖,𝑡_𝑗)²,𝑀𝑁𝐸𝑓=∑^𝑁_𝑗=1𝑓(𝑡_𝑗)−𝑓(𝑡_𝑗)²𝑁, (28) where 𝑢(𝑥_𝑖,𝑡_𝑗) and 𝑓(𝑡_𝑗) are the exact solution and 𝑢(𝑥_𝑖,𝑡_𝑗) and 𝑓(𝑡_𝑗) are the approximate numerical solution. The error results for different Δ𝑥 and Δ𝑡 with no noisy data are listed in Table 1.

Table 1 

The errors for different Δ𝑥 and Δ𝑡.

We will give the error results for different noisy parameter with Δ𝑥=0.1,  Δ𝑡=0.001 in Tables 2 and 3 and for different Δ𝑡 with the same noisy parameter in Table 4.

Table 2 

The errors for different 𝛾.

Table 3 

The errors for different 𝑞.

Table 4 

The errors for different Δ𝑡.

From Tables 2 and 3, we see that the errors vary with the noisy parameter. In Table 4, we find that the approximation effect varies with Δ𝑡 for the same noisy parameter, and these show that the length of Δ𝑡 has the effect of regularization, and our method in this paper is stable.

In the next example, we will compare the numerical solution using the finite difference method (FDM) and the meshless collocation method based on moving least squares (MLS), respectively.

Example 2. Consider the conditions 𝑢(𝑥,0)=𝜑(𝑥)=𝑥²,𝑢(0,𝑡)=ℎ₀(𝑡)(),𝑢()=2𝑡+sin𝜋𝑡𝑙,𝑡=ℎ_𝑙(𝑡)(),𝑢𝑥=4+2𝑡+sin𝜋𝑡_∗,𝑡=𝐸(𝑡)=1+2𝑡+sin(𝜋𝑡), (29) and we let 𝑙=2,  𝑇=2,  and  𝑥^∗=1.

The exact solutions are 𝑢(𝑥,𝑡)=𝑥²(),𝑓(𝑡)().+2𝑡+sin𝜋𝑡=𝜋cos𝜋𝑡 (30)

We plot the exact solution and numerical solution of 𝑓(𝑡) in three different cases using FDM and MLS, respectively.

Let Δ𝑥=0.1,  Δ𝑡=0.01; first, we plot the result in Figure 12 with no noisy data.

(a)

(b)

Figure 12 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) with no noisy data, (a) MLS and (b) FDM.

In Figures 13 and 14, we let 𝛾=0.001 and 𝑞=0.001, respectively.

(a)

(b)

Figure 13 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) with 𝛾=0.001, (a) MLS and (b) FDM.

(a)

(b)

Figure 14 

The exact solution 𝑓(𝑡) and the numerical solution of 𝑓(𝑡) with 𝑞=0.001, (a) MLS and (b) FDM.

From Figures 12, 13, and 14, we can see that when there is no noisy data, the two methods are stable relatively, and, when there exist noisy data, the method FDM is unstable, especially in the case that the noisy parameter is random; but the method MLS is always stable.

Finally, we give the errors of numerical solution with different cases using these two methods in Table 5.

Table 5 

The errors for different Δ𝑥, Δ𝑡 and noisy parameter.

From Table 5, we obtain that the method MLS is more stable and efficient than the method FDM in any case.

5. Conclusion

In this paper, the meshless method based on the moving least squares is used for solving the inverse problem of parabolic equation with the time-dependent source term. We use some transformations to change the problem, and, in order to overcome the ill-posedness, we use the appropriate length of time. From the experiments, we find that the numerical solution is stable. Moreover, compared with the method FDM, we see that our method is stable and highly efficient.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the Hunan Provincial Innovation Foundation for Postgraduate (CX2011B156) and NNSF of China (nos. 11271117 and 11201136).