A CCD-ADI Method for the Two-dimensional Parabolic Inverse Problem with an Overspecification at a point

Ling Weiwei Pan Kejia

(1. School of Mathematics and Statistics, Central South University, Changsha 410083, China;2. College of Social Management, Jiangxi College of Applied Technology,Ganzhou, Jiangxi 341000, China)

Abstract Many engineering problems can be modelled by parabolic differential equations with unknown parameters. It is often quite important to develop highly accurate numerical methods for solving these inverse problems. In this paper, we propose a linearized three-level combined compact difference (CCD) scheme with the alternating direction implicit (ADI) method for solving the two-dimensional unsteady reaction-diffusion equation with a control parameter. The proposed method is sixth-order accurate in space and second-order accurate in time. The resulting linear system at each ADI step is a block tridiagonal system which can be solved effectively by the block Thomas algorithm. In addition, we rigorously prove the existence and the uniqueness of the solution of the CCD-ADI method under the periodic boundary conditions. Finally, numerical tests are carried out to show the unconditional stability, the accuracy and efficiency of the proposed method compared with the existing spatial fourth-order method.

Key words Reaction-diffusion equation, Alternating direction implicit, Combined compact difference scheme, Solvability, Control parameter

1 Introduction

In this paper, the following inverse problem will be considered, i.e., find the control parameterp(t) and the solutionu(x,y,t) of the unsteady two-dimensional reaction-diffusion equation

ut=uxx+uyy+p(t)u+f(x,y,t),0≤x≤1,0≤y≤1,0

(1)

with initial condition

u(x,y,0)=φ(x,y),0≤x≤1,0≤y≤1

(2)

and boundary conditions

u(0,y,t)=gl(y,t),u(1,y,t)=gr(y,t),0≤t≤T,0≤y≤1,

u(x,0,t)=hl(x,t),u(x,1,t)=hr(x,t),0≤t≤T,0≤x≤1

(3)

subject to the over-specified data at a given point (x*,y*) in the spatial domain [1-5,13,26]

u(x*,y*,t)=E(t),0≤t≤T,

(4)

wheref(x,y,t),φ(x,y),gl(y,t),gr(y,t),hl(x,t),hr(x,t) andE(t) are sufficiently smooth given functions, |E(t)|≥E0>0, whileu(x,y,t) andp(t) are unknown functions which is assumed to be sufficiently smooth.

In applied science and engineering, the heat transfer process with a source parameter by given a measured dataE(t) at a fixed location (x*,y*) can be described by the above control problem (1)-(4). Many branches of physics, such as fluid dynamics, plasma dynamics, optics, field theory, condensed matter physics, are interested in the above kind of equations [14].

However, to our best knowledge, the spatial convergence rates of all previous methods in the literature are not exceed than 4. In this paper, we will propose a CCD-ADI method for solving the above two-dimensional problem (1)-(4). The CCD scheme is a three-point sixth-order implicit scheme, which is first proposed by [10], and the linear system from the CCD scheme is a block tridiagonal system which can be effectively solved by the block Thomas algorithm, see Remark 1 of this paper for details. Sengupta et al [38] proposed a new combined stable and dispersion relation preserving compact scheme and discussed the corresponding dissipation and de-aliasing properties in [39]. The CCD method has been successfully applied to solve many kinds of equations, which include the unsteady convection-diffusion equations [40], the linear second-order partial differential equations with mixed derivative [27], the nonlinear Schrödinger equations [28], the time-fractional advection-diffusion equations [11,20,21], the linear and nonlinear telegraph equations [23,8], the shallow water equation [35,44], the Camassa-Holmequation [45], the coupled nonlinear Burgers’ equation[7], and the Navier-Stokes equations [24,46,9].

The objective of this paper is to present a spatial high-order (sixth-order) and temporal second-order method for solving the parabolic inverse problem (1)-(4). Section 2 presents a CCD-ADI method for the above two-dimensional parabolic control problem, section 3 proves the existence and the uniqueness of the difference solution under periodic boundary conditions, section 4 gives some numerical results and comparisons with the existing fourth-order method, and the conclusions are provided in the final section.

2 Numerical methods

2.1 Combined compact difference scheme

We start with the second-order ordinary differential equation

γ(x)uxx+q(x)ux+r(x)u=s(x),a≤x≤b,

(5)

with boundary conditions

u(a)=ul,u(b)=ur,

(6)

whereγ(x),q(x),r(x),s(x) are given functions, anda,b,ul,urare constants.

(7)

