399
Views
2
CrossRef citations to date
0
Altmetric
Original Articles

A class of multistep numerical difference schemes applied in inverse heat conduction problem with a control parameter

&
Pages 887-942 | Received 06 Jun 2017, Accepted 12 Jun 2018, Published online: 01 Aug 2018

ABSTRACT

As this paper is concerned with a class of multistep numerical difference techniques to solve one-dimensional parabolic inverse problems with source control parameter p(t), we apply the linear multistep method combining with Lagrange interpolation to develop three different numerical difference schemes. The problem of numerical differentiation with noisy scattered data is mildly ill-posed, the smoothing spline model based on Tikhonov regularization method is developed to compute numerical derivative contaminated by noise error. Simultaneously, the truncation error estimations and the convergence conclusions are proposed for the above difference methods respectively. The results of numerical tests with different noise levels are given to show that the presented algorithms are accurate and effective.

2010 MSC CLASSIFICATIONS:

1. Introduction

Inverse heat conduction problems (IHCPs) arise in many fields of physics such as control modelling with heat propagation, mechanics, heat flux, the remote sensing technology, signal processing and industrial controlling. Many researchers have been attracted to study this field because their theories possess distinct novelty and are applicable to a broad range of applications.

It is known that heat conduction process is very smooth, but the studies of IHCPs are more difficult for its irreversible. The IHCP is ill-posed in some sense that any small errors for measured values on some specified points can result in a large error to the solution [Citation1], which implies that the resultant discretized matrix is ill-conditioned, i.e. that the solution is potentially sensitive to perturbations. So the key is to find a stable algorithm to overcome the ill-posedness of problems, using a reasonable regularization method generally [Citation2,Citation3]. Furthermore, Tikhonov regularization method plays an important role in efficiently solving the ill-posed problems [Citation4].

Many numerical algorithms for solving one-dimensional IHCP have been proposed during the last decade, including the finite difference method (FDM) [Citation5–10], the radial basis function (RBF) method [Citation11–15], the method of fundamental solutions (MFS) [Citation1,Citation16,Citation17], the collocation method (CM) [Citation18] and the boundary element method (BEM) used in the determination of a spacewise-dependent heat source [Citation19], etc. Meanwhile, some numerical technologies are applied in multidimensional parabolic equations [Citation20–23], Dehghan et al. [Citation24–31] proposed the meshless method of RBFs to solve a series of nonlinear high dimensional partial differential equations, and applied other schemes to solve inverse parabolic problems to obtain effective results. We consider that the main advantage of the forward FDM, whose form is simple and easy to compute. Moreover, when we handle the one-dimensional inverse parabolic equations, it is different between the spatial domain and temporal domain in the physical view [Citation14]. For an aspect of constructing tools with high accuracy, we apply the linear multistep technology to develop more accurate difference schemes. In practical problems, the overspecification node x cannot just be chosen a mesh point xi, so we use cubic or quartic Lagrange polynomial interpolation method to approximate function u~(x,t).

This paper is organized as follows. In Section 2, we present a class of IHCP with a source control parameter p(t) and introduce the main idea of transforming it into a direct problem to solve. In Section 3, first we describe the process of constructing numerical finite difference schemes by using linear multistep technologies. Second, we give three different FDMs. Furthermore, we deduce the order results of truncation error and the convergence estimates, including the case of E(t) with noise. The numerical results for the test used are given in Section 4. Finally we summarize this paper in Section 5.

2. Problem statement

We now consider a class of IHCP with the source control parameter p(t), which is stated as follows: (1) wt=2wx2+p(t)w+ϕ(x,t),0<x<1,0<tT,(1) with the initial condition (2) w(x,0)=f(x),0x1,(2) the boundary conditions (3) w(0,t)=g0(t),0tT,(3) (4) w(1,t)=g1(t),0tT,(4) and the overspecification at a node x(0,1) in the spatial domain (5) w(x,t)=E(t),0tT,(5) where ϕ(x,t),f(x),g0(t),g1(t) and E(t)(E(t)0) are the given functions, but the functions w(x,t) and p(t) are unknown. The function w(x,t) denotes temperature distribution which is the exact solution to the IHCP, and p(t) is a source control parameter, where (Equation5) produces a desired temperature at each time t and a given point x in the spatial domain. The system of equations (Equation1)–(Equation5) represents a heat transfer process, and the existence, uniqueness of the solution and some applications of this problem are introduced in [Citation32–35]. The aim of solving the problem is to identify the function w(x,t) and the parameter p(t). Since the formula (Equation1) includes two unknown functions w(x,t) and p(t), so we adopt the following treatment by J.R. Cannon [Citation33,Citation36], which can transform Equations (Equation1)–(Equation4) into another direct parabolic problem.

Let (6) r(t)=e0tp(s)ds(6) and (7) u(x,t)=r(t)w(x,t),(7) which can eliminate the nonlinear term p(t)w from Equation (Equation1) and to obtain a nonlocal parabolic equation, then the problems (Equation1)–(Equation4) can be transformed into the following equations: (8) ut=2ux2+r(t)ϕ(x,t),0<x<1,0<tT,(8) with the initial condition (9) u(x,0)=f(x),0x1,(9) and the boundary conditions (10) u(0,t)=r(t)g0(t),0tT,(10) (11) u(1,t)=r(t)g1(t),0tT.(11) By (Equation5) and (Equation7), r(t) can be obtained by (12) r(t)=u(x,t)E(t),(12) and w(x,t) can be presented by (13) w(x,t)=u(x,t)r(t),(13) then we can obtain the control parameter p(t) as (14) p(t)=r(t)r(t),(14) where the superscript represents the first derivative of some function in this paper. In addition, the general algorithms of this paper consisting of two subproblems are as follows:

A1. First, we compute the numerical solutions u~(xi,tj) (or u~δ(xi,tj)) of u(xi,tj) by multistep difference technologies and use Lagrange polynomial interpolation to compute u~(x,tj) (or u~δ(x,tj)), further to obtain r~(tj) (or r~δ(tj)) by (Equation12), then get w~(xi,tj) (or w~δ(xi,tj)) by (Equation13).

A2. For a set of given scattered data {r~(tj)}j=1n (or {r~δ(tj)}j=1n) above, we use the numerical differential scheme combining with Tikhonov regularization method to obtain the smooth regularized solution r~α(t) (or r~αδ(t)), further to compute p~α(tj)=((r~α(tj))/(r~(tj))) (or p~αδ(tj)=(((r~αδ)(tj))/(r~δ(tj)))) by (Equation14).

Remark 2.1

If the data in (Equation1)–(Equation5) are smooth and compatible, then the problem (Equation1)–(Equation5) is equivalent to the problem (Equation8)–(Equation12), and the solution of this parabolic problem exists uniquely, which is smooth from the standard regularization theory of parabolic equations. In addition, if u(x,t) is obtained numerically, the functions w(x,t) and p(t) can be computed by the forms (Equation13) and (Equation14) respectively, it means that the problem (Equation1)–(Equation5) is mildly ill-posed in [Citation8,Citation33,Citation37,Citation38].

Remark 2.2

When the measured data E(t) is exact in Section 3.1.3, the problem A1 above is well posed by Theorem 3.3. However, the problem A2 above is ill-posed by Theorem 3.4, in this case, the regularization parameter α1 depends on the errors from A1 with respect to the step sizes τ and h, i.e. if α1=τ+h2, the convergence can be obtained as τ,h0.

Remark 2.3

In Section 3.1.4, when E(t) has noise of level δ, i.e. |Eδ(t)E(t)|δ for all t[0,T], then the problems A1 and A2 above are both ill-posed. As shown in Theorem 3.6, the step sizes τ and h should be selected according to δ in order to obtain the order of convergence δ1/2 for w~δ(xi,tj) in A1. Furthermore, as shown in Theorem 3.8, the step sizes τ, h, and the regularization parameter α1δ are related to δ in A2, i.e. if τ=δ1/2 and α1δ=δ1/2, the order of convergence is δ(ν1)/4ν, ν2.

3. Multistep numerical difference schemes

We divide the domain [0,1]×[0,T] into an M×N uniform mesh with the spatial step-size h=1/M in the x direction and the temporal step-size τ=T/N, where M and N are positive integers. Let Ω={(xi,tj),xi=ih,tj=jτ,i=0,1,,M,j=0,1,,N}. We assume u~(xi,tj) is an approximation value of u(xi,tj), the values w~(xi,tj),r~(tj) and p~(tj) are the similar approximations of w(xi,tj),r(tj) and p(tj) respectively.

On the left side of Equation (Equation8), u/t can be substituted by (15) u(xi,tj)tu(xi,tj+1)u(xi,tj)τ(15) at the grid point (xi,tj), where i=0,1,,M,j=0,1,,N1. The formula 2u/x2 on the right side of Equation (Equation8) can be approximated by the linear multistep methods [Citation39], the key is to construct a k-step linear discretization equation as follows: (16) i=0kai2u(xn+i,tj)x21h2i=0kbiu(xn+i,tj),(16) where n=0,1,,Mk,j=0,1,,N, ai and bi are real constants. In order to determine the parameters ai and bi, we first define the following operator: (17) L[u(x,t);h]=i=0kbiu(x+ih,t)h2ai2u(x+ih,t)x2,(17) and expand u(x+ih,t) and (2u(x+ih,t))/x2 at point x by Taylor polynomials, then we can find (18) L[u(x,t);h]=C0u(x,t)+C1hu(x,t)x+C2h22u(x,t)x2++Cqhqqu(x,t)xq+,(18) where C0,C1,,Cq are constants and satisfy (19) C0=i=0kbi,C1=i=0kibi,C2=12!i=0ki2bii=0kai,Cq=1q!i=0kiqbi1(q2)!i=0kiq2ai,q=3,4,,(19) by choosing suitable values k, ai and bi such that C0=C1==Cq=0, and Cq+10, so we get the following truncation error: (20) L[u(x,t);h]=Cq+1hq+1q+1u(x,t)xq+1+O(hq+2).(20) By substituting u(x+ih,t) and (2u(x+ih,t))/x2 in Equation (Equation17) for u(xn+i,tj) and (2u(xn+i,tj))/x2, and we can obtain the k-step discrete equation  (Equation16) with its local truncated error Cq+1hq+1((q+1u(xn+i,tj))/(xq+1))+O(hq+2).

3.1. Three-point numerical difference scheme

3.1.1. Algorithm 1

Let k=2,j=1,2,,N, in Equations (Equation15)–(Equation16), and denote the grid ratio ρ=τ/h2. We choose an appropriate coefficient set {a0,a1,a2}, correspondingly obtain the derived set {b0,b1,b2} and the nodes-subscript set {i0,i1,i2} in Table , so we can get the general difference equation formula: (21) u~(xi,tj)=u~(xi,tj1)+ρ{b0u~(xi0,tj1)+b1u~(xi1,tj1)+b2u~(xi2,tj1)}+τr~(tj1)ϕ(xi,tj1).(21)

Table 1. The different configuration parameters in formula (Equation21).

Note that r~(t0)=1 and u~(xi,0)=f(xi),0iM by (Equation6) and (Equation9), respectively. In practice, the overspecification node x cannot just be chosen as a mesh point xi,i=0,1,,M, so we use cubic-Lagrange polynomial interpolation method to compute u~(x,tj), whose structure is simple, and to achieve the required accuracy of proposed methods. According to the location of x, there are three cases to be proposed as follows:

1) If x(x0,x1), then u~(x,tj) yields (22) u~(x,tj)=k=14u~(xk,tj)L3,k(x),1jN.(22) 2) If x[xi,xi+3], 1iM4, then u~(x,tj) yields (23) u~(x,tj)=k=ii+3u~(xk,tj)L3,k(x),1jN.(23) 3) If x(xM1,xM), then u~(x,tj) yields (24) u~(x,tj)=k=M4M1u~(xk,tj)L3,k(x),1jN,(24) the above L3,k(x) denotes a cubic-Lagrange basis function, and (25) L3,k(x)=ω4(x)(xxk)ω4(xk),(25) where ω4(x)=k=ii+3(xxk), xk indicates an interpolation node in the domain [x1,xM1]. Then according to the formula (Equation12), the approximation of r(tj) can be obtained by (26) r~(tj)=u~(x,tj)E(tj),1jN,(26) the boundary values u~(x0,tj) and u~(xM,tj) will be updated by Equations (Equation10) and (Equation11) as follows: (28) u~(x0,tj)=r~(tj)g0(tj)=g0(tj)E(tj)u~(x,tj),(28) (29) u~(xM,tj)=r~(tj)g1(tj)=g1(tj)E(tj)u~(x,tj),(29) where 1jN. According to the formula (Equation13), the approximation solution to w(xi,tj) can be computed by (29) w~(xi,tj)=u~(xi,tj)r~(tj),0iM,0jN.(29)

3.1.2. Algorithm 2

Since the values r~(tj) in Equation (Equation26) can only be derived by the discrete measured data with errors, we must apply the numerical differential technique to get the corresponding discrete data for the values p~(tj) by Equation (Equation14). However, the problem of numerical differentiation is well known to be ill-posed, it means that any small errors for the measurement values on some specified points can cause large errors to the implemented derivative [Citation40–43]. Thus it is the key aim to develop a stable scheme to deal with the ill-posed problem. The traditional FDM is simple and effective for precise data, but it is failed to compute accurate derivative of the measured values with noisy data, since the size of mesh depends on the noise level, i.e. which should not be too small, in other words, the number of measured points should not be chosen too big, otherwise, it may be very large for the final errors of the approximation solutions [Citation41,Citation44].

