Numerical solving the identification problem for the lower coefficient of parabolic equation1footnote 11footnote 1This work was supported by RFBR (project 13-01-00719)

# Numerical solving the identification problem for the lower coefficient of parabolic equation111This work was supported by RFBR (project 13-01-00719)

P.N. Vabishchevich22footnotemark: 2 V.I. Vasil’ev33footnotemark: 3 Nuclear Safety Institute, Russian Academy of Sciences, 52, B. Tulskaya, 115191 Moscow, Russia North-Eastern Federal University, 58, Belinskogo, 677000 Yakutsk, Russia
###### Abstract

In the theory and practice of inverse problems for partial differential equations (PDEs) much attention is paid to the problem of the identification of coefficients from some additional information. This work deals with the problem of determining in a multidimensional parabolic equation the lower coefficient that depends on time only. To solve numerically a nonlinear inverse problem, linearized approximations in time are constructed using standard finite element procedures in space. The computational algorithm is based on a special decomposition, where the transition to a new time level is implemented via solving two standard elliptic problems. The numerical results presented here for a model 2D problem demonstrate capabilities of the proposed computational algorithms for approximate solving inverse problems.

###### keywords:
Inverse problem, Control parameter, Parabolic partial differential equation, Finite element approximation, Difference scheme
###### Pacs:
02.30.Zz, 02.30.Jr
65J22, 65M32

## 1 Introduction

Mathematical modeling of many applied problems of science and engineering results in the numerical solution of inverse problems alifanov2011inverse (); tarantola1987inverse (); aster2011parameter (). Inverse problems often belong to the class of ill-posed (conditionally correct) problems, and therefore various regularization algorithms are employed to solve them numerically kirsch1996introduction (); tikhonov1977solutions (); engl2000regularization ().

Particular attention should be given to inverse problems for PDEs lavrentev1986ill (); isakov1998inverse (). In this case, a theoretical study includes the fundamental questions of uniqueness of the solution and its stability both from the viewpoint of the theory of differential equations express2000methods (); belov2002inverse () and from the viewpoint of the theory of optimal control for distributed systems maksimov2002dynamical (). Many inverse problems are formulated as non-classical problems for PDEs. To solve these problems approximately, emphasis is on the development of stable computational algorithms that take into account peculiarities of inverse problems vogel2002computational (); samarskii2007numerical ().

Among inverse problems for PDEs we distinguish coefficient inverse problems, which are associated with the identification of coefficients and/or the right-hand side of an equation using some additional information. When considering time-dependent problems, the identification of the coefficient dependences on space and on time is usually separated into individual problems isakov1998inverse (); express2000methods (). In some cases, we have linear inverse problems (e.g., identification problems for the right-hand side of an equation); this situation essentially simplify their study.

Much attention is paid to the problem of determining the lower coefficient of a parabolic equation of second order, where, in particular, the coefficient depends on time only. An additional condition is most often formulated as a specification of the solution at an interior point or as the average value integrated over the whole domain. The existence and uniqueness of the solution of such an inverse problem and well-posedness of this problem in various functional classes are examined, for example, in the works prilepko1985determination (); prilepko1987solvability (); macbain1986existence (); cannon1990inverse ().

Numerical methods for solving the problem of the identification of the lower coefficient of parabolic equations are considered in many works wang1989finite (); cannon1994numerical (); dehghan2005identification (); ye2007stability (); wang2012inverse (). In view of the practical use, we highlight separately studies dealing with numerical solving inverse problems for multidimensional parabolic equations dehghan2001numerical (); daoud2005splitting (). To construct computational algorithms for the identification of the lower coefficient of a parabolic equation, there is widely used the idea of transformation of the equation by introducing new unknowns that results in a linear inverse problem.

In this paper, for a multidimensional parabolic equation, we consider the problem of determining the lower coefficient that depends on time only. Approximation in space is performed using standard finite elements quarteroni2008numerical (); brenner2008mathematical (). The main features of the nonlinear inverse problem are taken into account via a proper choice of the linearized approximation in time. Linear problems at a particular time level are solved on the basis of a special decomposition into two standard elliptic problems. The paper is organized as follows. In Section 2, for a parabolic equation of second order, we formulate the inverse problem of the identification of the lower coefficient. The computational algorithm based on the linearization scheme is described in Section 3. Section 4 presents possibilities of the schemes with the second-order approximation in time. Numerical results for a model 2D inverse problem are discussed in Section 5.