(8)

At the two boundary pointsx0andxM, a pair of fifth-order one-sided CCD schemes are introduced in order to keep three-point structure[10],

(9)

(10)

In addition, the equation (5) and boundary conditions (6) are also used for the CCD scheme as follows：

γ(xi)U″i+q(xi)U′i+r(xi)Ui=s(xi), fori=0,1,…,M，

(11)

U0=ul,UM=ur.

(12)

Eqs. (7)-(12) gives a three-point CCD scheme. The previous numerical results show that, although the fifth-order boundary conditions are used at boundaries, sixth-order accuracy of the CCD method is guaranteed[10]. Furthermore, since the coefficient matrix of the CCD system possesses a block tridiagonal structure, it can be efficiently solved by the block Thomas algorithm, see the following Remark 1.

2.2 A linearized three-level CCD-ADI scheme

ut=Lxu+Lyu+p(t)u+f(x,y,t).

(13)

Using the three-level Crank-Nicolson scheme to discretize the time derivative att=tnin (13), we have

(14)

Collecting the terms forun+1andun-1in (14) yields

(1-Δtpn-ΔtLx-ΔtLy)un+1=(1+Δtpn+ΔtLx+ΔtLy)un-1+2Δtf(x,y,nΔt)+O(Δt3).

(15)

Eq. (15) can be rewritten in the factorial form

(16)

Since

(17)

Eq. (16) is indeed

(18)

Dropping the truncation error term in (18), evaluating at (xj,yk), and introducing an intermediate variableu*, Eq. (18) can be rewritten as

(19)

(20)

where

(21)

Since (19) and (20) are one-dimensional problems, the CCD method can be applied directly. Therefore, Eqs. (19) and (20) can be solved with the accuracy ofO(Δt2+Δx6+Δy6).

When computing the right-hand side of (19), all the second derivatives of the known valueun-1are required to be computed with sixth-order accuracy beforehand [40,23]. The method for calculating the first and second derivatives from a given function by the CCD method [10,40] is illustrated shortly as follows. Besides the internal formulas (7)-(8) and the boundary formulas (9)-(10), two additional formulas at boundary points[10]:

(22)

(23)

are required to compute the all 2(M+1) unknownsU′i,U″i,i=0,…,Mfrom the givenUi,i=0,…,M.

(24)

(25)

Andp0can be obtained by

(26)

Using the Taylor expansion and (1)-(2), we have

u(x,y,Δt)=u(x,y,0)+Δtut(x,y,0)+O(Δt2)

=φ(x,y)+Δt(φxx+φyy+p(0)φ+f(x,y,0))+O(Δt2),

(27)

and

uxxt(x,y,0)=uxxxx(x,y,0)+uyyxx(x,y,0)+p(0)uxx(x,y,0)+fxx(x,y,0)

=φxxxx+φxxyy+p(0)φxx+fxx(x,y,0),

(28)

uyyt(x,y,0)=uxxyy(x,y,0)+uyyyy(x,y,0)+p(0)uyy(x,y,0)+fyy(x,y,0)

=φyyxx+φyyyy+p(0)φyy+fyy(x,y,0).

(29)

Therefore,

uxx(x,y,Δt)=uxx(x,y,0)+Δtuxxt(x,y,0)+O(Δt2)

=φxx(x,y)+Δt(φxxxx+φxxyy+p(0)φxx+fxx(x,y,0))+O(Δt2),

(30)

uyy(x,y,Δt)=uyy(x,y,0)+Δtuyyt(x,y,0)+O(Δt2)

=φxx(x,y)+Δt(φyyxx+φyyyy+p(0)φyy+fyy(x,y,0))+O(Δt2).

(31)

By using the above equations, we can obtain the second-order temporal accurate approximations as follows:

(32)

(33)

(34)

forj=0,…,Mx,k=0,…,My.

From (1) we know that

(35)

From the above spatial and time discretization, the numerical approach is second-order accurate in time variable and sixth-order accurate in space variable, which will be confirmed by the numerical results in the later section.

The following algorithm 1 gives the details of the CCD-ADI method for solvingun+1andpn+1.