Recently, some numerical differential schemes based on the Tikhonov regularization technique and the smoothing spline have been proposed [Citation40–42,Citation44,Citation45], based on a suitable strategy of the regularization parameter, these differential methods with regularization are effective and stable. In this paper, we use numerical differentiation with regularization method presented in [Citation41,Citation42] to compute the values r~(tj).

For given a set of measured values {r~j}j=1n with the error data, and define Xν={rrCν1(R),r(ν)L2(R)}. Consider a smoothing functional on Xν: (30) Υν(r)=1nj=1nr(tj)r~j2+α1r(ν)L2(R)2,(30) where ν2 is a fixed integer, nν, let {tj}1n be a set of distinct points on [a,b], and α1 is a regularization parameter under Tikhonov regularization sense.

Lemma 3.1

[Citation41]

Define (31) r~α1(t)=j=1np¯j|ttj|2ν1+j=1νq¯jtj1,(31) where coefficients {p¯j}1n and {q¯j}1ν satisfy (32) r~α1(tj)+2(2ν1)!(1)να1np¯j=r~j,j=1,2,,n,(32) (33) j=1np¯jtji=0,i=0,1,,ν1,(33) then r~α1 is a unique minimizer of the smoothing functional (Equation30), i.e. Υν(r~α1)=minrXνΥν(r)Υν(y), yXν.

Lemma 3.2

[Citation41]

There is a unique solution for linear system of Equations (Equation32) and (Equation33).

By differentiating r~α1(t) of (Equation31) with respect to t, then we have (34) r~α1(t)=(2ν1)j=1np¯jsign(ttj)|ttj|2ν2+j=1ν(j1)q¯jtj2.(34) Furthermore, the approximation values p~α1(tj) of p(tj) can be obtained by (Equation14) (35) p~α1(tj)=r~α1(tj)r~(tj),0jN.(35) For convenience, the above approach is called 3PNDS in this paper.

3.1.3. The convergence estimates of 3PNDS for E(t) without noise

The proof of the convergence estimate for 3PNDS consists of the following details. In order to express convenience and simplicity, we define (36) Δh2u~ij=u~(xi+1,tj)+u~(xi1,tj)2u~(xi,tj)/h2,(36) and assume that G0j=g0(tj)/E(tj), GMj=g1(tj)/E(tj), Φij=ϕ(xi,tj)/E(tj), combining with Equations (Equation9)–(Equation12), (Equation21)–(Equation28) and (Equation36), then we have two corresponding systems for j=1 and j2 respectively as follows: (37) (u~(x0,t1)u~(x0,t0))/τ=Δh2u~10+ϕ(x0,t0),(u~(xi,t1)u~(xi,t0))/τ=Δh2u~i0+ϕ(xi,t0),1iM1,(u~(xM,t1)u~(xM,t0))/τ=Δh2u~M10+ϕ(xM,t0),u~(xi,0)=f(xi),0iM,u~(xi,t1)=Gi1u~(x,t1),i=0 or M (updateboundary),(37) and (38) (u~(x0,tj)u~(x0,tj1))/τ=Δh2u~1j1+Φ0j1u~(x,tj1),j2,(u~(xi,tj)u~(xi,tj1))/τ=Δh2u~ij1+Φij1u~(x,tj1),1iM1,j2,(u~(xM,tj)u~(xM,tj1))/τ=Δh2u~M1j1+ΦMj1u~(x,tj1),j2,u~(xi,0)=f(xi),0iM,u~(xi,tj)=Giju~(x,tj),i=0 or M, j2(updateboundary).(38) Suppose that eij=u(xi,tj)u~(xi,tj), ej=u(x,tj)u~(x,tj), according to Taylor's expansion, eij and ej satisfy the following formulas by Equations (Equation37) and (Equation38) (39) (e01e00)/τ=Δh2e10+ϵ~00,(ei1ei0)/τ=Δh2ei0+ϵ~i0,1iM1,(eM1eM0)/τ=Δh2eM10+ϵ~M0,ei0=0,0iM,ei1=Gi1e1,i=0 or M,(39) and (40) (e0je0j1)/τ=Δh2e1j1+Φ0j1ej1+ϵ0j1,j2,(eijeij1)/τ=Δh2eij1+Φij1ej1+ϵij1,1iM1,j2,(eMjeMj1)/τ=Δh2eM1j1+ΦMj1ej1+ϵMj1,j2,ei0=0,0iM,eij=Gijej,i=0 or M, j2,(40) where ϵ~00, ϵ~i0, ϵ~M0, ϵ0j1, ϵij1 and ϵMj1 are the truncation errors derived by the discretization of differential equation. In addition, if the given functions g0(t), g1(t), E(t) (E(t)0) and ϕ(x,t) are continuous on the interval [0,T] or [0,1]×[0,T], then the discretization values {Φij} and {Gij} are bounded, so there exists a positive constant K such that the following inequalities hold: max0iM,1jN1|Φij|K,maxi=0orM,1jN|Gij|K. The error bound of Lagrange polynomial is developed as follows:

Lemma 3.3

[Citation46]

Let x0,x1,,xn be distinct numbers in the interval [a,b] and P(x)Cn+1[a,b], and define ωn+1(x)=(xx0)(xx1)(xxn). Then, there exists a number ξ(x)(a,b) (generally unknown) between x0,x1,,xn, and satisfying the following formula: %P(x)=Ln(x)+P(n+1)(ξ(x))(n+1)!ωn+1(x), where Ln(x) is Lagrange interpolating polynomial of degree n given by %Ln(x)=P(x0)Ln,0(x)++P(xn)Ln,n(x)=k=0nP(xk)Ln,k(x), where Ln,k(x)=(ωn+1(x))/((xxk)ωn+1(xk)).

Theorem 3.1

Let u(x,t)C4,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,1], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] (E(t)0), then the truncation error order of 3PNDS can reach O(τ+h) at least.

Proof.

If the given point x[xs1,xs+2], where 1sM2, then we can obtain the following conclusion by Lemma 3.3: (41) u(x,t)=k=s1s+2u(xk,t)L3,k(x)+O((xxs1)(xxs)(xxs+1)(xxs+2)),(41) furthermore, x and t in Equation (Equation41) are replaced by x and tj, respectively, we have (42) u(x,tj)=k=s1s+2u(xk,tj)L3,k(x)+O((xxs1)(xxs)(xxs+1)(xxs+2)),=k=s1s+2u(xk,tj)L3,k(x)+O(h4),(42) where j=1,2,,N. We define the following linear operators: (43) L~u(x,t)=u(x,t)t2u(x,t)x2ϕ(x,t)E(t)u(x,t)(43) and (44) L~hu(x,t)=u(x,t+τ)u(x,t)τu(x+h,t)+u(xh,t)2u(x,t)h2ϕ(x,t)E(t)k=s1s+2u(xk,t)L3,k(x),(44) (45) L~0u(x,t)=u(x,t+τ)u(x,t)τu(x+2h,t)+u(x,t)2u(x+h,t)h2ϕ(x,t)E(t)k=s1s+2u(xk,t)L3,k(x),(45) (46) L~Mu(x,t)=u(x,t+τ)u(x,t)τu(x,t)+u(x2h,t)2u(xh,t)h2ϕ(x,t)E(t)k=s1s+2u(xk,t)L3,k(x).(46) Now we replace x and t in Equations (Equation43)–(Equation46) with xi and tj, respectively. In addition, we define the following truncation error operators Th, T0 and TM on u(xi,tj), here the truncation errors derive from two parts, i.e. one is caused by the truncation of infinite Taylor series, the other comes from the truncation of Lagrange interpolation, then we have (47) Thu(xi,tj)=L~hu(xi,tj)L~u(xi,tj),(47) (48) T0u(xi,tj)=L~0u(xi,tj)L~u(xi,tj),(48) (49) TMu(xi,tj)=L~Mu(xi,tj)L~u(xi,tj).(49) ♯ When j=0, since r(t0)=1 in Equation (Equation6). By the forms (Equation15)–(Equation20) and (Equation43)–(Equation46), combining with Taylor's expansion, we have (50) L~u(x,t)=u(x,t)t2u(x,t)x2ϕ(x,t)(50) and (51) L~hu(x,t)=u(x,t+τ)u(x,t)τu(x+h,t)+u(xh,t)2u(x,t)h2ϕ(x,t),(51) (52) L~0u(x,t)=u(x,t+τ)u(x,t)τu(x+2h,t)+u(x,t)2u(x+h,t)h2ϕ(x,t),(52) (53) L~Mu(x,t)=u(x,t+τ)u(x,t)τu(x,t)+u(x2h,t)2u(xh,t)h2ϕ(x,t).(53) By the formulas (Equation47)–(Equation49) and (Equation50)–(Equation53), we can obtain the truncation error orders as follows: (54) Thu(xi,tj)=O(τ+h2),1iM1,(54) (55) T0u(xi,tj)=O(τ+h),i=0,(55) (56) TMu(xi,tj)=O(τ+h),i=M,(56) from the above deduction, the forms (Equation54)–(Equation56) show that the truncation error order of 3PNDS can reach O(τ+h) at least for j=0.

♯ When 1jN1, let T be the truncation error operator on u(x,t), then we have the following truncation error estimation of u(x,tj) by Equation (Equation42): (57) Tu(x,tj)=u(x,tj)k=s1s+2u(xk,tj)L3,k(x)=O(h4),1jN.(57) According to the forms (Equation42)–(Equation46) and (Equation47)–(Equation49), we can obtain the following conclusions: (58) Thu(xi,tj)=O(τ+h2),1iM1,(58) (59) T0u(xi,tj)=O(τ+h),i=0,(59) (60) TMu(xi,tj)=O(τ+h),i=M,(60) where j=0,1,,N1. In short, by the results (Equation54)–(Equation56) and (Equation58)–(Equation60), the truncation error order of 3PDNS can reach O(τ+h) at least. According to the above arguments, the truncation error order is independent of the distribution interval of x.

The estimates of the truncation errors ϵ~i0 and ϵij are as follows:

Lemma 3.4

Suppose that u(x,t)C4,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, then there exists a positive constant Q~1 such that max1iM1|ϵ~i0|Q~1(τ+h2),max1iM1,1jN1|ϵij|Q~1(τ+h2),maxi={0,M}|ϵ~i0|Q~1(τ+h),maxi={0,M},1jN1|ϵij|Q~1(τ+h), where Q~1=max0x1,0tT{|2u/t2|,|3u/x3|,|4u/x4|} is independent of τ and h.

Proof.

The proof can be directly induced by the procedure of Theorem 3.1.

In practice, we can select a suitable small spatial-step h such that h<min{(|xx0|/2),(|xxM|/2)} for 3PDNS, thus the fixed point x will not be located in the intervals (x0,x1) and (xM1,xM), further, making it close to the middle of the interval (0,1). By (Equation39), (Equation40) and Lemma 3.4, it is shown that the truncated error order number can be improved from O(τ+h) to O(τ+h2), more details can be found in the following results.

Lemma 3.5

Let u(x,t)C4,2([0,1]×[0,T]), there exists a spatial-step h such that h<min{(|xx0|/2),(|xxM|/2)}, then the error ej satisfies the following estimate: (61) |ej|C~max1iM1|eij|+Q~2h4,j1,(61) where Q~2>0, C~1 are constants, and which are independent of h.

Proof.

By the given condition h<min{(|xx0|/2),(|xxM|/2)}, assume that x[xs1,xs+2], we note that s satisfies 2sM3, then ej satisfies the following form by (Equation22)–(Equation24) and (Equation42) (62) ej=u(x,tj)u~(x,tj)=k=s1s+2u(xk,tj)u~(xk,tj)L3,k(x)+O(h4)=k=s1s+2ekjL3,k(x)+O(h4),(62) further, we get the following estimation by the formula (Equation62): |ej|k=s1s+2L3,k(x)|ekj|+Q~2h4C~max1iM1|eij|+Q~2h4, where j1, Q~2>0 and C~=k=s1s+2|L3,k(x)|, according to the property of Lagrange interpolation basis functions, we note that C~1, and Q~2 and C~ are independent of h.

The range of stability for the formula (Equation21) (1iM1) is 0<ρ1/2 in [Citation10]. In conclusion, the convergence estimate of u(x,t) for 3PNDS is stated as follows.

Theorem 3.2

Suppose that u(x,t)C4,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, let 0<ττ0, 0<ρ1/2 and 0<KC~1, then there exists a constant Z~1>0, which is independent of τ and h, such that (63) max0iM,1jN|u(xi,tj)u~(xi,tj)|Z~1(τ+h2).(63)

Proof.

The proof of Theorem 3.2 is developed through the following two cases.

Case 1. (a) If j=1, 1iM1, there exists a constant τ0 such that 0<ττ0, combining with Equation (Equation39) and Lemma 3.4, we have (64) max1iM1|ei1|τmax1iM1|ϵ~i0|τ0Q~1(τ+h2).(64) (b) In addition, ek1 (k=0orM) in Equation (Equation39) satisfies the following estimation by Equation (Equation64) and Lemmas 3.4–3.5 (65) |ek1|=|Gk1e1|KC~max1iM1|ei1|+Q~2h4Kτ0C~Q~1+Q~2(τ+h2),(65) where k=0 or M. Thus, combining with Equations (Equation64) and (Equation65), we have (66) max0iM|ei1|Q~3(τ+h2),(66) where Q~3=max{τ0Q~1,K(τ0C~Q~1+Q~2)}>0 is a constant, and which is independent of τ and h.