## 2 Problem formulation

For simplicity, we restrict ourselves to a 2D problem. Generalization to the 3D case is trivial. Let and be a bounded polygon. The direct problem is formulated as follows. We search , such that it is the solution of the parabolic equation of second order:

The boundary and initial conditions are also specified:

 k(x)∂u∂n+g(x)u=0,x∈∂Ω,0
 u(x,0)=u0(x),x∈Ω, (3)

where is the normal to . The formulation (1)–(3) presents the direct problem, where the right-hand side, coefficients of the equation as well as the boundary and initial conditions are specified.

Let us consider the inverse problem, where in equation (1), the coefficient is unknown. An additional condition is often formulated as

 ∫Ωu(x,t)ω(x)dx=φ(t),0

where is a weight function. In particular, choosing (), where is the Dirac -function, from (4), we get

 u(x∗,t)=φ(t),0

We assume that the above inverse problem of finding a pair of from equations (1)–(3) and additional conditions (4) or (5) is well-posed. The corresponding conditions for existence and uniqueness of the solution are available in the above-mentioned works. In this paper, we consider only the numerical solution of these inverse problems omitting theoretical issues of the convergence of an approximate solution to the exact one.

From the nonlinear inverse problem we can proceed to the linear one. Suppose

 v(x,t)=χ(t)u(x,t),χ(t)=exp(∫t0p(θ)dθ).

Then from (1)–(3), we get

 k(x)∂v∂n+g(x)v=0,x∈∂Ω,0
 v(x,0)=u0(x),x∈Ω.

The additional conditions (4) and (5) to identify uniquely take the form

 ∫Ωv(x,t)ω(x)dx=χ(t)φ(t),0
 u(x∗,t)=χ(t)φ(t),0

The above transition from the nonlinear inverse problem to the linear one is in common use for numerical solving problems of identification. In our work, we focus on the original formulation of the inverse problem (1)–(4) (or (1)–(3), (5)) without going to the linear problem.

## 3 Computational algorithm

The inverse problem of determining the pair of is nonlinear. The standard approach is based on the simplest approximations in time and involves the iterative solution of the corresponding nonlinear problem for the evaluation of the approximate solution at a new level. In our work, we apply such approximations in time that lead to linear problems for evaluating the solution at the new time level.

Let us define a uniform grid in time

 ¯¯¯ωτ=ωτ∪{T}={tn=nτ,n=0,1,...,N,τN=T}

and denote . Finite element approximations in space are employed. In the polygon , we perform a triangulation and introduce for this computational grid a finite-dimensional space of finite elements.

Using the fully implicit scheme for approximation in time, we obtain the following variational problem:

 ∫Ωu0vdx=∫Ωu0vdx. (7)

The additional relations (4) and (5) take the form

 ∫Ωun+1ω(x)dx=φn+1, (8)
 un+1(x∗)=φn+1,n=0,1,...,N−1. (9)

To evaluate at the new time level the approximate solution from (6)–(8) or (6), (7), (9), some iterative procedures are necessary. In solving time-dependent problems, the solution slightly varies when it pass from the previous time level to the next one. This basic feature of time-dependent problems is widely used in numerical solving nonlinear problems through the application of linearization procedures. We use a similar approach for the numerical solution of the inverse problem that is concerned with the identification of the lower coefficient of a parabolic equation.

Instead of (6), we will solve the following equation:

In this case, the lower coefficient of the parabolic equation is taken at the upper time level, whereas the approximate solution is treated at the previous time level. Let us consider the solution procedure of the problem (7), (9), (10) in detail.

For the approximate solution at the new time level , we introduce the following decomposition samarskii2007numerical (); borukhov2000numerical ():

 un+1(x)=yn+1(x)+pn+1wn+1(x). (11)

To find , we employ the equation

The function is determined from

Using the decomposition (11)–(12), equation (10) holds automatically for any .

To evaluate , we apply the condition (9) (or (8)). The substitution of (11) into (9) yields

 pn+1=1wn+1(x∗)(φn+1−yn+1(x∗)). (14)

The fundamental point of applicability of this algorithm is associated with the condition . The auxiliary function is determined from the grid elliptic equation (13). The property of having fixed sign for is followed, in particular, from the same property of the solution at the previous time level . Such constraints on the solution can be provided by the corresponding restrictions on the input data of the inverse problem. In any case, this problem requires special and careful consideration. In this paper, we assume that the constraint is satisfied.