Algorithm 1 CCD-ADI algorithm for the parabolic inverse problem1:Compute u1j,k,(uxx)1j,k,(uyy)1j,k for all mesh points by (32)-(34) and p1 by (35). 2:n=1 to N-13:Step 1: Compute the right-hand side of (19) for all grid points.4:a) Compute gn-1j,k:=(1+Δtpn)un-1j,k+Δt(un-1yy)j,k for all mesh points.5:b) Compute (gx)n-1j,k and (gxx)n-1j,k (j=0,…,Mx) by (7)-(10) and (22)-(23) for k=1,…,My. 6:c) Compute hnj,k:=(1+Δtpn)gn-1j,k+Δt(gxx)n-1j,k for all mesh points.7:d) Compute dnj,k:=1-Δtpn1+Δtpnhnj,k+2Δt(1-Δtpn)fnj,kfor all mesh points. 8:Step 2: Compute the boundary conditionsu*j,k in (19) for j=0,Mx.9:a) Compute (uy)n+1j,k and (uyy)n+1j,k (k=0,…,My) by (7)-(10) and (22)-(23). 10:b) Compute u*j,k:=(1-Δtpn)un+1j,k-Δt(uyy)n+1j,k (k=0,…,My), for j=0,Mx. 11:Step 3: Solve u*j,k(j=0,…,Mx) from (1-Δtpn-ΔtLx)u*j,k=dnj,k and the CCD scheme (7)-(12) for k=0,…,My. 12:Step 4: Solve un+1j,k and (un+1yy)j,k(k=0,…,My) from (1-Δtpn-ΔtLy)un+1j,k=u*j,k and the CCD scheme (7)-(12) for j=0,…,Mx. 13:Step 5: Compute (ux)n+1j,k and (uxx)n+1j,k (j=0,…,Mx) by (7)-(10) and (22)-(23) for k=0,…,My. 14:Step 6: Solve pn+1 through (24).15:end for

3 Solvability

In this section, we provide a theoretical proof for the existence and the uniqueness of the numerical solutionun+1in (19)-(20) under the periodic boundary conditions.

Lemma 1

(36)

where

and

ProofThis lemma can be verified through direct computation.

Lemma 2[22] For any two given circulant matricesAandB, the sumA+Bis circulant, the productABis circulant, andAB=BA.

Lemma 3SupposeA,B,C,Dare alln×nmatrices, ifAC=CAandAis non-singular, then the following determinant can be simplified as

ProofThis lemma can be found in any elemental Linear algebra text book.

Lemma 4Under the periodic boundary conditions, the right-hand-side of (19) is uniquely determined through the CCD scheme (7)-(8) with the given valuesun,pnandun-1.

(37)

where

the unknown vectors

SinceA1,1is strictly diagonally dominant,A1,1is a non-singular matrix. Using Lemma 1-3 we have

Theorem 1Under the periodic boundary conditions, the CCD-ADI scheme (19)-(20) has a unique solutionun+1. In addition,pn+1is uniquely determined through (24).

ProofFrom Lemma 4, the right-hand-side of (19) is uniquely determined. We first prove that (19) has a unique solutionu*under the periodic boundary conditions.

(38)

where

the unknown vectors

IMxis theMx×Mxidentity matrix, anddnis the discrete vector from the right-hand-side of (19).

From Laplace’s expansion theorem and Lemma 3, we have

where we have used that

and

is strictly diagonally dominant and nonsingular. Then we know that

Thus, the linear system (38) has a unique solution, andu*is uniquely determined from (19).

Remark 1While implementing the CCD scheme to solve (5), we recommend not to write the CCD stencils given by Eqs. (7)-(12) as 3×3 block linear system with a relatively large bandwidth that involves large computations. The above Eq. (38) is used for analysis only. Instead, it is advisable to rewrite the CCD stencils given by Eqs. (7)-(12) as[39]

Au=b,

(39)

whereAis a block tridiagonal matrix with the form

(40)

and the unknown vectors

u=[U′1,U″1,U1,U′2,U″2,U2,…,U′M+1,U″M+1,UM+1]T.

(41)

The 3×3 sub-matricesAk,BkandCkare defined as

(42)

for 2≤k≤M. The block tridiagonal linear system (40) can be solved efficiently by the block Thomas algorithm, see [25].

Remark 2We note that ifpis a known bounded function (forward problem), or we make an assumption that the finite difference approximation ofpfor the inverse problem is uniformly bounded in the entire computational time, it can be shown that the CCD-ADI scheme is unconditionally stable following a similar Fourier stability analysis in [40] or [23]. When the termp(t)uis included, the unconditionally stability of the CCD-ADI scheme for the inverse problem, which is not easy to prove theoretically, is verified numerically in section 5.

4 Numerical results

In this section, numerical results for the CCD-ADI scheme presented in the section 2 are presented. For comparison propose, the method of [26] is also carried out. We wrote the program in Matlab and tested all numerical examples on a desktop equipped with Intel(R) Core(TM) i7-4790K CPU (4.00 GHz) and 8GB RAM using double-precision arithmetics.

Consider the two-dimensional control problem in [26]:

(43)

with initial condition

(44)

and boundary conditions

(45)

subject to the over-specification at a point in the spatial domain

(46)

The exact solution is given by

(47)

In the numerical computation, we setMx=Mysuch thatΔx=Δyh. Denote

and

Table 1 Maximum errors and the CPU seconds of the CCD-ADI scheme with time step Δt=25h3

SupposeΔt=βh3, whereβis a positive constant. If

E∞(h,βh3)≈c1hp1,F∞(h,βh3)≈c2hp2,

(48)

then

logE∞(h,βh3)≈logc1+p1(logh), logF∞(h,βh3)≈logc2+p2(logh). (49)

Table 3 Maximum errors and the CPU seconds of the CCD-ADI scheme with time step Δt=5h3

Table 1-3 giveE∞(h,Δt) andF∞(h,Δt), the corresponding convergent rates, and total computational time (CPU seconds required) using different spatial mesh sizes and time steps, where the convergent rates are defined as

(50)

Table 1, Table 2 and Table 3 correspond toβ=25,β=10 andβ=5, respectively. As we can see that the maximum errors of both numerical solutionsuandpconverge with rate 6 when the spatial mesh sizehis reduced by a factor of 2 and the temporal mesh sizeΔtis reduced by a factor of 8. This implies that the proposed CCD-ADI method in this paper is second-order accurate in time variable and sixth-order accurate in space variable.

Table 4 Maximum errors and the CPU seconds of the HOC-ADI scheme with time step Δt=h2

Figure 1 Numerical results and absolute errors of u at time T=1 with h=1/80 and Δt=1/51200

In order to investigate the convergence order of the method numericallly, using the data of Table 1, we obtain the following linear fitting functions

logE∞(h,βh3)≈11.0609+5.9766logh,

(51)

logF∞(h,βh3)≈12.8028+5.9777logh.

(52)

Figure 1 presents the numerical results and absolute errors ofuat timeT=1 withh=1/80 andΔt=1/51200. We also plot the absolute error ofpversustwithh=1/80 andΔt=1/51200 in Figure 2. As we can see, the proposed method produces satisfactory results.

In order the illustrate the unconditionally stability of the proposed method, we fixh=1/20 and varyΔtfrom 1/100000 to 1/10, numerical results ofF∞(h,Δt) versusN=1/Δtare plotted in Figure 3. As one can seethat these results clearly show that the time step is not related to the spatial meshsize, and when the temporal meshsize becomes small enough, the dominant error comes from the spatial part.

Figure 2 Numerical results and absolute errors of p with h=1/80 and Δt=1/51200

Figure 3 F∞(h, Δt) obtained by our CCD-ADI method with fixed h=1/20 but different Δt varing from 1/100000 to 1/10.

For comparison proposes, we now also give the results and show the performance using the method of [26]. Table 4 gives the maximum errorsE∞(h,Δt) andF∞(h,Δt) as well as the convergent rate with different spatial mesh sizes and time steps, where the convergent rate is defined as

(53)

As we can see that the maximum errors of both numerical solutionsuandpconverge with order 4 when the spatial mesh sizehis reduced by a factor of 2 and the temporal mesh sizeΔtis reduced by a factor of 4. This confirms that the method [26] is second-order convergent in time variable and fourth-order convergent in space variable.

Finally, from Table 2 and Table 4, one can find that both the CCD-ADI method and the fourth-order HOC-ADI method in [26] have the same starting meshsize (h=1/10,Δt=1/100), however, when refining to the same time stepsize (Δt=1/409600), the computing time used by the CCD-ADI method is only about one third of that by the HOC-ADI method. Moreover, the fourth-order method actually loses its theoretical order of convergence when reachingΔt=1/409600, and the errors for bothuandpare much larger than the corresponding errors obtained by the CCD-ADI method.

5 Conclusions

A linearized CCD-ADI method is developed for the two-dimensional parabolic control problem. The proposed method is second-order accurate in time variable and sixth-order accurate in space variable. The unique solvability of the difference scheme is obtained under the periodic boundary conditions. Comparing to the existing numerical methods in the literature, current CCD-ADI method improves the spatial accuracy greatly. Numerical results are presented to illustrate the unconditional stability, the accuracyand efficiency of the method.