Case 2. (a) If 2jN, 1iM1, 0<ρ12 and 0<KC~1, according to the formula (Equation40), Lemmas 3.4 and 3.5, let R0=KQ~2+Q~1, then eij satisfies the following estimation as (67) max1iM1|eij|max0iM|eij1|+τK|ej1|+τQ~1(τ+h2)max0iM|eij1|+τKC~max1iM1|eij1|+Q~2h4+τQ~1(τ+h2)(1+τKC~)max0iM|eij1|+τR0(τ+h2).(67) (b) In addition, if 2jN, k={0,M}, then ekj satisfies the following estimation by Equations (Equation40), (Equation67) and Lemmas 3.4 and 3.5 (68) maxk={0,M}|ekj|maxk={0,M}|Gkjej|KC~max1iM1|eij|+Q~2h4KC~(1+τKC~)max0iM|eij1|+τKC~R0(τ+h2)+KQ~2h4.(68) Since 0<ρ1/2, then there must exist a constant ρ0 such that 0<ρ0ρ, and let R1=R0+KQ~2/ρ0, by the forms (Equation67), (Equation68) and 0<KC~1, we have (69) max0iM|eij|=max|e0j|,|eMj|,max1iM1|eij|(1+τKC~)max0iM|eij1|+τR0(τ+h2)+τKQ~2ρ0h2(1+τKC~)max0iM|eij1|+τR1(τ+h2)(1+τKC~)2max0iM|eij2|+τ(1+τKC~)R1(τ+h2)+τR1(τ+h2)(1+τKC~)jmax0iM|ei0|+(1+τKC~)j1++1τR1(τ+h2)(1+τKC~)j1(1+τKC~)1τR1(τ+h2)1+TKC~NN1R1KC~(τ+h2)R1KC~exp(TKC~)1(τ+h2)=(1+ρ0)KQ~2+ρ0Q~1ρ0KC~exp(TKC~)1(τ+h2)=θ0(τ+h2),(69) where θ0=(((1+ρ0)KQ~2+ρ0Q~1)/ρ0KC~)(exp(TKC~)1)>0 is a constant, and which is independent of τ and h.

Combining with the forms (Equation66) and (Equation69), we have (70) max0iM,1jN|u(xi,tj)u~(xi,tj)|Z~1(τ+h2),(70) where Z~1=max{Q~3,θ0} is a positive constant, and which is independent of τ and h, the proof is completed.

Furthermore, based on Theorem 3.2, the error estimation for the original solution w(x,t) is as follows.

Theorem 3.3

Suppose that u(x,t)C4,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, there exists a positive constant P~1 such that the following error estimation holds (71) max0iM,1jN|w(xi,tj)w~(xi,tj)|P~1(τ+h2),(71) where P~1 is independent of τ and h.

Proof.

According to the condition u(x,t)C4,2([0,1]×[0,T]), we know that |u(x,t)| is bounded for any (x,t)[0,1]×[0,T], and assume that max0x1,0tT|u(x,t)|K2. By Equation (Equation6), it is clear that r(t)>0 for t[0,T], and suppose that r(t)[rmin,rmax] and r~0. By the formula (Equation13), we have (72) w(xi,tj)w~(xi,tj)=u(xi,tj)r(tj)u~(xi,tj)r~(tj)=|u(xi,tj)r~(tj)u(xi,tj)r(tj)+u(xi,tj)r(tj)r(tj)u~(xi,tj)||r(tj)r~(tj)||u(xi,tj)||r~(tj)r(tj)||r(tj)r~(tj)|+|u(xi,tj)u~(xi,tj)||r~(tj)|K2|r~(tj)r(tj)|rmin|r~(tj)|+|u(xi,tj)u~(xi,tj)||r~(tj)|.(72) Next, suppose that |E(t)|[Emin,Emax] for t[0,T], by the forms (Equation12), (Equation26), Lemma 3.5 and Theorem 3.2, we provide the error estimation of |r(tj)r~(tj)| as follows: (73) |r(tj)r~(tj)|=u(x,tj)E(tj)u~(x,tj)E(tj)1Emin|u(x,tj)u~(x,tj)|1EminC~max1iM1|eij|+Q~2h4C~Z~1Emin(τ+h2)+Q~2Eminh4G~(τ+h2),(73) where G~=((C~Z~1+Q~2)/Emin)>0 is a constant. By the above inequality (Equation73), we know that 0<|r~(tj)|G~(τ+h2)+rmax holds, then there must exist a constant η0 such that 0<η0|r~(tj)|, so we have (74) w(xi,tj)w~(xi,tj)K2G~+rminZ~1rminη0(τ+h2)=P~1(τ+h2),(74) where 0iM, 1jN, P~1=((K2G~+rminZ~1)/rminη0)>0, K2 and G~ are constants, which are independent of τ and h, the proof is completed.

In order to obtain the error estimation of p(t), we first present some results as necessary preparations. According to Lemma 3 and Theorem 4 in [Citation41], the following lemmas are provided.

Lemma 3.6

[Citation41]

Let {tˆ1,tˆ2,,tˆm} be any finite set of distinct points in [0,T], mν and {cˆ1,cˆ2,,cˆm} be any real multipliers such that (75) i=1mcˆitˆij=0,j=0,1,,ν1,(75) furthermore, assume gXν, then the following linear functional: (76) L(g)=i=1mcˆig(tˆi)(76) satisfies the conditions (77) L(g)(2(2ν1)!)1Cν(cˆ)1/2g(ν)L2([0,T]),(77) where (78) Cν(cˆ)=(1)ν2(2ν1)!i=1mj=1mcˆicˆjtˆitˆj2ν10.(78)

Similar to the proof in [Citation41], we have the following error estimation.

Lemma 3.7

Suppose that r(t) is an exact function in Xν, and r~α1(t) is an approximation of r(t) in Lemma 3.1, we have the following estimate of rr~α1 as (79) rr~α1L2([0,T])Wˆ0(τ+h2)+Wˆ1(τ+h2)α1+Wˆ2α11/2,(79) where Wˆ0, Wˆ1 and Wˆ2 are positive constants, and independent of τ, h and α1. In addition, if α1=τ+h2, the error estimate of rr~α1 can be obtained by rr~α1L2([0,T])Wˆ0α11/2+Wˆ1+Wˆ2α11/2, in this case, rr~α1L2([0,T])0, as α10.

Proof.

The proof of Lemma 3.7 consists of the following several steps.

Step 1. For any t[0,T], there exists p{0,1,,Nν} such that t[tp,tp+ν]. In Lemma 3.6, replacing g with the error function ζ=rr~α1, and we take cˆj=cj, tˆj=tp+j, j=0,1,,ν, and cˆν+1=1, tˆν+1=t, where ci is a basis function of Lagrange interpolation as follows: (80) ci=j=0,jiνttp+jtp+itp+j,i=0,1,,ν,(80) according to the property of Lagrange interpolation, it can be proved that (81) j=0ν+1cˆjtˆji=0,i=0,1,,ν.(81) In addition, since |ttp+j|ντ and |tp+itp+j|τ obtained by Equation (Equation80), so we have (82) |ci|ντντν=νν,i=0,1,,ν.(82) Step 2. By Lemma 3.6, we have (83) L(ζ)=i=0ν+1cˆiζ(tˆi)=c0ζ(tp)++cνζ(tp+ν)ζ(t),(83) further, the following estimate is developed by (84) |ζ(t)||L(ζ)|+i=0νciζ(tp+i),(84) by Lemma 3.6, we have (85) L(ζ)(2(2ν1)!)1Cν(cˆ)1/2ζ(ν)L2([0,T])=Lζ,(85) where Lζ=(2(2ν1)!)1Cν(cˆ)1/2ζ(ν)L2([0,T]), and L2([0,T]) is abbreviated below as .

Moreover, it follows from result (Equation78) that Cν(cˆ) is equivalent to (86) Cν(cˆ)=(1)ν2(2ν1)!i=0νj=0νcicjtp+itp+j2ν12i=0νcittp+i2ν12(2ν1)!(ντ)2ν1i=0ν|ci|2+2i=0ν|ci|2(2ν1)!ν2ν1(ν+1)νν+12τ2ν1=Aντ2ν1,(86) where Aν=2(2ν1)!ν2ν1[(ν+1)νν+1]2 and Lζ2=(2(2ν1)!)2Aντ2ν1ζ(ν)2. In addition, by the Cauchy–Schwarz inequality, we have (87) i=0νciζ(tp+i)2i=0νci2i=0νζ2(tp+i)(ν+1)ν2νi=0νζ2(tp+i)=Bνi=0νζ2(tp+i),(87) where Bν=(ν+1)ν2ν.

Step 3. According to the formula (Equation84), we get (88) |ζ(t)|2|L(ζ)|+Bνi=0νζ2(tp+i)1/222L2(ζ)+2Bνi=0νζ2(tp+i),(88) thus (89) tptp+νζ2(t)dt2Lζ2+Bνi=0νζ2(tp+i)(tp+νtp),(89) further, by the formulas (Equation77) and (Equation89), we have (90) 0Tζ2(t)dt((2ν1)!)2TAντ2ν1ζ(ν)2+4νBντi=0Nζ2(ti).(90) Since α1r~α1(ν)2Υν(r~α1)Υν(r) for the exact function rXν, and Υν(r) satisfies the following estimate by (Equation73): (91) Υν(r)G~2(τ+h2)2+α1r(ν)2,(91) thus we can obtain (92) r~α1(ν)2r(ν)2+G~2α1(τ+h2)2,(92) further (93) ζ(ν)2=r(ν)r~α1(ν)22r(ν)2+r~α1(ν)24r(ν)2+2G~2α1(τ+h2)2.(93) On the other hand, we get (94) τi=0Nζ2(ti)2TN+1i=0Nζ2(ti)2TN+1i=0Nr(ti)r~(ti)+r~(ti)r~α1(ti)24TN+1i=0Nr(ti)r~(ti)2+r~(ti)r~α1(ti)24TG~2(τ+h2)2+Υν(r~α1)4TG~2(τ+h2)2+Υν(r)8TG~2(τ+h2)2+4Tα1r(ν)2.(94) Since ν2 in (Equation30), combining with (Equation90), (Equation92)–(Equation94), we have (95) rr~α1L2([0,T])2(2ν1)!2TAντ2ν1ζ(ν)2+16νBνT2G~2(τ+h2)2+α1r(ν)2W~0τ2+W~1(τ+h2)2+W~2(τ+h2)2α1+W~3α1,(95) where W~0=4(2ν1)!2T2(ν1)Aνr(ν)2,W~1=32νBνTG~2,W~2=2(2ν1)!2T2νAνG~2,W~3=16νBνTr(ν)2. Furthermore, by (Equation95), we can obtain (96) rr~α1L2([0,T])W~01/2τ+W~11/2(τ+h2)+W~21/2(τ+h2)α1+W~31/2α11/2W~01/2+W~11/2(τ+h2)+W~21/2(τ+h2)α1+W~31/2α11/2=Wˆ0(τ+h2)+Wˆ1(τ+h2)α1+Wˆ2α11/2,(96) where Wˆ0=W~01/2+W~11/2, Wˆ1=W~21/2 and Wˆ2=W~31/2 are positive constants, and independent of τ, h and α1. Specially, if α1=τ+h2, the error estimate of rr~α1 can be obtained by rr~α1L2([0,T])Wˆ0α11/2+Wˆ1+Wˆ2α11/2, in this case, rr~α1L2([0,T])0, as α10, the proof is completed.

Lemma 3.8

[Citation47]

Let <a<b<, and let 0<ε0<. There exists a finite constant K~ such that for 0<εε0, and for every function f which is mth continuously differentiable on the open interval (a,b), we have (97) ab|f(t)|2dtK~εab|f(m)(t)|2dt+K~ε1/(m1)ab|f(t)|2dt.(97)

By Lemmas 3.6–3.8, we obtain an error estimation of p(t)p~α1(t) in the sense of L2-norm as follows.

Theorem 3.4

There exist positive constants {V~i}i=06 such that the following inequality holds: pp~α1L2([0,T])V~0(τ+h2)(ν1)/ν+V~1(τ+h2)+V~2(τ+h2)α12ν+V~3τ+h2(2ν1)/να1+V~4τ+h2α1(ν1)/ν+V~5τ+h2α1(2ν1)/ν+V~6α1(ν1)/2ν, where {V~i}i=06 are independent of τ, h and α1. In addition, if α1=τ+h2, the error estimate of pp~α1 can be obtained by pp~α1L2([0,T])V~0α1(ν1)/2ν+V~1α1(ν+1)/2ν+V~3α1(2ν1)/2ν+(V~2+V~5)α11/2+V~4+V~6α1(ν1)/2ν where ν2, in this case, pp~α1L2([0,T])0, as α10.

Proof.

In Lemma 3.8, we take a=0, b=T and m=ν, respectively. Suppose that r~C[0,T] (r~0) and |r|r¯max on [0,T], we know that |r~| is bounded, then there must exist a constant η1 such that 0<η1|r~|, and the assumption r[rmin,rmax] in Theorem 3.3. Let ζ=rr~α1, by Lemma 3.7, for selecting a set of suitable small τ, h and α1 such that (98) ζL2([0,T])Wˆ0(τ+h2)+Wˆ1(τ+h2)α1+Wˆ2α11/21.(98) In the sense of L2-norm, the error function p(t)p~α1(t) satisfies the following estimate: (99) pp~α1L2([0,T])2=0Trr+r~α1r~2dt1rminη120Trr~α1rr~2dt=1rminη120Trr~α1rr+rrrr~2dt(rmax+r¯max)rminη12rmaxrr~α12+r¯maxrr~2,(99) further, by (Equation99), we can get (100) pp~α1L2([0,T])C¯1rr~α1+C¯2rr~,(100) where C¯1=([rmax(rmax+r¯max)]1/2)/rminη1, C¯2=([r¯max(rmax+r¯max)]1/2)/rminη1 are two positive constants.

If maxtitti+1{|r(t)r~(t)|}max{|r(ti)r~(ti)|,|r(ti+1)r~(ti+1)|}, i=0,1,,N1 holds, by formula (Equation73), we have (101) rr~2=0T(rr~)2dt=i=0N1titi+1(rr~)2dt2TG~2(τ2+h4).(101) By Lemma 3.8, assume that ε=ζL2([0,T])2(ν1)/ν, and the formula (Equation93), the following inequality yields (102) rr~α12=0Trr~α12dtK~ζ(ν)L2([0,T])2+1ζL2([0,T])2(ν1)/νK~4r(ν)2+2G~2α1(τ+h2)2+1ζL2([0,T])2(ν1)/ν.(102) In addition, we note that (a+b)θ~aθ~+bθ~ for a,b>0, θ~(0,1). By (Equation98), and 0<(ν1)/ν<1, we have (103) ζ(ν1)/νWˆ0(ν1)/ν(τ+h2)(ν1)/ν+Wˆ1(ν1)/ντ+h2α1(ν1)/ν+Wˆ2(ν1)/να1(ν1)/2ν.(103) Combining formulas (Equation100)–(Equation103) yields pp~α1L2([0,T])K~C¯12r(ν)+2G~α1(τ+h2)+1ζ(ν1)/ν+2TC¯2G~(τ+h2)K~C¯1Wˆ0(ν1)/ν2r(ν)+1(τ+h2)(ν1)/ν+2TC¯2G~(τ+h2)+2K~C¯1G~Wˆ2(ν1)/ν(τ+h2)α12ν+Wˆ0(ν1)/ντ+h2(2ν1)/να1+K~C¯1Wˆ1(ν1)/ν2r(ν)+1τ+h2α1(ν1)/ν+2G~τ+h2α1(2ν1)/ν+K~C¯1Wˆ2(ν1)/ν2r(ν)+1α1(ν1)/2ν=V~0(τ+h2)(ν1)/ν+V~1(τ+h2)+V~2(τ+h2)α12ν+V~3τ+h2(2ν1)/να1+V~4τ+h2α1(ν1)/ν+V~5τ+h2α1(2ν1)/ν+V~6α1(ν1)/2ν, where V~0=K~C¯1Wˆ0(ν1)/ν2r(ν)+1,V~1=2TC¯2G~,V~2=2K~C¯1G~Wˆ2(ν1)/ν,V~3=2K~C¯1G~Wˆ0(ν1)/ν,V~4=K~C¯1Wˆ1(ν1)/ν2r(ν)+1,V~5=2K~C¯1G~Wˆ1(ν1)/ν,V~6=K~C¯1Wˆ2(ν1)/ν2r(ν)+1 are positive constants, and independent of τ, h and α1. Specially, if α1=τ+h2, the error estimate of pp~α1 can be obtained by pp~α1L2([0,T])V~0α1(ν1)/2ν+V~1α1(ν+1)/2ν+V~3α1(2ν1)/2ν+(V~2+V~5)α11/2+V~4+V~6α1(ν1)/2ν, where ν2, in this case, pp~α1L2([0,T])0, as α10, the proof is complete.

3.1.4. The convergence estimates of 3PNDS for E(t) with noise

In practice, since the measured data of E(t), in the forms (Equation5) and (Equation12), always contain some noise, here we can assume that the given functions f(x),g0(t),g1(t) and ϕ(x,t) are exact data, but Eδ(t) is an approximation of the exact function E(t) with |Eδ(t)E(t)|δ for all t[0,T], where δ denotes a small noise error, then the formula (Equation12) is replaced by (104) rδ(t)=u(x,t)Eδ(t).(104) Similarly, combining with the formula (Equation21), we can use Eδ(t) in (Equation104) to get an stable numerical solution u~δ(xi,tj), and further to obtain r~δ(tj). Thus we can use (r~δ)(tj) and r~δ(tj) to compute p~δ(tj) by the following formula: (105) p~δ(tj)=(r~δ)(tj)r~δ(tj),0jN.(105) Since the values r~δ(tj) in (Equation105) include errors induced by computation and the noise of E(t), it is the key that we use the regularization method to approximate the exact r(t). Similar to (Equation30), for the set of {r~jδ}j=1n with noise errors, and r~jδ=(u~δ(x,tj))/(Eδ(tj)), we define the following smoothing functional on Xν: (106) Υν(r)=1nj=1nr(tj)r~jδ2+α1δr(ν)L2(R)2,(106) where α1δ is a regularization parameter, which is dependent on δ. Furthermore, according to the formulas (Equation31)–(Equation34) and (Equation105), we can obtain the regularized solutions r~α1δ(tj) of the exact values r(tj), moreover, using (r~α1δ)(tj) to compute the approximation values p~α1δ(tj)=(((r~α1δ)(tj))/(r~δ(tj))) of the exact p(tj).

Then we have the following estimation between an stable approximation solution u~δ(x,t) and the exact u(x,t).

Theorem 3.5

Assume that u(x,t)C4,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, let |Eδ(t)E(t)|δ for all t[0,T], 0<ττ0, 0<ρ1/2 and 0<KC~1, then there exist positive constants γ1, γ2 such that max0iM,1jN|u(xi,tj)u~δ(xi,tj)|γ1(τ+h2)+γ2δτ, where γ1 and γ2 are independent of τ, h and δ. Specially, when the temporal step τ=δ1/2, 0<ρ0ρ, and h2=τ/ρτ/ρ0, the following error estimation can be obtained by max0iM,1jN|u(xi,tj)u~δ(xi,tj)|κγ1+γ2δ1/2, where κ=1+(1/ρ0) is a positive constant, in this case, |u(xi,tj)u~δ(xi,tj)|0, as δ0 for 0iM, 1jN.

Proof.

Similar to the proof of Theorem 3.2, suppose |g0(t)|g0max, |g1(t)|g1max for all t[0,T], |E(t)|[Emin,Emax], |Eδ(t)|[Eminδ,Emaxδ] (Eδ(t)0), |u~δ(x,t)|Kˆδ and |u~δ(x,t)|u~maxδ, where the four bounds Eminδ, Emaxδ, Kˆδ and u~maxδ are independent of δ. For simplicity, we define eδij=u(xi,tj)u~δ(xi,tj), eδj=u(x,tj)u~δ(x,tj). In order to prove the related result, we should discuss the following two cases, i.e. j=1 and 2jN.

Case 1. If j=1, 1iM1, there exists a constant τ0 such that 0<ττ0, by Lemma 3.4 and Equation (Equation39), we have (107) max1iM1|eδi1|τmax1iM1|ϵ~i0|τ0Q~1(τ+h2).(107) In addition, according to Lemma 3.5, and Equations (Equation61), (Equation107), and the previous assumption maxi=0,M|Gij|K for 1jN(Gij=(gi(tj))/(E(tj)),i=0,M), let CE=max{(g0max/EminδEmin),(g1max/EminδEmin)}>0, thus eδi1 (i=0,M) satisfies the following estimation: (108) |eδ01|=g0(t1)E(t1)u(x,t1)g0(t1)Eδ(t1)u~δ(x,t1)g0(t1)E(t1)u(x,t1)u~δ(x,t1)+g0(t1)u~δ(x,t1)EminEminδEδ(t1)E(t1)K|eδ1|+u~maxδCEδKC~max1iM1|eδi1|+Q~2h4+u~maxδCEδKτ0C~Q~1(τ+h2)+Q~2h4+u~maxδCEδKτ0C~Q~1+Q~2(τ+h2)+u~maxδCEδ.(108) Analogously, we have (109) |eδM1|=g1(t1)E(t1)u(x,t1)g1(t1)Eδ(t1)u~δ(x,t1)K|eδ1|+u~maxδCEδKτ0C~Q~1+Q~2(τ+h2)+u~maxδCEδ.(109) Combining with the formulas (Equation107), (Equation108) and (Equation109), we have (110) max0iM|eδi1|=max|eδ01|,|eδM1|,max1iM1|eδi1|γ~1(τ+h2)+γ~2δ,(110) where γ~1=max{τ0Q~1,K(τ0C~Q~1+Q~2)} and γ~2=u~maxδCE are positive constants.

Case 2. (a) If 2jN, 1iM1 and 0<ρ12, similar to (Equation67), then eδij satisfies the following estimation as (111) max1iM1|eδij|max0iM|eδij1|+τ|ϕ(xi,tj1)|k=s1s+2u(xk,tj1)E(tj1)u~δ(xk,tj1)Eδ(tj1)|L3,k(x)|+τQ~1(τ+h2)max0iM|eδij1|+τKEmaxEminEminδk=s1s+2|eδkj1|Emaxδ+δKˆδ|L3,k(x)|+τQ~1(τ+h2)max0iM|eδij1|+τKC~EmaxEmaxδEminEminδmax0iM|eδij1|+τKC~EmaxKˆδδEminEminδ+τQ~1(τ+h2)(1+τCP)max0iM|eδij1|+τQ~1(τ+h2)+τCQδ,(111) where CP=KC~EmaxEmaxδ/EminEminδ, CQ=KC~EmaxKˆδ/EminEminδ are positive constants, which are independent of τ, h and δ.

(b) If 2jN, k={0,M}, by 0<KC~1 and the formulas (Equation40), (Equation108), (Equation109), (Equation111), and Lemmas 3.4, 3.5, then eδkj satisfies the following estimation as (112) maxk={0,M}|eδkj|maxi={0,1}gi(tj)E(tj)u(x,tj)gi(tj)Eδ(tj)u~δ(x,tj)max1iM1|eδij|+KQ~2h4+u~maxδCEδ(1+τCP)max0iM|eδij1|+τCQδ+τQ~1(τ+h2)+KQ~2h4+u~maxδCEδ.(112) In the above (Equation112), based on the condition 0<ρ0ρ=τ/h2, thus h2τ/ρ0. Let R2=Q~1+KQ~2/ρ0, R3=τ0CQ+u~maxδCE, by (Equation111), (Equation112) and |eδi0|=0, we get (113) max0iM|eδij|=max|eδ0j|,|eδMj|,max1iM1|eδij|(1+τCP)max0iM|eδij1|+τCQδ+τQ~1(τ+h2)+τKQ~2ρ0h2+u~maxδCEδ(1+τCP)max0iM|eδij1|+τR2(τ+h2)+R3δ(1+τCP)2max0iM|eδij2|+τR2(1+τCP)+1(τ+h2)+R3(1+τCP)+1δτR2k=0j1(1+τCP)k(τ+h2)+R3k=0j1(1+τCP)kδR2CPexp(TCP)1(τ+h2)+R3CPexp(TCP)1δτ=γˆ1(τ+h2)+γˆ2δτ,(113) where γˆ1=(R2/CP)(exp(TCP)1) and γˆ2=(R3/CP)(exp(TCP)1) are positive constants.

Combining with the forms (Equation110) and (Equation113), we have the following estimate of eδij as follows: (114) max0iM,1jN|eδij|γ1(τ+h2)+γ2δτ,(114) where γ1=max{γ~1,γˆ1}, γ2=max{τ0γ~2,γˆ2} are two positive constants, and independent of τ, h and δ.

In addition, when the temporal step τ=δ1/2, h2τ/ρ0, by (Equation114), we can obtain the estimate of eδij as follows: (115) max0iM,1jN|eδij|γ1(τ+τρ0)+γ2δτκγ1+γ2δ1/2,(115) where κ=1+(1/ρ0) is a positive constant, in this case, |eδij|0, as δ0 for 0iM, 1jN, the proof is completed.

Next, we present the error estimation between the original solution w(x,t) and the approximation w~δ(x,t) as follows.

Theorem 3.6

Suppose that u(x,t)C4,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, there exist positive constants {ω~i}i=02 such that the following error estimation holds: max0iM,1jN|w(xi,tj)w~δ(xi,tj)|ω~0(τ+h2)+ω~1δτ+ω~2δ, where {ω~i}i=02 are independent of τ, h and δ. In addition, when the temporal step τ=δ1/2, the error estimate of |w(xi,tj)w~δ(xi,tj)| can be obtained by w(xi,tj)w~δ(xi,tj)ω~0(τ+τρ0)+ω~1δ1/2+ω~2δ=ω~2δ1/2+κω~0+ω~1δ1/2, in this case, |w(xi,tj)w~δ(xi,tj)|0, as δ0 for 0iM, 1jN.

Proof.

The proof can be similarly proved as Theorem 3.3. According to some assumptions in Theorem 3.3, we know that r(t)[rmin,rmax]. In addition, suppose that |r~δ|rmaxδ, and r~δ0 for all t[0,T], where rmaxδ is independent of δ. Then there must exist a constant ηˆ0 such that 0<ηˆ0|r~δ| for t[0,T].

In addition, by Lemma 3.5, we have (116) |eδj|C~max1iM1|eδij|+Q~2h4,(116) furthermore, according to (Equation114) and (Equation116), the error r(tj)r~δ(tj) for 1jN satisfies the following estimate: (117) |r(tj)r~δ(tj)|=u(x,tj)E(tj)u~δ(x,tj)Eδ(tj)1Emin|eδj|+u~maxδEminEminδδ1EminC~max0iM|eδij|+Q~2h4+u~maxδEminEminδδC~γ1Emin(τ+h2)+Q2~Eminh4+C~γ2Eminδτ+u~maxδEminEminδδ(C~γ1+Q2~)Emin(τ+h2)+C~γ2Eminδτ+u~maxδEminEminδδ=γ~3(τ+h2)+γ~4δτ+γ~5δ,(117) where γ~3=(C~γ1+Q2~)/Emin, γ~4=C~γ2/Emin and γ~5=u~maxδ/EminEminδ are positive constants, which are independent of τ, h and δ. For simplicity, here we denote ωr(τ,h,δ):=γ~3(τ+h2)+γ~4δτ+γ~5δ in this paper.

Combining with the formulas (Equation13), (Equation114) and (Equation117), we get (118) w(xi,tj)w~δ(xi,tj)=u(xi,tj)r(tj)u~δ(xi,tj)r~δ(tj)K2|r(tj)r~δ(tj)|rmin|r~δ(tj)|+|u(xi,tj)u~δ(xi,tj)||r~δ(tj)|K2γ~3+rminγ1rminηˆ0(τ+h2)+K2γ~4+rminγ2rminηˆ0δτ+K2γ~5rminηˆ0δ=ω~0(τ+h2)+ω~1δτ+ω~2δ,(118) where ω~0=(K2γ~3+rminγ1)/rminηˆ0, ω~1=(K2γ~4+rminγ2)/rminηˆ0, and ω~2=K2γ~5/rminηˆ0 are positive constants, which are independent of τ, h and δ.

Specially, when the temporal step τ=δ1/2, the error estimate of |w(xi,tj)w~δ(xi,tj)| can be obtained by w(xi,tj)w~δ(xi,tj)ω~0(τ+τρ0)+ω~1δ1/2+ω~2δ=ω~2δ1/2+κω~0+ω~1δ1/2, in this case, |w(xi,tj)w~δ(xi,tj)|0, as δ0 for 0iM, 1jN, the proof is completed.

Furthermore, we have the following error estimation between the Tikhonov regularized solution r~α1δ(t) and the exact r(t).

Theorem 3.7

Suppose that r(t) is an exact function in Xν, and r~α1δ(t) is the Tikhonov regularized solution to r(t) in (Equation106), in the sense of L2-norm, we have the following estimate of rr~α1δ as rr~α1δL2([0,T])Y~0(τ+h2)+Y~1(τ+h2)α1δ+Y~2α1δ+Y~3δτα1δ+Y~4δα1δ+Y~5δτ+Y~6δ, where {Y~i}i=06 are positive constants, and which are independent of τ, h, δ and α1δ.

In addition, if taking the case of τ=δ in Theorem 3.5, where 0<ρ0ρ, and α1δ=δ, the following error estimation can be obtained by rr~α1δL2([0,T])Pˆ0δ1/4+Pˆ4δ1/2+Pˆ1δ3/4+Pˆ2+Pˆ3δ1/4, in this case, rr~α1δL2([0,T])0, as δ0, where {Pˆi}i=04 are positive constants, which are independent of τ, h, δ and α1δ.

Proof.

The proof can be similarly proved as Lemma 3.7. Suppose that ζ=rr~α1δ, by the forms (Equation80)–(Equation89) and Lemma 3.6, we have the following results: (119) 0Tζ2(t)dt((2ν1)!)2TAντ2ν1ζ(ν)2+4νBντi=0Nζ2(ti),(119) where Aν=2(2ν1)!ν2ν1[(ν+1)νν+1]2 and Bν=(ν+1)ν2ν.

Since α1δ(r~α1δ)(ν)2Υν(r~α1δ)Υν(r) for the exact function rXν, where r~α1δ is a regularized solution of r, by (Equation106) and (Equation117), Υν(r) satisfies (120) Υν(r)ωr2(τ,h,δ)+α1δr(ν)2,(120) thus we can obtain (121) r~α1δ(ν)2r(ν)2+1α1δωr2(τ,h,δ),(121) further, by (Equation121), we get (122) ζ(ν)2=r(ν)r~α1δ(ν)22r(ν)2+r~α1δ(ν)24r(ν)2+2α1δωr2(τ,h,δ).(122) On the other hand, similar to (Equation94), we have (123) τi=0Nζ2(ti)2TN+1i=0Nζ2(ti)2TN+1i=0Nr(ti)r~δ(ti)+r~δ(ti)r~α1δ(ti)24TN+1i=0Nr(ti)r~δ(ti)2+r~δ(ti)r~α1δ(ti)24Tωr2(τ,h,δ)+Υν(r~α1δ)4Tωr2(τ,h,δ)+Υν(r)8Tωr2(τ,h,δ)+4Tα1δr(ν)2.(123) Thus, combining with (Equation119), (Equation122) and (Equation123), we can get (124) rr~α1δL2([0,T])2(2ν1)!2TAντ2ν1ζ(ν)2+16νBνT2ωr2(τ,h,δ)+α1δr(ν)2W~5α1δ+W~6γ~3(τ+h2)+γ~4δτ+γ~5δ2+W~4τ2+W~7α1δ,(124) where W~4=4(2ν1)!2AνT(2ν2)r(ν)2,W~5=2(2ν1)!2AνT2ν,W~6=32νBνT,W~7=16νBνTr(ν)2, are positive constants, which are independent of τ, h, δ and α1δ. Furthermore, by (Equation124), we have (125) rr~α1δL2([0,T])W~41/2+W~61/2γ~3(τ+h2)+W~51/2γ~3(τ+h2)α1δ+W~7α1δ1/2+W~51/2γ~4δτα1δ+W~51/2γ~5δα1δ+W~61/2γ~4δτ+W~61/2γ~5δ=Y~0(τ+h2)+Y~1(τ+h2)α1δ+Y~2α1δ+Y~3δτα1δ+Y~4δα1δ+Y~5δτ+Y~6δ,(125) where Y~0=W~41/2+W~61/2γ~3, Y~1=W~51/2γ~3, Y~2=W~71/2, Y~3=W~51/2γ~4, Y~4=W~51/2γ~5, Y~5=W~61/2γ~4 and Y~6=W~61/2γ~5, which are positive constants, and independent of τ, h, δ and α1δ.

In addition, if we take the case of τ=δ, 0<ρ0ρ in Theorem 3.5, and α1δ=δ for discussing the formula (Equation125), thus (126) rr~α1δL2([0,T])Pˆ0δ+Pˆ1δ+Pˆ2α1δ+Pˆ3δα1δ+Pˆ4δα1δPˆ0δ1/4+Pˆ4δ1/2+Pˆ1δ3/4+Pˆ2+Pˆ3δ1/4,(126) where Pˆ0=κW~41/2+W~61/2γ~3+W~61/2γ~4,Pˆ1=W~61/2γ~5,Pˆ2=W~71/2,Pˆ3=κW~51/2γ~3+W~51/2γ~4,Pˆ4=W~51/2γ~5, which are independent of τ, h, δ and α1δ. If τ=δ and α1δ=δ, then rr~α1δL2([0,T])0, as δ0, the proof is complete.

Furthermore, we have the convergence estimate between the approximation solution p~α1δ(t) and the exact p(t).

Theorem 3.8

In the sense of L2-norm, we have the following estimate of pp~α1δ as pp~α1δL2([0,T])α1δ(ν1)/2νD2+D19δτ+D26δ+τ+h2(2ν1)/νD10+D111α1δ(ν1)/2ν+δτ(ν1)/νD5+D22δτ+D29δ+δ(2ν1)/νD30+D231τ+τ+h2(ν1)/νD0+D24δ+D17δτ+1α1δ(ν1)/νD1+D25δ+D18δτ+τ+h2D7+D12α1δ(ν1)/2ν+D14δα1δ(ν1)/ν+D15δτ(ν1)/ν+D16δ(ν1)/ν+D13δτα1δ(ν1)/ν+δD9+D81τ+D6δ(ν1)/ν+δα1δ(ν1)/νD4+D21δτ+D28δ+1τ(ν1)/νD3+D20δτ+D27δ, where {Di}i=030 are positive constants, and which are independent of τ, h, δ and α1δ. Specially, if taking the case of τ=δ, h2τ/ρ0 and α1δ=δ, the following error estimation can be obtained by pp~α1δL2([0,T])D2+D19δ1/2+D26δ+κ(2ν1)/νD10δ(3ν1)/4ν+D11δ1/2+D5+D22δ1/2+D29δδ(ν1)/4ν+D30δ(7ν3)/4ν+D23δ(5ν3)/4ν+κ(ν1)/νD1+D25δ+D18δ1/2+D0+D24δ+D17δ1/2δ(ν1)/4ν+κD7+D12δ(ν1)/4ν+D14δ(3ν3)/4ν+D15δ(ν1)/2ν+D16δ(ν1)/ν+D13δ(ν1)/4νδ(ν+1)/4ν+D9δ(3ν+1)/4ν+D8δ(ν+1)/4ν+D6δ(3ν3)/4ν+D4+D21δ1/2+D28δδ(ν1)/2ν+D3+D20δ1/2+D27δδ(ν1)/4ν, where ν2. In this case, pp~α1δL2([0,T])0, as δ0.

Proof.

The proof can be similarly proved as Theorem 3.4. Suppose that r~δC[0,T] (r~δ0) and |r|r¯max on [0,T], we know that |r~δ| is bounded, then there must exist a constant η1δ such that 0<η1δ|r~δ|, and the assumption r[rmin,rmax] in Theorem 3.3. Let ζ=rr~α1δ, and assume that ε=ζL2([0,T])2(ν1)/ν. In the sense of L2-norm, by the forms (Equation14) and (Equation105), the error function p(t)p~α1δ(t) satisfies the following estimate: (127) pp~α1δL2([0,T])2=0Trr+r~α1δr~δ2dt1rminη1δ20Trr~α1δrr~δ2dt=1rminη1δ20Trr~α1δrr+rrrr~δ2dt(rmax+r¯max)rminη1δ2rmaxrr~α1δ2+r¯maxrr~δ2,(127) furthermore, by (Equation127), we can get (128) pp~α1δL2([0,T])Cˆ1rr~α1δ+Cˆ2rr~δ,(128) where Cˆ1=([rmax(rmax+r¯max)]1/2)/rminη1δ, Cˆ2=([r¯max(rmax+r¯max)]1/2)/rminη1δ are two positive constants, and independent of τ, h, δ and α1δ.

If maxtitti+1{|r(t)r~δ(t)|}max{|r(ti)r~δ(ti)|,|r(ti+1)r~δ(ti+1)|}, i=0,1,,N1 holds, by formula (Equation117), we have (129) rr~δ2=0T(rr~δ)2dt=i=0N1titi+1(rr~δ)2dtTωr2(τ,h,δ).(129) By Lemma 3.8 and (Equation122), the following inequality satisfies (130) rr~α1δ2=0Trr~α1δ2dtK~ζ(ν)L2([0,T])2+1ζL2([0,T])2(ν1)/νK~4r(ν)2+2α1δωr2(τ,h,δ)+1ζL2([0,T])2(ν1)/ν.(130) In addition, by Theorem 3.7, and 0<(ν1)/ν<1, we have (131) ζ(ν1)/νY~0(ν1)/ν(τ+h2)(ν1)/ν+Y~1(τ+h2)α1δ(ν1)/ν+Y~2(ν1)/να1δ(ν1)/2ν+Y~3(ν1)/νδτα1δ(ν1)/ν+Y~4δα1δ(ν1)/ν+Y~5δτ(ν1)/ν+Y~6(ν1)/νδ(ν1)/ν.(131) Combining with the formulas (Equation127)–(Equation131), we have pp~α1δL2([0,T])K~Cˆ12r(ν)+2α1δωr(τ,h,δ)+1ζ(ν1)/ν+TCˆ2ωr(τ,h,δ)=K~Cˆ12r(ν)+1ζ(ν1)/ν+2K~Cˆ1γ~3(τ+h2)+γ~4δτ+γ~5δα1δζ(ν1)/ν+TCˆ2γ~3(τ+h2)+γ~4δτ+γ~5δα1δ(ν1)/2νD2+D19δτ+D26δ+τ+h2(2ν1)/νD10+D111α1δ(ν1)/2ν+δτ(ν1)/νD5+D22δτ+D29δ+δ(2ν1)/νD30+D231τ+τ+h2(ν1)/νD0+D24δ+D17δτ+1α1δ(ν1)/νD1+D25δ+D18δτ+τ+h2D7+D12α1δ(ν1)/2ν+D14δα1δ(ν1)/ν+D15δτ(ν1)/ν+D16δ(ν1)/ν+D13δτα1δ(ν1)/ν+δD9+D81τ+D6δ(ν1)/ν+δα1δ(ν1)/νD4+D21δτ+D28δ+1τ(ν1)/νD3+D20δτ+D27δ, where D0=K~Cˆ1(2r(ν)+1)Y~0(ν1)/ν, D1=K~Cˆ1(2r(ν)+1)Y~1(ν1)/ν, D2=K~Cˆ1(2r(ν)+1)Y~2(ν1)/ν, D3=K~Cˆ1(2r(ν)+1)Y~3(ν1)/ν, D4=K~Cˆ1(2r(ν)+1)Y~4(ν1)/ν, D5=K~Cˆ1(2r(ν)+1)Y~5(ν1)/ν, D6=K~Cˆ1(2r(ν)+1)Y~6(ν1)/ν, D7=TCˆ2γ~3, D8=TCˆ2γ~4, D9=TCˆ2γ~5, D10=2K~/α1δCˆ1γ~3Y~0(ν1)/ν, D11=2K~/α1δCˆ1γ~3Y~1(ν1)/ν, D12=2K~/α1δCˆ1γ~3Y~2(ν1)/ν, D13=2K~/α1δCˆ1γ~3Y~3(ν1)/ν, D14=2K~/α1δCˆ1γ~3Y~4(ν1)/ν, D15=2K~/α1δCˆ1γ~3Y~5(ν1)/ν, D16=2K~/α1δCˆ1γ~3Y~6(ν1)/ν, D17=2K~/α1δCˆ1γ~4Y~0(ν1)/ν, D18=2K~/α1δCˆ1γ~4Y~1(ν1)/ν, D19=2K~/α1δCˆ1γ~4Y~2(ν1)/ν, D20=2K~/α1δCˆ1γ~4Y~3(ν1)/ν, D21=2K~/α1δCˆ1γ~4Y~4(ν1)/ν, D22=2K~/α1δCˆ1γ~4Y~5(ν1)/ν, D23=2K~/α1δCˆ1γ~4Y~6(ν1)/ν, D24=2K~/α1δCˆ1γ~5Y~0(ν1)/ν, D25=2K~/α1δCˆ1γ~5Y~1(ν1)/ν, D26=2K~/α1δCˆ1γ~5Y~2(ν1)/ν, D27=2K~/α1δCˆ1γ~5Y~3(ν1)/ν, D28=2K~/α1δCˆ1γ~5Y~4(ν1)/ν, D29=2K~/α1δCˆ1γ~5Y~5(ν1)/ν and D30=2K~/α1δCˆ1γ~5Y~6(ν1)/ν are positive constants, which are independent of τ, h, δ and α1δ.

In addition, if we take the case of τ=δ, h2τ/ρ0 and α1δ=δ, thus we have pp~α1δL2([0,T])D2+D19δ1/2+D26δ+κ(2ν1)/νD10δ(3ν1)/4ν+D11δ1/2+D5+D22δ1/2+D29δδ(ν1)/4ν+D30δ(7ν3)/4ν+D23δ(5ν3)/4ν+κ(ν1)/νD1+D25δ+D18δ1/2+D0+D24δ+D17δ1/2δ(ν1)/4ν+κD7+D12δ(ν1)/4ν+D14δ(3ν3)/4ν+D15δ(ν1)/2ν+D16δ(ν1)/ν+D13δ(ν1)/4νδ(ν+1)/4ν+D9δ(3ν+1)/4ν+D8δ(ν+1)/4ν+D6δ(3ν3)/4ν+D4+D21δ1/2+D28δδ(ν1)/2ν+D3+D20δ1/2+D27δδ(ν1)/4ν, where ν2. If τ=δ and α1δ=δ, then pp~α1δL2([0,T])0, as δ0, the proof is completed.

3.2. Five-point numerical difference scheme

Let k=4,j=1,2,,N, in Equations (Equation15) and (Equation16), we choose an appropriate coefficient set {a0,a1,a2,a3,a4}, then we can obtain the derived set {b0,b1,b2,b3,b4} and the nodes-subscript set {i0,i1,i2,i3,i4} in Table , so we can get general difference equation scheme (132) u~(xi,tj)=u~(xi,tj1)+ρ{b0u~(xi0,tj1)+b1u~(xi1,tj1)+b2u~(xi2,tj1)+b3u~(xi3,tj1)+b4u~(xi4,tj1)}+τr~(tj1)ϕ(xi,tj1).(132)

Table 2. The different configuration parameters in formula (Equation132).

Similar to 3PNDS, where r~(t0)=1 and u~(xi,0)=f(xi),0iM, 1jN. We use quartic-Lagrange polynomial interpolation method to calculate u~(x,tj), assume that x[xi,xi+4], i=0,1,,M4 (but xx0 and xM), then we have the following formula: (133) u~(x,tj)=k=ii+4u~(xk,tj)L4,k(x),(133) where L4,k(x) denotes a quartic-Lagrange basis function, and (134) L4,k(x)=ω5(x)(xxk)ω5(xk),(134) where ω5(x)=k=ii+4(xxk), xk are nodes in the spatial domain. Similar to the implementation for the formulas (Equation26)–(Equation35) of 3PNDS, we can obtain the approximation values r~(tj),w~(xi,tj),r~(tj) and p~(tj) to the values r(tj),w(xi,tj),r(tj) and p(tj), respectively. For convenience, the above approach is called 5PNDS in this paper. Similar to 3PNDS, we have the following truncation error estimation for 5PNDS.

Theorem 3.9

Let u(x,t)C5,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,1], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] (E(t)0), then the truncation error order of 5PNDS can reach O(τ+h3) at least.

The convergence estimate of u(x,t) for 5PNDS is presented as follows.

Theorem 3.10

Suppose that u(x,t)C5,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, let 0<ττ0, 0<ρ3/8 and 0<KC~1, then there exists a constant Z~2>0, which is independent of τ and h, such that max0iM,0jN|u(xi,tj)u~(xi,tj)|Z~2(τ+h3).

Proof.

The proof can be similarly proved as Theorem 3.2.

The error estimation of the original solution w(x,t) for 5PNDS is stated as follows.

Theorem 3.11

There exists a positive constant P~2 such that the following error estimation holds: %max0iM,1jN|w(xi,tj)w~(xi,tj)|P~2(τ+h3), where P~2 is independent of τ and h.

Similarly, we obtain an error estimation of p(t)p~α2(t) for 5PNDS in the sense of L2-norm as follows.

Theorem 3.12

There exist positive constants {Vˆi}i=06 such that the following inequality holds: pp~α2L2([0,T])Vˆ0(τ+h3)(ν1)/ν+Vˆ1(τ+h3)+Vˆ2(τ+h3)α22ν+Vˆ3τ+h3(2ν1)/να2+Vˆ4τ+h3α2(ν1)/ν+Vˆ5τ+h3α2(2ν1)/ν+Vˆ6α2(ν1)/2ν, where {Vˆi}i=06 are independent of τ, h and α2. In addition, if α2=τ+h3, the error estimate of pp~α2 can be obtained by pp~α2L2([0,T])V~0α2(ν1)/2ν+V~1α2(ν+1)/2ν+V~3α2(2ν1)/2ν+(V~2+V~5)α21/2+V~4+V~6α2(ν1)/2ν, where ν2, in this case, pp~α2L2([0,T])0, as α20.

3.3. Including the case of E(t) with noise for 5PNDS

Similarly, we have the following estimation between an stable approximation solution u~δ(x,t) and the exact u(x,t).

Theorem 3.13

Assume that u(x,t)C5,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, let |Eδ(t)E(t)|δ for all t[0,T], 0<ττ0, 0<ρ3/8 and 0<KC~1, then there exist positive constants γ¯1, γ¯2 such that max0iM,1jN|u(xi,tj)u~δ(xi,tj)|γ¯1(τ+h3)+γ¯2δτ, where γ¯1 and γ¯2 are independent of τ, h and δ. Specially, when the temporal step τ=δ1/2, 0<ρ0ρ and h2=τ/ρτ/ρ0, the following error estimation can be obtained by max0iM,1jN|u(xi,tj)u~δ(xi,tj)|κγ¯1+γ¯2δ1/2, in this case, |u(xi,tj)u~δ(xi,tj)|0, as δ0 for 0iM, 1jN.

Next, we present the error estimation between the original solution w(x,t) and the approximation w~δ(x,t) as follows.

Theorem 3.14

Suppose that u(x,t)C5,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, there exist positive constants {ω~¯i}i=02 such that the following error estimation holds: %max0iM,1jN|w(xi,tj)w~δ(xi,tj)|ω~¯0(τ+h3)+ω~¯1δτ+ω~¯2δ, where {ω~¯i}i=02 are independent of τ, h and δ. In addition, when the temporal step τ=δ1/2, the error estimate of |w(xi,tj)w~δ(xi,tj)| can be obtained by w(xi,tj)w~δ(xi,tj)ω~¯0(τ+τρ0)+ω~¯1δ1/2+ω~¯2δ=ω~¯2δ1/2+κω~¯0+ω~¯1δ1/2, in this case, |w(xi,tj)w~δ(xi,tj)|0, as δ0 for 0iM, 1jN.

Finally, we have the convergence estimate between the approximation solution p~α2δ(t) and the exact p(t).

Theorem 3.15

In the sense of L2-norm, we have the following estimate of pp~α2δ as pp~α2δL2([0,T])α2δ(ν1)/2νD¯2+D¯19δτ+D¯26δ+τ+h3(2ν1)/νD¯10+D¯111α2δ(ν1)/2ν+δτ(ν1)/νD¯5+D¯22δτ+D¯29δ+δ(2ν1)/νD¯30+D¯231τ+τ+h3(ν1)/νD¯0+D¯24δ+D¯17δτ+1α2δ(ν1)/νD¯1+D¯25δ+D¯18δτ+τ+h3D¯7+D¯12α2δ(ν1)/2ν+D¯14δα2δ(ν1)/ν+D¯15δτ(ν1)/ν+D¯16δ(ν1)/ν+D¯13δτα2δ(ν1)/ν+δD¯9+D¯81τ+D¯6δ(ν1)/ν+δα2δ(ν1)/νD¯4+D¯21δτ+D¯28δ+1τ(ν1)/νD¯3+D¯20δτ+D¯27δ, where {D¯i}i=030 are positive constants, and which are independent of τ, h, δ and α2δ. Specially, if taking the case of τ=δ, h2τ/ρ0 and α2δ=δ, the following error estimation can be obtained by pp~α2δL2([0,T])D¯2+D¯19δ1/2+D¯26δ+κ(2ν1)/νD¯10δ(3ν1)/4ν+D¯11δ1/2+D¯5+D¯22δ1/2+D¯29δδ(ν1)/4ν+D¯30δ(7ν3)/4ν+D¯23δ(5ν3)/4ν+κ(ν1)/νD¯1+D¯25δ+D¯18δ1/2+D¯0+D¯24δ+D¯17δ1/2δ(ν1)/4ν+κD¯7+D¯12δ(ν1)/4ν+D¯14δ(3ν3)/4ν+D¯15δ(ν1)/2ν+D¯16δ(ν1)/ν+D¯13δ(ν1)/4νδ(ν+1)/4ν+D¯9δ(3ν+1)/4ν+D¯8δ(ν+1)/4ν+D¯6δ(3ν3)/4ν+D¯4+D¯21δ1/2+D¯28δδ(ν1)/2ν+D¯3+D¯20δ1/2+D¯27δδ(ν1)/4ν, where ν2. In this case, pp~α2δL2([0,T])0, as δ0.

3.4. Seven-point numerical difference scheme

Let k=6,j=1,2,,N, in the formulas (Equation15) and (Equation16), we choose the appropriate coefficient set {a0,a1,a2,a3,a4,a5,a6}, then we obtain the derived set {b0,b1,b2,b3,b4,b5,b6} and the nodes-subscript set {i0,i1,i2,i3,i4,i5,i6} in Table , so we can get the following difference equation scheme: (135) u~(xi,tj)=u~(xi,tj1)+ρ{b0u~(xi0,tj1)+b1u~(xi1,tj1)+b2u~(xi2,tj1)+b3u~(xi3,tj1)+b4u~(xi4,tj1)+b5u~(xi5,tj1)+b6u~(xi6,tj1)}+τr~(tj1)ϕ(xi,tj1).(135)

Table 3. The different configuration parameters in formula (Equation135).

Similar to 5PNDS, we use quartic-Lagrange polynomial interpolation method to approximate u~(x,tj) according to Equation (Equation133), the rest of approximation values r~(tj),w~(xi,tj),r~(tj) and p~(tj) to the values r(tj),w(xi,tj),r(tj) and p(tj) are computed, respectively, by the forms (Equation26)–(Equation35) in the 3PNDS. For convenience, the above approach is called 7PNDS in this paper. Similar to 3PNDS, we have the following truncation error estimation for 7PNDS.

Theorem 3.16

Let u(x,t)C5,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,1], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] (E(t)0), then the truncation error order of 7PNDS can reach O(τ+h5) at least.

The convergence estimate of u(x,t) for 7PNDS is proposed as follows.

Theorem 3.17

Suppose that u(x,t)C5,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, let 0<ττ0, 0<ρ18/49 and 0<KC~1, then there exists a constant Z~3>0, which is independent of τ and h, such that max0iM,0jN|u(xi,tj)u~(xi,tj)|Z~3(τ+h5).

Proof.

The proof can be similarly proved as Theorem 3.2.

The error estimation of the original solution w(x,t) for 7PNDS is proposed as follows.

Theorem 3.18

There exists a positive constant P~3 such that the following error estimation holds: max0iM,1jN|w(xi,tj)w~(xi,tj)|P~3(τ+h5), where P~3 is independent of τ and h.

Analogously, we obtain an error estimation of p(t)p~α3(t) for 7PNDS in the sense of L2-norm as follows.

Theorem 3.19

There exist positive constants {V¯i}i=06 such that the following inequality holds: pp~α3L2([0,T])V¯0(τ+h5)(ν1)/ν+V¯1(τ+h5)+V¯2(τ+h5)α32ν+V¯3τ+h5(2ν1)/να3+V¯4τ+h5α3(ν1)/ν+V¯5τ+h5α3(2ν1)/ν+V¯6α3(ν1)/2ν, where {V¯i}i=06 are independent of τ, h and α3. In addition, if α3=τ+h5, the error estimate of pp~α3 can be obtained by pp~α3L2([0,T])V~0α3(ν1)/2ν+V~1α3(ν+1)/2ν+V~3α3(2ν1)/2ν+(V~2+V~5)α31/2+V~4+V~6α3(ν1)/2ν, where ν2, in this case, pp~α3L2([0,T])0, as α30.

3.5. Including the case of E(t) with noise for 7PNDS

Similarly, we have the following estimation between an stable approximation solution u~δ(x,t) and the exact u(x,t).

Theorem 3.20

Assume that u(x,t)C5,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, let |Eδ(t)E(t)|δ for all t[0,T], 0<ττ0, 0<ρ18/49 and 0<KC~1, then there exist positive constants γ¯ˆ1, γ¯ˆ2 such that max0iM,1jN|u(xi,tj)u~δ(xi,tj)|γ¯ˆ1(τ+h5)+γ¯ˆ2δτ, where γ¯ˆ1 and γ¯ˆ2 are independent of τ, h and δ. Specially, when the temporal step τ=δ1/2, 0<ρ0ρ, and h2=τ/ρτ/ρ0, the following error estimation can be obtained by max0iM,1jN|u(xi,tj)u~δ(xi,tj)|κγ¯ˆ1+γ¯ˆ2δ1/2, in this case, |u(xi,tj)u~δ(xi,tj)|0, as δ0 for 0iM, 1jN.

Next, we present the error estimation between the original solution w(x,t) and the approximation w~δ(x,t) as follows.

Theorem 3.21

Suppose that u(x,t)C5,2([0,1]×[0,T]), ϕ(x,t)C[0,1]×[0,T], f(x)C[0,T], g0(t)C[0,T], g1(t)C[0,T], E(t)C[0,T] and E(t)0, there exist positive constants {ω~ˆi}i=02 such that the following error estimation holds: %max0iM,1jN|w(xi,tj)w~δ(xi,tj)|ω~ˆ0(τ+h5)+ω~ˆ1δτ+ω~ˆ2δ, where {ω~ˆi}i=02 are independent of τ, h and δ. In addition, when the temporal step τ=δ1/2, the error estimate of |w(xi,tj)w~δ(xi,tj)| can be obtained by w(xi,tj)w~δ(xi,tj)ω~ˆ0(τ+τρ0)+ω~ˆ1δ1/2+ω~ˆ2δ=ω~ˆ2δ1/2+κω~ˆ0+ω~ˆ1δ1/2, in this case, |w(xi,tj)w~δ(xi,tj)|0, as δ0 for 0iM, 1jN.

Finally, we have the convergence estimate between the approximation solution p~α3δ(t) and the exact p(t).

Theorem 3.22

In the sense of L2-norm, we have the following estimate of pp~α3δ as pp~α3δL2([0,T])α3δ(ν1)/2νD~2+D~19δτ+D~26δ+τ+h5(2ν1)/νD~10+D~111α3δ(ν1)/2ν+δτ(ν1)/νD~5+D~22δτ+D~29δ+δ(2ν1)/νD~30+D~231τ+τ+h5(ν1)/νD~0+D~24δ+D~17δτ+1α3δ(ν1)/νD~1+D~25δ+D~18δτ+τ+h5D~7+D~12α3δ(ν1)/2ν+D~14δα3δ(ν1)/ν+D~15δτ(ν1)/ν+D~16δ(ν1)/ν+D~13δτα3δ(ν1)/ν+δD~9+D~81τ+D~6δ(ν1)/ν+δα3δ(ν1)/νD~4+D~21δτ+D~28δ+1τ(ν1)/νD~3+D~20δτ+D~27δ, where {D~i}i=030 are positive constants, and which are independent of τ, h, δ and α3δ. Specially, if taking the case of τ=δ, h2τ/ρ0, and α3δ=δ, the following error estimation can be obtained by pp~α3δL2([0,T])D~2+D~19δ1/2+D~26δ+κ(2ν1)/νD~10δ(3ν1)/4ν+D~11δ1/2+D~5+D~22δ1/2+D~29δδ(ν1)/4ν+D~30δ(7ν3)/4ν+D~23δ(5ν3)/4ν+κ(ν1)/νD~1+D~25δ+D~18δ1/2+D~0+D~24δ+D~17δ1/2δ(ν1)/4ν+κD~7+D~12δ(ν1)/4ν+D~14δ(3ν3)/4ν+D~15δ(ν1)/2ν+D~16δ(ν1)/ν+D~13δ(ν1)/4νδ(ν+1)/4ν+D~9δ(3ν+1)/4ν+D~8δ(ν+1)/4ν+D~6δ(3ν3)/4ν+D~4+D~21δ1/2+D~28δδ(ν1)/2ν+D~3+D~20δ1/2+D~27δδ(ν1)/4ν, where ν2. In this case, pp~α3δL2([0,T])0, as δ0.

4. Numerical examples

In order to test the efficiency of our method to solve the one-dimensional inverse parabolic differential equations with a control parameter, we present three different numerical examples. In numerical computation, we consistently choose T=1,h=1/10,1/15, τ=(h/4)2 and ν=3, respectively, in three numerical examples.

In practical applications, the inherent noisy errors are usually contained the measured values, so we should replace the original exact data in the initial and boundary conditions of problem (Equation1)–(Equation4) by the following formula: (136) fδ(x)=f(x)(1+δ(2rand()1)),(136) where δ denotes the noise level, rand() is a pseudo-random value, which can be generated by using MATLAB rand. Similarly, other given data g0δ(t), g1δ(t), ϕδ(x,t) and Eδ(t) are obtained by the formula (Equation136) in the numerical simulation of this paper. The numerical results are developed in the following details, and which are derived from the measured data with four different noise levels. The noise level δ is added to all given functions f(x), g0(t), g1(t), ϕ(x,t) and E(t) in three numerical tests.

In addition, the ratios of the absolute error of w(x,t) with respect to its convergence order at point (h1,τ1) are denoted by w1, w2 and w3 for 3PNDS, 5PNDS and 7PNDS, respectively, as follows: (137) w1=|w(h1,τ1)w~δ(h1,τ1)|τ1+h12,(137) (138) w2=|w(h1,τ1)w~δ(h1,τ1)|τ1+h13,(138) (139) w3=|w(h1,τ1)w~δ(h1,τ1)|τ1+h15,(139) based on the above description, we take the spatial-step h1=1/(2i+1), 1/(5i+1) and 1/(6i+1) for 3PNDS, 5PNDS and 7PNDS, respectively, let τ1=(h1/4)2, where i=1,2,,20.

In the following figures and tables described, E2w and Ew denote the Root-Mean-Square (RMS) error and the maximum error between the exact solution w(xi,tj) and the numerical solution w~δ(xi,tj), respectively, where tj=0.1,0.2,,1. E2w(tj) and Ew(tj) are defined by (140) E2w(tj)=1M+1i=0Mw~δ(xi,tj)w(xi,tj)21/2(140) and (141) Ew(tj)=max0iMw~δ(xi,tj)w(xi,tj).(141)

Similar to Equations (Equation140) and (Equation141), E2p and Ep denote the RMS error and the maximum error between the exact solution p(tj) and the numerical solution p~αiδ(tj), respectively, where j=0,1,2,,N. E2p and Ep are defined by (142) E2p=1N+1j=0Np~αiδ(tj)p(tj)21/2(142) and (143) Ep=max0jNp~αiδ(tj)p(tj).(143)

The regularization parameters {αiδ}13 above are important to the numerical results obtained by using the regularization scheme. As shown in Theorems 3.8, 3.15 and 3.22, when the parameters αiδ=δ1/2, we can get the convergence order δ(ν1)/4ν, ν2. However, since the final error estimations include some constants which cannot be known, such prior information is unlikely to be used, so the appropriate values αiδ are more difficult to directly obtain, and the case of αiδ=δ1/2 is not available in practice. In this paper, we select the parameters αiδ by another heuristic method, which is similar to see literatures [Citation42,Citation45,Citation48], as follows: (144) αiδ=M12M22n~iδ1/2,1n~iN~i,(144) where the matrices M1=(|titj|2ν1)n×n, M2=2(1)ν(2ν1)!ndiag(1,1,,1), 2 denotes the 2-norm of a matrix, n~i and N~i are positive integers. The optimal value n~i of n~i can be obtained by minimizing an error model, more details for the computation of values n~i may be found in [Citation41,Citation42]. In addition, with the external data containing different noise levels, the numerical results of the proposed schemes in this article are compared with the RBF method (hereinafter abbreviated as RBFs) given in [Citation14].

Example 4.1

We consider the first problem (Equation1)–(Equation5): ϕ(x,t)=(π2(t+1)2)exp(t2)(cos(πx)+sin(πx)),f(x)=cos(πx)+sin(πx),g0(t)=exp(t2),g1(t)=exp(t2),x=0.25,E(t)=2exp(t2), and Example 4.1 has the exact solution as follows: w(x,t)=exp(t2)(cos(πx)+sin(πx)),p(t)=1+t2.

In Figures  and , the results present the values of E2w and Ew obtained by 3PNDS, 5PNDS and 7PNDS with respect to different time levels tj, respectively, the measured data are added to four various noise levels, where δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.1. When the noise level δ is smaller, it is shown that the numerical results of approximation solutions w~δ(x,t) get more accurate among the techniques 3PNDS, 5PNDS and 7PNDS for the following two cases: (a) the spatial-step h becomes smaller; (b) the points number k of the k-step discretization schemes (such as 3PNDS (k=2), 5PNDS (k=4) and 7PNDS (k=6)) becomes bigger, but this effect is not obvious as the error level δ is growing. According to Figure  and Table , it is observed that the approximation source p~αiδ(t) becomes more accurate, as the points number k increases. In addition, the numerical results basically keep unchanged to different noise level δ for taking the same numerical approach and spatial-step h, moreover, when δ becomes smaller, the accuracy of numerical results is generally increasing.

Figure 1. The RMS error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.1. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 1. The RMS error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.1. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 2. The maximum error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.1. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 2. The maximum error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.1. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 3. The ratio plots of w1, w2 and w3 for 3PNDS, 5PNDS and 7PNDS, respectively, in Example 4.1, and the numbers ‘1,2,’ indicates the expressed sequence for different h1 above. (a) 3PNDS, (b) 5PNDS and (c) 7PNDS.

Figure 3. The ratio plots of w1, w2 and w3 for 3PNDS, 5PNDS and 7PNDS, respectively, in Example 4.1, and the numbers ‘1,2,…’ indicates the expressed sequence for different h1 above. (a) 3PNDS, (b) 5PNDS and (c) 7PNDS.

Figure 4. The exact function p(t) and its approximation for three difference schemes as 3PNDS, 5PNDS and 7PNDS with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.1. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 4. The exact function p(t) and its approximation for three difference schemes as 3PNDS, 5PNDS and 7PNDS with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.1. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Table 4. The results of E2p and Ep obtained by 3PNDS, 5PNDS, 7PNDS and RBFs, respectively, taking the case of αiδ in (Equation144) for Example 4.1.

Example 4.2

We consider the second problem (Equation1)–(Equation5): ϕ(x,t)=exp(t)(2+t3+(1+t3)cos(x)),f(x)=cos(x)+1,g0(t)=2exp(t),g1(t)=exp(t)(cos1+1),x=π/4,E(t)=exp(t)(2/2+1), and Example 4.2 has the exact solution w(x,t)=exp(t)(cos(x)+1),p(t)=1+t3.

In Figures  and , the numerical results display the values of E2w and Ew derived by 3PNDS, 5PNDS and 7PNDS with respect to different time levels tj, respectively, the measured data are added to four different error levels, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=0.5% in Example 4.2. It is shown that the accuracy of approximation w~δ(x,t) becomes higher, when the noise level δ is smaller. We show that the accuracy of approximation p~αiδ(t) is getting more accurate and more stable, when the noise level δ becomes smaller by and Table . Furthermore, the illustrated results of p~αiδ(t) are also satisfactory and effective with the reduction of spatial-step h.

Figure 5. The RMS error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=0.5% in Example 4.2. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 5. The RMS error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=0.5% in Example 4.2. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 6. The maximum error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=0.5% in Example 4.2. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 6. The maximum error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=0.5% in Example 4.2. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 7. The ratio plots of w1, w2 and w3 for 3PNDS, 5PNDS and 7PNDS, respectively, in Example 4.2, and the numbers ‘1,2,’ indicate the expressed sequence for different h1 above. (a) 3PNDS, (b) 5PNDS and (c) 7PNDS.

Figure 7. The ratio plots of w1, w2 and w3 for 3PNDS, 5PNDS and 7PNDS, respectively, in Example 4.2, and the numbers ‘1,2,…’ indicate the expressed sequence for different h1 above. (a) 3PNDS, (b) 5PNDS and (c) 7PNDS.

Figure 8. The exact function p(t) and its approximation for three difference schemes as 3PNDS, 5PNDS and 7PNDS with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=0.5% in Example 4.2. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 8. The exact function p(t) and its approximation for three difference schemes as 3PNDS, 5PNDS and 7PNDS with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=0.5% in Example 4.2. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Example 4.3

We consider the third problem (Equation1)–(Equation5): ϕ(x,t)=(π2+2t)exp(t)cos(πx)+2xtexp(t),f(x)=x+cos(πx),g0(t)=exp(t),g1(t)=0,x=2/3,E(t)=exp(t)/6, and Example 4.3 has the exact solution w(x,t)=exp(t)(cos(πx)+x),p(t)=12t.

The numerical results display the values of E2w and Ew derived by 3PNDS, 5PNDS and 7PNDS with respect to different time levels tj by Figures  and , respectively; the measured data are added to four different error levels, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.3. It is observed that the results of approximation solution w~δ(x,t) become more accurate, when the spatial-step h is smaller and the points number k ( k-step nodes) increases. Besides, the accuracy of w~δ(x,t) depends on the noise level δ according to the above three tests, if the selection of δ is a smaller value, the numerical solutions work effectively in this paper. In Figure  and Table , the numerical results of approximation p~αiδ(t) are stable by the choice of some appropriate factors, such as the noise level δ, the regularization parameters αiδ, the spatial-step h and the points number k.

Figure 9. The RMS error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.3. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 9. The RMS error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.3. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 10. The maximum error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.3. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 10. The maximum error of function w(x,t) obtained by 3PNDS, 5PNDS and 7PNDS, respectively, with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.3. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

In the theoretical analysis, when the regularization parameters αiδ=δ1/2, we can obtain the corresponding orders of convergence. As mentioned before, the rule is not a suitable choice for practical applications. We also give the numerical results in Tables and , which are obtained by the choice rule αiδ=δ1/2. In comparison, it is clear that the results obtained by the heuristic choice (Equation144) are better, as shown in Tables  and .

Table 5. The results of E2p and Ep obtained by 3PNDS, 5PNDS and 7PNDS, respectively, taking the case of αiδ=δ1/2 for Example 4.1.

Table 6. The results of E2p and Ep obtained by 3PNDS, 5PNDS, 7PNDS and RBFs, respectively, taking the case of αiδ in (Equation144) for Example 4.2.

Table 7. The results of E2p and Ep obtained by 3PNDS, 5PNDS and 7PNDS, respectively, taking the case of αiδ=δ1/2 for Example 4.2.

Table 8. The results of E2p and Ep obtained by 3PNDS, 5PNDS, 7PNDS and RBFs, respectively, taking the case of αiδ in (Equation144) for Example 4.3.

Table 9. The results of E2p and Ep obtained by 3PNDS, 5PNDS and 7PNDS, respectively, taking the case of αiδ=δ1/2 for Example 4.3.

In addition, according to Figures  and , we know that the inclination of the images, which are obtained by computing the ratios w1, w2 and w3 at the point (h1,τ1), gets bigger through three numerical approaches as 3PNDS, 5PNDS and 7PNDS, it means that the convergence rate of the presented methods should be generally improved with increasing of the points number k ( k-step nodes).

Figure 11. The ratio plots of w1, w2 and w3 for 3PNDS, 5PNDS and 7PNDS, respectively, in Example 4.3, and the numbers ‘1,2,’ indicate the expressed sequence for different h1 above. (a) 3PNDS, (b) 5PNDS and (c) 7PNDS.

Figure 11. The ratio plots of w1, w2 and w3 for 3PNDS, 5PNDS and 7PNDS, respectively, in Example 4.3, and the numbers ‘1,2,…’ indicate the expressed sequence for different h1 above. (a) 3PNDS, (b) 5PNDS and (c) 7PNDS.

Figure 12. The exact function p(t) and its approximation for three difference schemes as 3PNDS, 5PNDS and 7PNDS with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.3. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

Figure 12. The exact function p(t) and its approximation for three difference schemes as 3PNDS, 5PNDS and 7PNDS with four different noise levels added to the measured data, namely δ=0.001%,δ=0.01%,δ=0.1% and δ=1% in Example 4.3. (a) 3PNDS (h=1/10), (b) 5PNDS (h=1/10), (c) 7PNDS (h=1/10), (d) 3PNDS (h=1/15), (e) 5PNDS (h=1/15) and (f) 7PNDS (h=1/15).

5. Conclusion

In this paper, we construct a family of multistep numerical difference schemes, which are applied in the one-dimensional parabolic inverse problem with a control parameter p(t). Specially, we construct three difference schemes 3PNDS, 5PNDS and 7PNDS, and develop their corresponding truncation error estimations and convergence conclusions. It is well known that the problem of numerical differentiation with noisy errors is mildly ill-posed, we use the smoothing spline based on the Tikhonov regularization technique to obtain stable approximation solutions. The numerical results are given to observe that the accuracy of w~δ(x,t) and p~αiδ(t) depends on the selection of some factors, such as the noise level δ, the regularization parameters αiδ, the time-step τ, the spatial-step h and the points number k (k-step nodes). If the step sizes τ, h, and the noise level δ are decreased, thus the approximation solutions w~δ(x,t) and p~αiδ(t) can get more accurate results. If the step sizes h, τ, and the noise level δ are chosen much more smaller, thus the approximation solution p~αiδ(t) can get more accurate results, however, it would require more computational cost in this case, we will improve this in our further work. In conclusion, the presented numerical schemes are stable and effective.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China (Nos. 11471066, 11290143, 11572081), the Fundamental Research of Civil Aircraft (No. MJ-F-2012-04), the Fundamental Research Funds for the Central Universities (DUT15LK44), the Scientific Research Funds of Inner Mongolia University for the Nationalities (NMDYB17154) and PhD research startup foundation of Inner Mongolia University for Nationalities (BS424).

References

  • Hon YC, Wei T. A fundamental solution method for inverse heat conduction problem. Eng Anal Bound Elem. 2004;28:489–495. doi: 10.1016/S0955-7997(03)00102-4
  • Hansen PC. Regularization tools: a Matlab package for analysis and solution of discrete ill-posed problems. Numer Algorithms. 1994;6:1–35. doi: 10.1007/BF02149761
  • Cheng J, Yamamoto M. One new strategy for a priori choice of regularizing parameters in Tikhonov's regularization. Inverse Probl. 2000;16:31–38. doi: 10.1088/0266-5611/16/4/101
  • Tikhonov AN, Arsenin VY. Solutions of ill-posed problems. Washington (DC), New York: V.H. Winston & Sons, John Wiley & Sons; 1977.
  • Dehghan M. Numerical solution of one-dimensional parabolic inverse problem. Appl Math Comput. 2003;136:333–344.
  • Dehghan M. Finding a control parameter in one-dimensional parabolic equations. Appl Math Comput. 2003;135:491–503.
  • Dehghan M. An inverse problem of finding a source parameter in a semilinear parabolic equation. Appl Math Model. 2001;25:743–754. doi: 10.1016/S0307-904X(01)00010-5
  • Cannon JR, Lin YP, Xu SZ. Numerical procedures for the determination of an unknown coefficient in semi-linear parabolic differential equations. Inverse Probl. 1994;10:227–243. doi: 10.1088/0266-5611/10/2/004
  • Dehghan M. Parameter determination in a partial differential equation from the overspecified data. Math Comput Model. 2005;41:196–213. doi: 10.1016/j.mcm.2004.07.010
  • Mitchell AR, Griffiths DF. The finite difference method in partial differential equations. New York: Wiley; 1980.
  • Tatari M, Dehghan M. A method for solving partial differential equations via radial basis functions: application to the heat equation. Eng Anal Bound Elem. 2010;34:206–212. doi: 10.1016/j.enganabound.2009.09.003
  • Chen W, Tanaka M. A meshless, integration-free, and boundary-only RBF technique. Comput Math Appl. 2002;43:379–391. doi: 10.1016/S0898-1221(01)00293-0
  • Dehghan M, Tatari M. Determination of a control parameter in a one-dimensional parabolic equation using the method of radial basis function. Math Comput Model. 2006;44:1160–1168. doi: 10.1016/j.mcm.2006.04.003
  • Ma LM, Wu ZM. Radial basis functions method for parabolic inverse problem. Int J Comput Math. 2011;88(2):384–395. doi: 10.1080/00207160903452236
  • Arghand M, Amirfakhrian M. A meshless method based on the fundamental solution and radial basis function for solving an inverse heat conduction. Adv Math Phys. 2015;2015:1–8. Article ID 256726. doi: 10.1155/2015/256726
  • Hon YC, Wei T. The method of fundamental solution for solving multidimensional inverse heat conduction problems. Comput Model Eng Sci. 2005;7(2):119–132.
  • Dong CF, Sun FY, Meng BQ. A method of fundamental solutions for inverse heat conduction problems in an anisotropic medium. Eng Anal Bound Elem. 2007;31:75–82. doi: 10.1016/j.enganabound.2006.04.007
  • Shen SY. A numerical study of inverse heat conduction problems. Comput Math Appl. 1999;38:173–188. doi: 10.1016/S0898-1221(99)00248-5
  • Johansson T, Lesnic D. Determination of a spacewise dependent heat source. J Comput Appl Math. 2007;209:66–80. doi: 10.1016/j.cam.2006.10.026
  • Dehghan M. Implicit solution of a two-dimensional parabolic inverse problem with temperature overspecification. J Comput Anal Appl. 2001;3(4):383–398.
  • Dehghan M. Finite difference schemes for two-dimensional parabolic inverse problem with temperature overspecification. Int J Comput Math. 2000;75:339–349. doi: 10.1080/00207160008804989
  • Chantasiriwan S. An algorithm for solving multidimensional inverse heat conduction problem. Int J Heat Mass Transf. 2001;44:3823–3832. doi: 10.1016/S0017-9310(01)00037-0
  • Dehghan M. Determination of a control function in three-dimensional parabolic equations. Math Comput Simul. 2003;61:89–100. doi: 10.1016/S0378-4754(01)00434-7
  • Dehghan M, Abbaszadeh M, Mohebbi A. The numerical solution of nonlinear high dimensional generalized Benjamin–Bona–Mahony–Burgers equation via the meshless method of radial basis function. Comput Math Appl. 2014;68:212–237. doi: 10.1016/j.camwa.2014.05.019
  • Dehghan M, Mohammadi V. The numerical solution of Cahn–Hilliard (CH) equation in one, two and three-dimensions via globally radial basis functions (GRBFs) and RBFs-differential quadrature (RBFs-DQ) methods. Eng Anal Bound Elem. 2015;51:74–100. doi: 10.1016/j.enganabound.2014.10.008
  • Lakestani M, Dehghan M. The use of Chebyshev cardinal functions for the solution of a partial differential equation with an unknown time-dependent coefficient subject to an extra measurement. J Comput Appl Math. 2010;235:669–678. doi: 10.1016/j.cam.2010.06.020
  • Mohebbi A, Dehghan M. High-order scheme for determination of a control parameter in an inverse problem from the over-specified data. Comput Phys Comm. 2010;181:1947–1954. doi: 10.1016/j.cpc.2010.09.009
  • Dehghan M. Fourth-order techniques for identifying a control parameter in the parabolic equations. Int J Eng Sci. 2002;40(4):433–447. doi: 10.1016/S0020-7225(01)00066-0
  • Dehghan M, Shakeri F. Method of lines solutions of the parabolic inverse problem with an overspecification at a point. Numer Algorithms. 2009;50(4):417–437. doi: 10.1007/s11075-008-9234-3
  • Shamsi M, Dehghan M. Determination of a control function in three-dimensional parabolic equations by Legendre pseudospectral method. Numer Methods Partial Differ Equ. 2012;28(1):74–93. doi: 10.1002/num.20608
  • Saadatmandia A, Dehghan M. Computation of two time-dependent coefficients in a parabolic partial differential equation subject to additional specifications. Int J Comput Math. 2010;87(5):997–1008. doi: 10.1080/00207160802253958
  • Cannon JR, Duchateau P. An inverse problem for an unknown source in heat equation. J Math Anal Appl. 1980;75:465–485. doi: 10.1016/0022-247X(80)90095-5
  • Cannon JR, Lin Y, Wng S. Determination of source parameter in parabolic equation. Meccanica. 1992;27:85–94. doi: 10.1007/BF00420586
  • Macbain JA, Bendar JB. Existence and uniqueness properties for the one-dimensional magnetotellurics inversion problem. J Math Phys. 1986;27:645–649. doi: 10.1063/1.527219
  • Macbain JA. Inversion theory for a parametrized diffusion problem. SIAM J Appl Math. 1987;47:1386–1391. doi: 10.1137/0147091
  • Cannon JR, Lin Y. An inverse problem of finding a parameter in a semi-linear heat equation. J Math Anal Appl. 1990;145:470–484. doi: 10.1016/0022-247X(90)90414-B
  • Cannon JR, Lin YP. Determination of a parameter p(t) in some quasi-linear parabolic differential equations. Inverse Probl. 1988;4:35–45. doi: 10.1088/0266-5611/4/1/006
  • Cannon JR, Lin YP. Determination of parameter p(t) in Ho¨lder classes for some semilinear parabolic equations. Inverse Probl. 1988;4:595–606. doi: 10.1088/0266-5611/4/3/005
  • Li CJ. A kind of multistep finite difference methods for arbitrary order linear boundary value problem. Appl Math Comput. 2008;196:858–865.
  • Hanke M, Scherzer O. Inverse problems light: numerical differentiation. Am Math Mon. 2001;108(6):512–521. doi: 10.1080/00029890.2001.11919778
  • Wei T, Li M. High order numerical derivatives for one-dimensional scattered noisy data. Appl Math Comput. 2006;175:1744–1759.
  • Zhang YF, Li CJ. A Gaussian RBFs method with regularization for the numerical solution of inverse heat conduction problems. Inverse Probl Sci Eng. 2016;24(9):1606–1646. doi: 10.1080/17415977.2015.1131825
  • Engl HW, Hanke M, Neubauer A. Regularization of inverse problems. Dordrecht/Boston/London: Kluwer Academic Publishers; 1996.
  • Wang YB, Jia XZ, Cheng. J. A numerical differentiation method and its application to reconstruction of discontinuity. Inverse Probl. 2002;18:1461–1476. doi: 10.1088/0266-5611/18/6/301
  • Hào DN, Chuong LH, Lesnic. D. Heuristic regularization methods for numerical differentiation. Comput Math Appl. 2012;63:816–826. doi: 10.1016/j.camwa.2011.11.047
  • Süli E, Mayers D. An introduction to numerical analysis. Cambridge (UK): Cambridge University Press; 2003.
  • Adams RA. Sobolev spaces, pure and applied mathematics. Vol. 65. New York/London: Academic Press; 1975.
  • de Boor C. Spline toolbox for use with matlab. Natick: The MathWorks, Inc.; 1992.

Reprints and Corporate Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

To request a reprint or corporate permissions for this article, please click on the relevant link below:

Academic Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

Obtain permissions instantly via Rightslink by clicking on the button below:

If you are unable to obtain permissions via Rightslink, please complete and submit this Permissions form. For more information, please visit our Permissions help page.