In solving problem (8)–(10), instead of (14), we have

 pn+1=1∫Ωwn+1ω(x)dx(φn+1−∫Ωyn+1ω(x)dx) (15)

under the condition that

 ∫Ωwn+1ω(x)dx≠0.

In this case, additional restrictions are formulated on the function , e.g., its fixed sign in .

Thus, the computational algorithm for solving the inverse problem (1)–(4) (or (1)–(3), (5)) based on the linearized scheme (7), (8), (10) (or (7),(9), (10)) involves the solution of two standard grid elliptic equations for the auxiliary functions (equation (12)) and (equation (13)), the further evaluation of from (15) (or (14)), and the final calculation from the relation (11).

## 4 Scheme of the second-order accuracy

The nonlinear inverse problem (1)–(4) is characterized by a quadratic nonlinearity. When using the scheme with linearization (7), (8), (10), the nonlinear term is approximated with the first order with respect to . It is possible to apply the linearized scheme of second order. Let us consider the approximation

 a(tn+1/2)b(tn+1/2)=12a(tn+1)b(tn)+12a(tn)b(tn+1)+O(τ2).

Approximation of equation (1) with the boundary conditions (2) using the Crank-Nicolson scheme yields the linearized scheme

The scheme (7), (8), (16) belongs to the class of linearized schemes. In comparison with the scheme (7), (8), (10), it has a higher order of accuracy in time.

To implement (7), (8), (16), we again use the decomposition (11). In this case, for , we have

The auxiliary function is defined as the solution of the equation

Further, as in the case of the first-order scheme, we employ (15).

The Crank-Nicholson scheme for numerical solving the direct problems for parabolic equations is not very often used in computational practice. It is inferior to the fully implicit scheme in sense of conservation of monotonicity (fulfilment of the maximum principle for the grid problem), it has poor asymptotic properties for solving problems with large integration time, and it is not unconditionally SM-stable scheme samarskiaei1995computational (); vabishchevich2012sm (). For this reason, it is appropriate to consider another variant of linearization of this inverse problem, where the second-order approximation is applied only for the nonlinear term. In this case, instead of (10), (16), we put

The numerical implementation of the scheme (7), (8), (19) is performed in the standard way using the decomposition (11).

## 5 Numerical examples

To demonstrate possibilities of the above linearization schemes for solving the problem of the identification of the lower coefficient of the parabolic equation, we consider a 2D model problem. In the examples below, we put

 k(x)=1,f(x,t)=0,u0(x)=1,x∈Ω,g(x)=10,x∈∂Ω.

The problem is considered on a triangular grid, which consists of 1,180 nodes (2,230 triangles) and is shown in Fig.1. Here is the trapezoid with the vertices coordinates . The calculations were carried out for . The coefficient is taken in the form

 p(t)={1000t,0

The solution of the direct problem (1)–(3) at the observation point is depicted in Fig.2. It was otained using the fully implicit scheme with different time steps. The solution at the final time moment is presented in Fig.3.

The results of solving the inverse problem with variuos grids in time are shown in Fig.4. The solution of the direct problem obtained with is used as the input data (the function in the condition (5)). It is easy to see that the approximate solution of the inverse problem converges with decreasing the time step. These results were obtained using the first-order scheme (7), (9), (10).

Numerical results obtained for the above problem using the second-order scheme (7), (9), (16) are shown in Fig.5. For the discontinuous right-hand side (20), we observe characteristic wiggles of the identified coefficient. Such oscillations of the approximate solution are typical for the scheme (7), (9), (19).

If the desired solution (the coefficient ) is smooth, then the effect of using the second-order approximation is clearly expressed. As an example, we present the results of numerical solving the inverse problem, where the lower coefficient (the exact solution) has the form

 p(t)=1000t1+500t2.

The approximate solution obtained via the first-order scheme (7), (9), (10) is shown in Fig.6, whereas Fig.7 demonstrates the computations conducted by means of the second-order scheme (7), (9), (16).

## References

• (1) O. M. Alifanov, Inverse Heat Transfer Problems, Springer, 2011.
• (2) A. Tarantola, Inverse problem theory: methods for data fitting and model parameter estimation, Elsevier, 1987.
• (3) R. C. Aster, B. Borchers, C. H. Thurber, Parameter Estimation and Inverse Problems, Elsevier Science, 2011.
• (4) A. Kirsch, An introduction to the mathematical theory of inverse problems, Springer, 1996.
• (5) A. N. Tikhonov, V. Y. Arsenin, Solutions of ill-posed problems, Winston, 1977.
• (6) H. W. Engl, M. Hanke, A. Neubauer, Regularization of inverse problems, Kluwer Academic Publishers, 2000.
• (7) M. M. Lavrent’ev, V. G. Romanov, S. P. Shishatskii, Ill-posed Problems of Mathematical Physics and Analysis, American Mathematical Society, 1986.
• (8) V. Isakov, Inverse problems for partial differential equations, Springer, 1998.
• (9) A. I. Prilepko, D. G. Orlovsky, I. A. Vasin, Methods for Solving Inverse Problems in Mathematical Physics, Marcel Dekker, Inc, 2000.
• (10) Y. Y. Belov, Inverse problems for partial differential equations, VSP, 2002.
• (11) V. I. Maksimov, Dynamical inverse problems of distributed systems, VSP, 2002.
• (12) C. R. Vogel, Computational Methods for Inverse Problems, no. 10, Society for Industrial and Applied Mathematics, 2002.
• (13) A. A. Samarskii, P. N. Vabishchevich, Numerical Methods for Solving Inverse Problems of Mathematical Physics, De Gruyter, 2007.
• (14) A. I. Prilepko, D. G. Orlovskii, Determination of evolution parameter of an equation, and inverse problems in mathematical physics, part I and II, Differ. Equat. 21 (1, 4) (1985) 119–129, 694–701, in Russian.
• (15) A. I. Prilepko, V. V. Soloev, Solvability of the inverse boundary value problem of finding a coefficient of a lower order term in a parabolic equation, Differ. Equat. 23 (1) (1987) 136–143, in Russian.
• (16) J. A. MacBain, J. B. Bednar, Existence and uniqueness properties for the one-dimensional magnetotellurics inversion problem, Journal of mathematical physics 27 (1986) 645–649.
• (17) J. R. Cannon, Y. Lin, An inverse problem of finding a parameter in a semi-linear heat equation, Journal of Mathematical Analysis and Applications 145 (2) (1990) 470–484.
• (18) S. Wang, Y. Lin, A finite-difference solution to an inverse problem for determining a control function in a parabolic partial differential equation, Inverse problems 5 (4) (1989) 631–640.
• (19) J. R. Cannon, Y. Lin, S. Xu, Numerical procedures for the determination of an unknown coefficient in semi-linear parabolic differential equations, Inverse Problems 10 (2) (1994) 227–243.
• (20) M. Dehghan, Identification of a time-dependent coefficient in a partial differential equation subject to an extra measurement, Numerical Methods for Partial Differential Equations 21 (3) (2005) 611–622.
• (21) C. Ye, Z. Sun, On the stability and convergence of a difference scheme for an one-dimensional parabolic inverse problem, Applied mathematics and computation 188 (1) (2007) 214–225.
• (22) W. Wang, B. Han, M. Yamamoto, Inverse heat problem of determining time-dependent source parameter in reproducing kernel space, Nonlinear Analysis: Real World Applications 14 (1) (2013) 875–887.
• (23) M. Dehghan, Numerical methods for two-dimensional parabolic inverse problem with energy overspecification, International journal of computer mathematics 77 (3) (2001) 441–455.
• (24) D. S. Daoud, D. Subasi, A splitting up algorithm for the determination of the control parameter in multi dimensional parabolic problem, Applied mathematics and computation 166 (3) (2005) 584–595.
• (25) A. Quarteroni, A. Valli, Numerical approximation of partial differential equations, Springer, 2008.
• (26) S. C. Brenner, L. R. Scott, The mathematical theory of finite element methods, Springer, 2008.
• (27) V. T. Borukhov, P. N. Vabishchevich, Numerical solution of the inverse problem of reconstructing a distributed right-hand side of a parabolic equation, Computer physics communications 126 (1) (2000) 32–36.
• (28) A. A. Samarskii, P. N. Vabishchevich, Computational heat transfer, Wiley, 1995.
• (29) P. N. Vabishchevich, SM-stability of operator-difference schemes, Computational Mathematics and Mathematical Physics 52 (6) (2012) 887–894.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters   