Solving of the coefficient inverse problems for a nonlinear singularly perturbed reaction-diffusion-advection equation with the final time data

Accepted Manuscript

Solving of the coefficient inverse problems for a nonlinear singularly perturbed reaction-diffusion-advection equation with the final time data D.V. Lukyanenko, M.A. Shishlenin, V.T. Volkov PII: DOI: Reference:

S1007-5704(17)30215-0 10.1016/j.cnsns.2017.06.002 CNSNS 4222

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Received date: Revised date: Accepted date:

13 August 2016 24 December 2016 3 June 2017

Please cite this article as: D.V. Lukyanenko, M.A. Shishlenin, V.T. Volkov, Solving of the coefficient inverse problems for a nonlinear singularly perturbed reaction-diffusion-advection equation with the final time data, Communications in Nonlinear Science and Numerical Simulation (2017), doi: 10.1016/j.cnsns.2017.06.002

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • Asymptotic method for solving coefficient inverse problem is proposed.

CR IP T

• Strategy to generate a dynamic adapted mesh is described.

AC

CE

PT

ED

M

AN US

• A priori information for the construction is provided by the asymptotic analysis.

1

ACCEPTED MANUSCRIPT

CR IP T

Solving of the coefficient inverse problems for a nonlinear singularly perturbed reaction-diffusion-advection equation with the final time data✩ D.V. Lukyanenkoa,∗, M.A. Shishleninb,c,d , V.T. Volkova a

AN US

Department of Mathematics, Faculty of Physics, Lomonosov Moscow State University, Moscow 119991, Russia b Sobolev Institute of Mathematics, Novosibirsk 630090, Russia c Institute of Computational Mathematics and Mathematical Geophysics, Novosibirsk 630090, Russia d Novosibirsk State University, Novosibirsk 630090, Russia

Abstract

CE

PT

ED

M

We propose the numerical method for solving coefficient inverse problem for a nonlinear singularly perturbed reaction-diffusion-advection equation with the final time observation data based on the asymptotic analysis and the gradient method. Asymptotic analysis allows us to extract a priory information about interior layer (moving front), which appears in the direct problem, and boundary layers, which appear in the conjugate problem. We describe and implement the method of constructing a dynamically adapted mesh based on this a priory information. The dynamically adapted mesh significantly reduces the complexity of the numerical calculations and improve the numerical stability in comparison with the usual approaches. Numerical example shows the effectiveness of the proposed method.

AC

Keywords: singularly perturbed problem, interior and boundary layers, ✩

The work was supported by RFBR (projects NNo. 17-01-00519, 17-01-00670, 17-0100159, 16-01-00755, 16-29-15120 and 16-01-00437), International Mathematical Center of Novosibirsk State University and the Ministry of Education and Science of the Russian Federation. ∗ Corresponding author Email addresses: [email protected] (D.V. Lukyanenko), [email protected] (M.A. Shishlenin), [email protected] (V.T. Volkov) Preprint submitted to Communications in Nonlinear Science and Numerical SimulationJune 5, 2017

ACCEPTED MANUSCRIPT

CR IP T

dynamically adapted mesh, reaction-diffusion-advection equation, coefficient inverse problem, final time observed data. 2000 MSC: 65M32, 65L04, 65L12, 65L20, 65M20, 35G31 1. Introduction

AC

CE

PT

ED

M

AN US

In this paper the numerical method for solving the coefficient inverse problem for singularly perturbed reaction-diffusion-advection (RDA) equations with the final time observation data is developed. The method is based on the asymptotic analysis and optimization method. Due to the nonlinearity of the considered direct and inverse problems we minimize the cost functional by the gradient method to find approximate solution of the coefficient inverse problem. On each step of the gradient method the direct and corresponding conjugate problems are solved. The statements of these problems have the small parameter at higher derivative. Therefore the solutions of direct and conjugate RDA problems often feature narrow boundary and interior layers (stationary or moving fronts) and are extremely difficult for a numerical treatment. The main idea of our approach is that the asymptotic analysis of the direct problems allows to reduce significantly the complexity and optimize the numerical calculations, save the computational resources and improve the numerical stability of solving of corresponding singularly perturbed inverse problems. The asymptotic analysis allows us to locate boundary and interior layers and to develop an efficient algorithm due to the following reasons: the smaller the parameter, the more rough and unstable numerical solution we obtain and more precise a priori information about the exact solution we get from asymptotic analysis. This fact gives the possibility for a productive combination of asymptotic and numerical approaches. Asymptotic analysis also allows to prove the existence of the solution of certain class for the direct problem and highlight a priori information for solving of the inverse problem. Moreover, if consider the solutions with the internal layers or moving fronts, the spatial dimension of some problems for numerical calculations can be decreased. For example, to obtain a priori information about the location or the speed of the internal layer (moving front) we solve the equations with less spatial dimension than the original partial differential equation (PDE) problem. Namely, if the original PDE problem is n–dimensional, the equation for the moving front location is (n − 1)–dimensional (this equation can be even ordinary differential or algebraic). 3

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

CR IP T

We demonstrate our approach on solving of initial-boundary problem for the following singularly perturbed reaction-diffusion-advection equation (ε  1): ∂u ∂ 2 u ∂u = −u + q(x) u. (1.1) ε 2− ∂x ∂t ∂x This type of equation or systems are used in mathematical models of biology, chemical kinetics, theory of combustion e.t.c. The asymptotic analysis of the direct problem for equation (1.1) has been already performed in [1, 2]. Later it was applied to construct the dynamically adapted mesh based on the a priory information about the solution [3]. Different types of using a priori information in coefficient inverse problems were considered in [4, 5, 6] but in this paper we represent completely new approach that allows to solve inverse problems for singularly perturbed equations more effectively. The paper is structured as follows. In Section 2 we discuss the general statement of the coefficient inverse problem which solution usually can not have any special features, however during the solution of relevant direct and conjugate problems, the difficulties connected with the interior and boundary layers may arise. Then, in Section 3, we formally mention the general method of its solution. In Section 4 we describe the main ideas of asymptotic theory that allows to get a priori information for the dynamically adapted mesh construction. In Section 5 we explain some nuances of the dynamically adapted mesh constructing. In Section 6 we perform a numerical experiment for some particular example to demonstrate the effectiveness of the proposed asymptotic-numerical approach. 2. Problem formulation

AC

CE

Let us consider the direct problem  2 ∂u ∂ u ∂u  ε 2 − = −u + q(x) u, x ∈ (0, 1), t ∈ (0, T ),  ∂x ∂t ∂x u(0, t) = ul (t), u(1, t) = ur (t), t ∈ (0, T ),    u(x, 0) = u0 (x), x ∈ (0, 1).

(2.1)

Inverse problem consists of finding function q(x) by known additional information: u(x, T ) = fobs (x), x ∈ (0, 1). (2.2) 4

ACCEPTED MANUSCRIPT

AN US

CR IP T

We apply the optimal control problem for the numerical solution (2.1)– (2.2). Ill-posedness analysis Let us show that the inverse problem (2.1)–(2.2) is ill-posed [7]. We will analyse two formulations of inverse problem: linear and nonlinear. Let us investigate ill-posedness of the problem of recovering the initial data by known final time measured data (reverse time problem). We will construct the example of instability. The problem of recovering the initial data is more ”simple” in comparison with the original problem formulation. Therefore, if we investigate the ill-posedness of reverse time problem, we can expect that the instability effects will influence harder in the numerical algorithm of the solution of the nonlinear coefficient inverse problem. The linear analog of the problem (2.1) is:  2 ∂u ∂ u ∂u   x ∈ (0, 1), t ∈ (0, T ),  ε ∂x2 − ∂t = − ∂x , (2.3) u(0, t) = 0, u(1, t) = 0, t ∈ (0, T ),    u(x, 0) = q(x), x ∈ (0, 1).

ED

M

Inverse problem consists of finding function q(x) by known additional information: u(x, T ) = fobs (x), x ∈ (0, 1). (2.4) The solution of the direct problem (2.3) is given by formula:

PT

u(x, t) =

∞ X

x

qn e 2ε −

4ε2 π 2 n2 +1 t 4ε

sin(πnx).

(2.5)

n=1

AC

CE

Here qn is the Fourier coefficients of the function q(x). Let us suppose that fobs (x) =

∞ X

x

fobs n e 2ε sin(πnx).

n=1

Then putting t = T in (2.5) we obtain that fobs n = qn e−

Therefore qn = fobs n e 5

4ε2 π 2 n2 +1 T 4ε

4ε2 π 2 n2 +1 T 4ε

.

. (2.6)

ACCEPTED MANUSCRIPT

CR IP T

It is clear that if we have a noisy data the inverse problem (2.3)–(2.4) becomes unstable. The nonlinear analog of the problem (2.1) is:  2 ∂u ∂ u ∂u   x ∈ (0, 1), t ∈ (0, T ),  ε ∂x2 − ∂t = −u ∂x , (2.7) u(0, t) = u(1, t) = 0, t ∈ (0, T ),    u(x, 0) = q(x), x ∈ (0, 1).

AN US

Inverse problem consists of finding function q(x) by known additional information: u(x, T ) = fobs (x), x ∈ (0, 1). (2.8) For solution of the direct problem (2.7) we will apply the Florin-Hopf-Cole transformation [8, 9, 10] u(x, t) = 2ε

∂ ln ϕ(x, t). ∂x

ED

M

We suppose that ϕ(x, t) 6= 0. Therefore we obtain the following problem for the function ϕ: ( ϕt = εϕxx , x ∈ (0, 1), t ∈ (0, T ), ϕx (0, t) = ϕx (1, t) = 0, t ∈ (0, T ).

(2.9)

The particular solution of the problem (2.9) is

PT

ϕn (x, t) = b + e−επ

2 n2 t

cos(πnx).

CE

Here, the constant b > 1. Therefore the solution of the problem (2.7)–(2.8) can be written as un (x, t) =

2 2

πne−επ n t sin(πnx) . b + e−επ2 n2 t cos(πnx)

(2.10)

AC

The initial condition will be qn (x) = un (x, 0) =

πn sin(πnx) . b + cos(πnx)

We see that for any T > 0 the data un (x, t) → 0 for n → ∞ and the solution of inverse problem qn (x) → ∞ for n → ∞. Therefore inverse problem (2.7)– (2.8) is unstable. 6

ACCEPTED MANUSCRIPT

CR IP T

In both cases it follows from (2.5) and (2.10) that for small ε we can construct the stable algorithm to find the initial data. Much research and effort has gone into solving and analyzing Burgers equation. Beginning with the seminal work of Burgers [11], many analytic and numerical methods have been developed to investigate and solve different initial and initial-boundary value problems for Burgers equation. Fletcher has an excellent, if by now somewhat outdated, survey of the various methods and applications of Burgers equation [12]. More recent surveys of numerical methods can be found in [13, 14, 15, 16, 17, 18]. As far as theoretical results we mention [19, 20, 21, 22, 23, 24] and references there.

AN US

3. Numerical solution

Let us find the solution of the inverse problem (2.1)—(2.2) using the gradient method of minimizing of the cost functional J[q] =

Z1 

0

M

0

Z1  2 2 u(x, T ) − fobs (x) dx + α q(x) dx.

ED

Here α is the regularization parameter. In numerical experiments we will put α = 0 and use the iteration number s as the regularization parameter [25, 26, 7]. The algorithm of numerical solving of the inverse problem is following.

AC

CE

PT

1. Set s := 0 and q (0) (x) as an initial guess. We have to note that the choice of the function q (0) (x) must be done according to the conditions that are denoted by asymptotic analysis described in the following Section. This guarantee the existence of the solution for the direct problem with internal layer and allows to determine the front location correctly. 2. Find the solution of the direct problem  2 (s) (s) ∂u(s) ∂ u  (s) ∂u   ε − = −u + q (s) (x) u(s) , x ∈ (0, 1), t ∈ (0, T ],  ∂x2 ∂t ∂x u(s) (0, t) = ul (t), u(s) (1, t) = ur (t),     (s) u (x, 0) = u0 (x). 7

ACCEPTED MANUSCRIPT

4. Find the gradient of the functional 

(s)

ZT

u(s) (x, t)ψ (s) (x, t) dt + 2α q (s) (x).

0

AN US

 J q (x) = 0

t ∈ [0, T ),

CR IP T

3. Find the solution of the conjugate (adjoint) problem  (s) ∂ 2 ψ (s) ∂ψ (s)  (s) ∂ψ   ε = u + q (s) (x) ψ (s) , x ∈ (0, 1), +   ∂x2 ∂t ∂x ψ (s) (0, t) = 0, ψ (s) (1, t) = 0,       (s)  ψ (x, T ) = −2 u(s) (x, T ) − fobs (x) .

5. Find the approximate solution on the next iteration step   q (s+1) (x) = q (s) (x) − βs J 0 q (s) (x),

(3.1)

ED

M

where βs is the descent parameter. 6. Check the stopping criterion of the iterative process. If it is true, put q (s+1) (x) as the solution of the inverse problem. Otherwise, set s := s+1 and go to item 2. (a) In the case of noisy data fobs (kf¯−fobs k ≤ δ, f¯ — is exact  right-side  hand of the equation (2.4)) the stopping criterion is J q (s+1) ≤ δ 2 . (b) In the case of the exact data the iterative process stops when  J q (s+1) is less than the approximation error.

AC

CE

PT

Theoretical results about convergence of the iteration process (3.1) can be found in [27, 28, 29, 30]. The convergence of Landweber iteration for 2D coefficient inverse problem was proved in [31]. The solutions of these direct and conjugate problems have narrow boundary and interior layers (stationary or moving fronts) [1, 2]. Thus, these problems are extremely difficult for numerical treatment [3]. We remind that the main idea of our paper is to show how asymptotic analysis can help to solve this numerical problem. Here, asymptotic analysis is able to extract a priori information about solutions of both problems. This allows to construct so called dynamically adapted mesh that is able to reduce significantly the complexity of numerical calculations, save the computational resources and improve the numerical stability of solution of corresponding nonlinear singularly perturbed inverse problems. The main ideas of asymptotic analysis are described in the following Section. 8

ACCEPTED MANUSCRIPT

4. Asymptotic analysis

CR IP T

We illustrate the main ideas of using asymptotic analysis to obtain a priori information for our approach by considering the following class of problems  2 ∂u ∂ u ∂u    ε ∂x2 − ∂t = A(u, x) ∂x + B(u, x), x ∈ (0, 1), t ∈ (0, T ], (4.1) u(0, t) = ul , u(1, t) = ur ,    u(x, 0) = u0 (x),

M

AN US

where functions A(u, x) and B(u, x) are sufficiently smooth for (x, u) ∈ [0; 1] × (−∞; +∞). We are interested in solution of the front type which is close to two different levels in different parts of x ∈ [0; 1] with thin moving transition layer between these levels. The rigorous treatment of the problem (4.1) was done in [1]. It is important that under certain conditions asymptotic analysis allows to prove the existence of the solution (4.1), to study it’s properties and extract the following a priori information for the inverse problem solution: 1. speed and location of the transition layer; 2. width of the transition layer.

ED

The details of asymptotic analysis of (4.1) are developed in [1]. The main ideas of the asymptotic procedure are briefly represented below. According to [1] let consider two initial value problems:

PT

du + B(u, x) = 0, dx du A(u, x) + B(u, x) = 0, dx

CE

A(u, x)

u(0) = ul ,

x ∈ (0, 1],

(4.2)

u(1) = ur ,

x ∈ [0, 1).

(4.3)

AC

We suppose, that the solutions of problems (4.2) and (4.3) exist on x ∈ [0, 1], and denote them as ϕl (x) and ϕr (x) respectively. These functions define two levels and the internal layer arise between these levels. Assume, that the following inequalities are fulfilled for all x ∈ [0, 1]: ϕl (x) < ϕr (x),

(4.4)

A(ϕl (x), x)) > 0,

A(ϕr (x), x) < 0.

9

(4.5)

ACCEPTED MANUSCRIPT

AN US

CR IP T

These and some further assumptions were used in [1] for the construction of the formal asymptotic of the moving front type solution of (4.1) and proof of its existence. The problem of formation of the moving fronts is discussed in [32] and [33]. Particularly, for the reaction-diffusion problem the fast formation stage is described in [32]. Assume that the front is already formed at the moment t = 0 and the initial function uinit (x) has a thin transition layer between levels ϕl (x) and ϕr (x) located at the point x∗0 ∈ [0, 1]. Let define the moving front location by the function xt.p. (t, ε), that is the point of the intersection of the solution of (4.1) u(x, t, ε) and the level 1 ϕ(x) = (ϕl (x) + ϕr (x)). Note that this definition is conventionally: for 2 example, in [33] another definition for the front location is used. We put xt.p. (0, ε) = x∗0 and will describe the location and speed of the transition point in the form of power series of ε xt.p. (t, ε) = x0 (t) + εx1 (t) + . . . vt.p. (t, ε) = v0 (t) + εv1 (t) + . . . ,

(4.6) (4.7)

M

dxi . where vi (t) = dt Assume, that the Cauchy problem

ED

r (x) ϕR

x(0) = x∗0

(4.8)

PT

A(u, x)du dx ϕl (x) = r , dt ϕ (x) − ϕl (x)

CE

has solution x(t) such that x(t) ∈ (0, 1) for all t ∈ [0; T ],

AC

and the following inequalities are fulfilled ϕZr (x)

A(u, x)du > 0 for all x ∈ [0, 1],

(4.9)

ϕl (x)

Zs

ϕl (x)


(A(u, x) − V (x))du > 0 for all s ∈ ϕl (x), ϕr (x) , 10

(4.10)

ACCEPTED MANUSCRIPT

V (x) =

r (x) ϕR

A(u, x)du

ϕl (x)

ϕr (x) − ϕl (x)

.

(4.11)

CR IP T

where

AN US

Condition (4.9) guarantee that (4.1) has no stationary solutions, and (4.10) provides solvability of some equations in further asymptotic procedure. Note that both conditions (4.4) and (4.10) provides that the solution of (4.8) is strongly monotone which is important for further mesh construction. ¯ l and D ¯ r which are located at the left and right We define the domains D T T sides with respect to the point xt.p. (t, ε) respectively and will find the solution of (4.1) separately in these parts ( ¯l ul for (x, t) ∈ D T u= (4.12) ¯r . ur for (x, t) ∈ D T The functions ul and ur can be written in the form

ur = u¯r (x, ε) + Qr (ξ, t, ε) .

(4.13)

M

ul = u¯l (x, ε) + Ql (ξ, t, ε),

ED

Here u¯l,r (x, ε) – regular functions, which represent the solution far from transition point xt.p. (t, ε); function Ql,r (ξ, t, ε), where ξ = (x − xt.p. (t, ε))/ε, describe transition layer (moving front) near this point (ξ ≤ 0 for the function with the index l and ξ ≥ 0 for the function with the index r). The functions u¯l,r (x, ε) and Ql,r (ξ, t, ε) we write as power series of ε: (4.14)

l,r n l,r Ql,r (ξ, t, ε) = Ql,r 0 (ξ, t) + εQ1 (ξ, t) + . . . + ε Qn (ξ, t) + . . . .

(4.15)

CE

PT

n l,r u¯l,r (x, ε) = u¯l,r ul,r ¯n (x) + . . . , 0 (x) + ε¯ 1 (x) + . . . + ε u

AC

Terms of series (4.14) and (4.15) can be built by the standard procedure (see [1]) separately at two sides of the point xt.p. (t, ε). The main idea is that the functions ul and ur are joint with the continuous first derivatives (C (1) –matching conditions): u¯l (xt.p. (t, ε), ε) + Ql (0, t, ε) = = u¯r (xt.p. (t, ε), ε) + Qr (0, t, ε) = ϕ(xt.p. (t, ε)),

ε

(4.16)

∂Ql d¯ ur ∂Qr d¯ ul (xt.p. (t, ε), ε) + (0, t, ε) = ε (xt.p. (t, ε), ε) + (0, t, ε). (4.17) dx ∂ξ dx ∂ξ 11

ACCEPTED MANUSCRIPT

CR IP T

Substituting (4.14) into (4.1), at zero order of ε we obtain equations (4.2) and (4.3) . So we have u¯l0 (x) = ϕl (x), u¯r0 (x) = ϕr (x) and ( u¯l0 = ϕl (x) for x ∈ [0, xt.p. (t, ε)) (4.18) u¯0 (x, t) = u¯r0 = ϕr (x) for x ∈ (xt.p. (t, ε), 1] .

CE

PT

ED

M

AN US

Functions u¯l,r k (x), k ≥ 1 are the solutions of the linear Cauchy problems and can be written explicitly. In order to obtain the transition layer functions (4.15) we must do the change of variables (x, t) by (ξ, t) in (4.1), where ξ = (x − xt.p. (t, ε))/ε. Op∂2 ∂ 1 ∂2 1 ∂ ∂ erator ε 2 − takes the form + vt.p. (t, ε) − where vt.p. (t, ε) 2 ∂x ∂t ε ∂ξ ε ∂ξ ∂t ∂ 1 ∂ is defined in (4.6); operator changes to . ∂x ε ∂ξ Substituting these expressions and (4.15) into (4.1) and equating terms with different powers of ε, we obtain equations for functions Ql,r i (ξ, t) , i = 1, 2, . . .. For Ql,r 0 (ξ, t) we have the following problem:  2 l,r   ∂Ql,r ∂ Q0  l,r 0 l,r   − A(ϕ (x (t, ε)) + Q , x (t, ε)) − v (t, ε) = 0, t.p. t.p. t.p. 0  2  ∂ξ ∂ξ     Ql (0, t) + ϕl (xt.p. (t, ε), ε) = ϕ(xt.p. (t, ε)), 0 Q  r0 (0, t) + ϕr (xt.p. (t, ε), ε) = ϕ(xt.p. (t, ε)),      Ql0 (ξ, t) → 0 for ξ → −∞,    r Q0 (ξ, t) → 0 for ξ → +∞. (4.19) The exponential estimates for the solutions Ql,r (ξ, t) of (4.19) satisfy 0 ξ ≤ 0, ξ ≥ 0,

t ∈ [0, T ], t ∈ [0, T ],

(4.20) (4.21)

AC

|Ql0 (ξ, t)| ≤ Ceκξ , |Qr0 (ξ, t)| ≤ Ce−κξ ,

where C and κ are positive constants, were established in [1]. Using the continuous function ( l l ˜ xt.p. (t, ε), vt.p. (t, ε)) = ϕ (xt.p. (t, ε)) + Q0 (ξ, t) for ξ ≤ 0 Q(ξ, ϕr (xt.p. (t, ε)) + Qr0 (ξ, t) for ξ ≥ 0, 12

(4.22)

ACCEPTED MANUSCRIPT

ξ ∈ (−∞, +∞),

CR IP T

we can rewrite (4.19) as  2  ˜ ˜  ∂ Q   ˜ xt.p. (t, ε)) − vt.p. (t, ε) ∂ Q = 0,  − A( Q,  2  ∂ξ   ∂ξ ˜ xt.p. (t, ε), vt.p. (t, ε)) = ϕ(xt.p. (t, ε)), Q(0,   ˜  Q(−∞, xt.p. (t, ε), vt.p. (t, ε)) = ϕl (xt.p. (t, ε)),    ˜ Q(+∞, xt.p. (t, ε), vt.p. (t, ε)) = ϕr (xt.p. (t, ε)).

ϕrR(x0 )

AN US

(4.23) ˜ ∂Q give the folThe C (1) –matching conditions and explicit formulas for ∂ξ lowing problem (see (4.8)) for zero order term of the moving front location A(u, x0 )du dx0 ϕl (x0 ) = r , dt ϕ (x0 ) − ϕl (x0 )

x0 (0) = x∗0 .

(4.24)

ED

M

Note that for general cases the initial value problem (4.24) has to be solved numerically. In the example, considered in Section 5, it can be solved explicitly. Continuing this procedure the linear equations for functions Ql,r k (ξ, t), k ≥ l,r 1 can be obtained. Thus, functions Qk (ξ, t) can be written explicitly and satisfy the estimates

PT

|Qlk (ξ, t)| ≤ Ceκξ

|Qrk (ξ, t)| ≤ Ce−κξ

(ξ ≤ 0),

(4.25)

(ξ ≥ 0),

(4.26)

AC

CE

where C and κ are some positive constants. Using the C (1) –matching conditions (4.16)–(4.17) at first order of ε and ∂Ql,r l,r 1 , we obtain the problem for first order explicit formulas for Q1 (ξ, t) and ∂ξ term of the moving front position dx1 K(x0 (t)) = r · x1 + G1 (x0 (t), t), dt ϕ (x0 ) − ϕl (x0 )

x1 (0) = 0,

(4.27)

where K(x0 (t)) and G1 (x0 (t), t) are known functions (see [1]). The main terms of moving front speed v0 (t) + εv1 (t) are defined in (4.24) and (4.27). So the location of the moving front at two higher orders of asymptotic is defined as the solutions of the Cauchy problems (4.24) and (4.27). 13

ACCEPTED MANUSCRIPT

CR IP T

The width of the transition layer is defined by the estimates (4.20)–(4.21). It means that interior layer exponentially tends to functions ϕr,l (x) and its thickness is h = C|ε ln ε|. (4.28) Similar asymptotic analysis can be also done for the conjugate problem. Corresponding investigation shows that solution of the conjugate problem has two boundary layers near the points x = 0 and x = 1. The thickness of these layers is defined by (4.28). In the next Section we demonstrate effectiveness of the proposed approach by the particular case of problem (4.1).

AN US

5. Numerical schemes

ED

M

The small parameter at higher derivatives leads to narrow boundary and interior layers in the solutions of the direct and conjugate problems. For effective numerical solving of the inverse problem we propose to use special dynamically adapted meshes. These meshes must be refined at the areas with high gradients of the solutions of the direct and conjugate problems according the information which can be obtained from the asymptotic analysis described in . The main steps of dynamically adapted mesh constructing are presented below.

PT

5.1. Dynamically adapted mesh algorithm The main step is to determine the transition point location x(t) in a neighbourhood of which the mesh must be refined. It can be done by the following way using the results the previous Section.

AC

CE

1. Find the solution ϕl (x) and ϕr (x) of the equations (4.2) and (4.3), which in the case of problem (2.1) take the form dϕl (x) = q(x), dx dϕr (x) = q(x), dx

u(0) = ul ,

x ∈ (0, 1],

(5.1)

u(1) = ur ,

x ∈ [0, 1).

(5.2)

Note, that the function q(x) and boundary conditions ul , ur should satisfy to the conditions (4.4)–(4.5) and (4.9)–(4.10). In our numerical example we put ul = −8 and ur = 4. 14

ACCEPTED MANUSCRIPT

2. Find the speed of the front V (x) from the equation (4.11), which in the case of problem (2.1) take the form r (x) ϕR

V (x) =

ϕr (x) − ϕl (x)

.

3. Find the solution x(t) of the problem (4.24): dx = V (x), dt

x(0) = x∗0 .

CR IP T

(−u)du

ϕl (x)

(5.3)

(5.4)

AN US

4. Construct the mesh which has some refining in neighbourhoods of the transition point and the boundaries for each time layer.

In the following subsection we describe in details how to realize this algorithm numerically.

AC

CE

PT

ED

M

5.2. How to generate numerically the dynamically adapted mesh At the beginning we introduce uniform mesh XN0 only on x–dimension that has number of nodes N0 + 1 (that equals to N0 intervals): XN0 =
{xn , 0 ≤ n ≤ N0 : xn = 0 + (1 − 0)/N0 n}. This mesh will be used for calculation of the solution q(x). Because the solution q(x) does not have any features the number N0 could be quite small (in comparison with number of intervals that will be used for solving both direct and conjugate problems). Then it is possible to solve equations (5.1)–(5.2) and find numerically l ϕ (x) and ϕr (x) on the mesh XN0 . For numerical solving of (5.1)–(5.2) any suitable scheme can be applied. Then right-hand side (5.3) of the equation (5.4) can be calculated on the same mesh XN0 . After that we have to solve the problem (5.4). At first, we introduce uniform mesh TM0 only on t–dimension that has number of nodes M0 + 1 (that equals to M0 intervals): TM0 = {tm , 0 ≤ m ≤ M0 : tm = 0 + (T −
0)/M0 m}. We suggest to use Rosenbrock scheme with complex coefficient (CROS1), which is monotone, stable and has the order of accuracy O(τ 2 ) (see [34]):     V x(tm )   . (5.5) x(tm+1 ) = x(tm ) + (tm+1 − tm ) Re 1 − 1+i (t − t )V x(t ) m+1 m x m 2 15

ACCEPTED MANUSCRIPT

AN US

CR IP T

Remark: It is very important to mention that we have to perform interpolation of the tabulated function V (x) from basic mesh XN0 on each point x(tm ). It is need not to perform interpolation with order of accuracy greater than 3. Thus, we have dependence of transition point position x on time t as the solution of (5.4) (see Figure 1(a)). After that we are able to interpolate corresponding function t(x) on the uniform mesh X, which has intervals equal to |ε log ε|/C, C > 1 (where |ε log ε| is thickness of the interior layer). This interpolation is possible to perform because of strongly monotonic of the function x(t) (see Figure 1(b)). Therefor, we have grid TM which is quasi-uniform grid [35]. And finally, for each time node m of the mesh TM we are able to construct the piecewise uniform mesh XTm on x–dimension that has some refining in a neighborhood of the transition point x(tm ) and the bounds:

ED

M

XTm = {xn , 0 ≤ n ≤ N + Kint N + 2Kbound N :  xn = 0 for n = 0,      xn = xn−1 + hbound for n = 1, Kbound N ,    x = x for n = Kbound N + 1, Kbound N + Nlef t , n n−1 + hlef t  xn = xn−1 + hint for n = Kbound N + Nlef t + 1, Kbound N + Nlef t + Kint N ,      xn = xn−1 + hright for n = Kbound N + Nlef t + Kint N + 1, Kbound N + N + Kint N ,    xn = xn−1 + hbound for n = Kbound N + N + Kint N + 1, N + Kint N + 2Kbound N

PT

},

AC

CE

where N is a control parameter which determines the number of intervals allocated inside interior and boundary layers (intervals outside of the layers

16

ACCEPTED MANUSCRIPT

CR IP T

rationed out, so hlef t ∼ = hright ),       x(tm ) − Cint |ε ln ε| − 0 + Cbound |ε ln ε| N , Nlef t =  1 − 0 − 2Cint |ε ln ε| − 2Cbound |ε ln ε|

AN US

Nright = N − Nlef t ,     x(tm ) − Cint |ε ln ε| − 0 + Cbound |ε ln ε| hlef t = , Nlef t     1 − Cbound |ε ln ε| − x(tm ) + Cint |ε ln ε| hright = , Nright Cbound |ε ln ε| 2Cint |ε ln ε| hbound = , hint = , Kbound N Kint N

ED

M

Cint = 1, Cbound = 1 are control parameters that are able to adjust thickness of the interior and boundary layers respectively; Kint = 2, Kbound = 1 are control parameters that denote relative density of the refining inside interior and boundary layers respectively. As a result we have constructed the dynamically adapted mesh XTM (see for example Figure 2), which has on each time layer Tm number of intervals equal to N + Kint N + 2Kbound N . We want to remind that this mesh is ε-independent on spatial variable.

(5.6)

CE

PT

5.3. Solving of the forward problem For numerical computations we rewrite the problem (2.1) as  2 ∂u ∂ u ∂u    ε ∂x2 − ∂t = −u ∂x + u q(x), x ∈ (0, 1), t ∈ (0, T ], u(0, t) = ul , u(1, t) = ur ,    u(x, 0) = uinit (x),

AC

where ul (t) = −8, ur (t) = 4, and the initial condition uinit (x) has the following form:

x − 0.25 . ε For numerical solution of the system (5.6) we apply the stiff method of lines (SMOL) in order to reduce the PDE (5.6) to the system of ODEs that can be solved by Rosenbrock scheme with complex coefficient [34, 36, 37]. uinit (x) = (x2 − x − 2) − 6 tanh(−3ξ),

17

ξ=

ACCEPTED MANUSCRIPT

T

a)

t

t

x0

0

1

x

x0

x

AN US

0

b)

CR IP T

T

1

M

Figure 1: a) Solution x(t) of (5.4) for M0 = 6, N0 = 20. b) Quasiunidform mesh TM on t constructed by interpolating procedure.

AC

CE

t

PT

ED

T

0

x0

x

Figure 2: Example of the dynamically adapted mesh XTM for N = 4.

18

1

ACCEPTED MANUSCRIPT

AN US

This system can be rewritten as   du = f (u, q, t), dt  u(0) = uinit ,

CR IP T

At first, we apply finite-difference approximations of the derivatives with the second order of accuracy in (5.6). After that we obtain the following system of ODEs from which we should determine N − 1 unknowns functions un ≡ un (t) ≡ u(xn , t), n = 1, N − 1, (u0 and uN given as the boundary conditions):    dun 2 un+1 − un un − un−1 un+1 − un−1     dt = ε xn+1 − xn−1 xn+1 − xn − xn − xn−1 + un xn+1 − xn−1 − un qn , u0 (t) = −8, uN (t) = 4,     un (0) = uinit (xn ). t ∈ (0, T ],

(5.7)

PT

for n = 2, N − 2:

ED

M


T
T where u = u1 , u2 , u3 , . . . , uN −1 , q = q1 , q2 , q3 , . . . , qN −1 , f =
T
T f1 , f2 , f3 , . . . , fN −1 and uinit = uinit (x1 ), uinit (x2 ), uinit (x3 ), . . . , uinit (xN −1 ) . The vector-function f has the following structure. For n = 1:   u1 + 8 u2 − u1 u2 + 8 2ε − − u1 q1 ; + u1 f1 = x2 − x0 x2 − x 1 x1 − x0 x2 − x0

CE

2ε fn = xn+1 − xn−1



un+1 − un un − un−1 − xn+1 − xn xn − xn−1



+ un

un+1 − un−1 − un qn ; xn+1 − xn−1

AC

and for n = N − 1: fN −1

  2ε 4 − uN −1 uN −1 − uN −2 = − + xN − xN −2 xN − xN −1 xN −1 − xN −2 4 − uN −2 + uN −1 − uN −1 qN −1 . xN − xN −2

For numerical solution of this system of ODEs (5.7) we use the Rosenbrock scheme with complex coefficient (CROS1), which is monotone, stable and has 19

ACCEPTED MANUSCRIPT

the order of accuracy O(τ 2 ) (see [34]):

CR IP T

u (tm+1 ) = u (tm ) + (tm+1 − tm ) Re w , where w is a solution of the system    1+i (tm+1 − tm ) fu u (tm ), q (tm ), tm w = I− 2  tm + tm+1  = f u (tm ), q (tm ), . 2

AN US

Here I is the identity matrix, and the components of the Jacobian matrix fu are given by: for n = 1:   2ε 1 1 u2 + 8 ∂f1 = − − + − q1 , ∂u1 x2 − x0 x2 − x1 x1 − x0 x2 − x0   2ε 1 u1 ∂f1 = ; + ∂u2 x2 − x0 x2 − x1 x2 − x0 for n = 2, N − 2:

PT

ED

M

  2ε 1 ∂fn un = , − ∂un−1 xn+1 − xn−1 xn − xn−1 xn+1 − xn−1   ∂fn 2ε 1 1 un+1 − un−1 = − − qn , − + ∂un xn+1 − xn−1 x − xn xn − xn−1 xn+1 − xn−1  n+1  ∂fn 2ε 1 un = ; + ∂un+1 xn+1 − xn−1 xn+1 − xn xn+1 − xn−1

and for n = N − 1:

AC

CE

  ∂fN −1 2ε 1 uN −1 = − , ∂uN −2 xN − xN −2 xN −1 − xN −2 xN − xN −2   ∂fN −1 2ε 1 1 4 − uN −2 = − − + − qN −1 . ∂uN −1 xN − xN −2 xN − xN −1 xN −1 − xN −2 xN − xN −2

The other components of fu for the considered equation equal to zero. Note, that in the numerical scheme we use for calculations u (tm ) that has been determined on the XTm mesh. So, after applying the Rosenbrock scheme, we have to interpolate u (tm+1 ) on XTm+1 mesh. We note that q (tm ) should be also taken on XTm mesh. 20

ACCEPTED MANUSCRIPT

obs

t ∈ [0, T ),

CR IP T

5.4. Solving of the conjugate problem Now we rewrite the conjugate problem  ∂ 2 ψ ∂ψ ∂ψ    ε =u + ψ q(x), x ∈ (0, 1), +  2  ∂x ∂t ∂x ψ(0, t) = 0, ψ(1, t) = 0,        ψ(x, T ) = −2 u(x, T ) − f (x) ,

(5.8)

M

AN US

where u(x, t) is the known solution of the direct problem (5.6), fobs (x) — input data of the inverse problem (2.1)–(2.2). For numerical solution of the system (5.8) we also apply the stiff method of lines (SMOL) in order to reduce the PDE (5.8) to the system of ODEs from which we should determine N −1 unknowns functions ψn ≡ ψn (t) ≡ ψ(xn , t), n = 1, N − 1, (ψ0 and ψN given as the boundary conditions):    dψn 2 ψn+1 − ψn ψn − ψn−1 ψn+1 − ψn−1   = −ε − + u(xn , t) + ψn qn ,   xn+1 − xn−1 xn+1 − xn xn − xn−1 xn+1 − xn−1  dt ψ0 (t) = 0, ψN (t) = 0,        ψn (T ) = −2 u(xn , T ) − fobs (xn ) .

PT

ED

This system can be rewritten as  dΨ   = f (Ψ, u , q , t), t ∈ [0, T ), dt     Ψ(T ) = −2 u (T ) − fobs ,

(5.9)

AC

CE


T
T where Ψ = (ψ1 , ψ2 , ψ3 , . . . , ψN −1 , q = q1 , q2 , q3 , . . . , qN −1 , f =
T
T f1 , f2 , f3 , . . . , fN −1 , u (t) = u1 (t), u2 (t), u3 (t), . . . , uN −1 (t) and
T fobs = fobs (x1 ), fobs (x2 ), fobs (x3 ), . . . , fobs (xN −1 ) (note, that fobs should be preliminarily interpolated on XTM mesh). The vector-function f in the case of conjugate problem has the following structure. For n = 1:   ψ2 − ψ1 ψ1 − 0 ψ2 − 0 2ε − + u1 + ψ1 q1 ; f1 = − x2 − x0 x2 − x1 x1 − x0 x2 − x0 21

ACCEPTED MANUSCRIPT

2ε fn = − xn+1 − xn−1



ψn+1 − ψn ψn − ψn−1 − xn+1 − xn xn − xn−1



+ un

ψn+1 − ψn−1 + ψn qn ; xn+1 − xn−1

CR IP T

for n = 2, N − 2:

and for n = N − 1:

  2ε 0 − ψN −1 ψN −1 − ψN −2 fN −1 = − − + xN − xN −2 xN − xN −1 xN −1 − xN −2 0 − ψN −2 + uN −1 + ψN −1 qN −1 . xN − xN −2

AN US

For numerical solution of (5.9) we use the Rosenbrock scheme with complex coefficient (CROS1) in the following form:

M

Ψ (tm−1 ) = Ψ (tm ) + (tm−1 − tm ) Re w , where w is a solution of the SLAE    1+i (tm−1 − tm ) fΨ Ψ (tm ), u (tm ), q (tm ), tm w = I− 2  tm−1 + tm  = f Ψ (tm ), u (tm ), q (tm ), . 2

CE

PT

ED

Here I is the identity matrix, and the components of the Jacobian matrix fΨ are given by: where for n = 1:   2ε 1 1 ∂f1 =− − − + q1 , ∂ψ1 x2 − x0 x2 − x1 x1 − x0   ∂f1 2ε 1 1 =− + u1 ; ∂ψ2 x2 − x0 x2 − x1 x2 − x0

AC

for n = 2, N − 2:

  ∂fn 2ε 1 1 =− − un , ∂ψn−1 xn+1 − xn−1 xn − xn−1 xn+1 − xn−1   ∂fn 2ε 1 1 =− − − + qn , ∂ψn xn+1 − xn−1 xn+1 − xn xn − xn−1   ∂fn 2ε 1 1 + un ; =− ∂ψn+1 xn+1 − xn−1 xn+1 − xn xn+1 − xn−1 22

ACCEPTED MANUSCRIPT

∂fN −1 =− ∂ψN −2 xN ∂fN −1 =− ∂ψN −1 xN

  2ε 1 1 − uN −1 , − xN −2 xN −1 − xN −2 xN − xN −2   2ε 1 1 − − + qN −1 . − xN −2 xN − xN −1 xN −1 − xN −2

CR IP T

and for n = N − 1:

AN US

The other components of fΨ for the considered equation equal to zero. Note, that in the numerical scheme we use for calculations Ψ (tm ) that has been determined on the XTm mesh. So, after applying the Rosenbrock scheme, we have to interpolate Ψ (tm−1 ) on XTm−1 mesh. We also have to remind that q (tm ) should be also taken on XTm mesh. 6. Numerical example

CE

PT

ED

M

Let us consider some examples of numerical calculations in order to demonstrate the efficiency of the proposed algorithm. The examples of solving of the direct problem (see Figure 3 and corresponding video-file attached to this paper) and the conjugate problem (see Figure 4 and corresponding video-file attached to this paper) have been calculated for the following set of parameters: ε = 10−1.5 ' 0.03, Cint = 0.5, Cbound = 0.5, 2N = N0 = 40, Kint = 2, Kbound = 1, T = 0.2, M = 100, βs = 0.03. We use interpolation for simulated data in the solution of the conjugate problem. Solution of corresponding inverse problem for the model function q(x) = sin(3πx) and initial approximation q (0) (x) ≡ 0 is represented on Figures 5 and 6 (corresponding iterative processes are also represented in the videofiles attached to this paper) for different noise-levels δ which are typical for corresponding practical applications. 7. Conclusion

AC

Using the information based on the rigorous asymptotic analysis we propose the effective numerical algorithm for solving of the coefficient inverse problem for singularly perturbed parabolic equation. The narrow boundary and interior layers which appears in the solutions of singularly perturbed reaction-diffusion-advection equations extremely complicate the process of the numerical solving of the inverse problems for such type of equations and often require to use special meshes. Otherwise, the 23

ACCEPTED MANUSCRIPT

4

CR IP T

2 0 uinitial

t = 0.09

u −2 −4

AN US

−6 −8 0

0.2

0.4

x

0.6

0.8

1

M

Figure 3: Solution of the direct problem u(x, t) for q(x) = 2x − 1 + 2 sin(5πx) + 0.35.

ED

4 2

PT

0

ψinitial

CE

ψ −2 −4

t = 0.09

AC

−6 −8

0

0.2

0.4

x

0.6

0.8

1

Figure 4: Solution of the conjugate problem u(x, t) for q(x) = 2x − 1 + 2 sin(5πx) + 0.35.

24

ACCEPTED MANUSCRIPT

1.5

CR IP T

1

0.5

q

0

AN US

−0.5

−1

−1.5

0

0.2

0.4

x

0.6

0.8

1

M

Figure 5: Model function q(x) has been restored for δ ' 0.01%: dotted line — model function, squared line — approximate solution. 1.5

ED

1

PT

0.5

0

CE

q

AC

−0.5

−1

−1.5

0

0.2

0.4

x

0.6

0.8

1

Figure 6: Model function q(x) has been restored for δ ' 0.1%: dotted line — model function, squared line — approximate solution.

25

ACCEPTED MANUSCRIPT

AN US

CR IP T

small parameter allows to get a priori information about the exact solution, which gives the possibility for productive combinations of asymptotic and numerical approaches for special mesh construction. The dynamically adapted meshes which we use for the numerical solving are built on the basis of the asymptotic analysis and the main properties of the solutions (boundary and interior layers) of direct and conjugate problems. It allows to speed up the numerical calculations, significantly reduce the instability of the numerical solution and substantially save computing resources. It is important to mention that the proposed method can be applied for wide class of the nonlinear inverse problems for singularly perturbed reactiondiffusion-advection equations. The authors kindly thank all the reviewers for all remarks and comments which helped to improve the article. References

M

[1] E. A. Antipov, N. T. Levashova, N. N. Nefedov, Asymptotics of the front motion in the reaction-diffusion-advection problem, Computational Mathematics and Mathematical Physics 54 (10) (2014) 1536– 1549.

ED

[2] V. T. Volkov, N. N. Nefedov, E. A. Antipov, Asymptotic-numerical investigation of generation and motion of fronts in phase transition models, Lecture Notes in Computer Science 9045 (2015) 408–416.

CE

PT

[3] D. V. Lukyanenko, V. T. Volkov, N. N. Nefedov, L. Recke, K. Schneider, Analytic-numerical approach to solving singularly perturbed parabolic equations with the use of dynamic adapted meshes, Modeling and Analysis of Information Systems 23 (3) (2016) 334–341.

AC

[4] S. I. Kabanikhin, M. A. Shishlenin, Quasi-solution in inverse coefficient problems, Journal of Inverse and Ill-Posed Problems 16 (7) (2008) 705– 713. [5] S. I. Kabanikhin, M. A. Shishlenin, Using the a priori information in the coefficient inverse problems for hyperbolic equations, Proceedings of Institute of mathematics and mechanics UrB RAS 18 (1) (2012) 147– 164. 26

ACCEPTED MANUSCRIPT

CR IP T

[6] S. I. Kabanikhin, N. S. Novikov, I. V. Oseledets, M. A. Shishlenin, Fast toeplitz linear system inversion for solving two-dimensional acoustic inverse problem, Journal of Inverse and Ill-Posed Problems 23 (6) (2015) 687–700. [7] S. I. Kabanikhin, Inverse and Ill-posed Problems Theory and Applications, de Gruyter, 2011.

AN US

[8] V. A. Florin, Some the simplest nonlinear problems of consolidation of water-saturated earth environment, Proceedings of the Academy of Sciences of the USSR. Department of technical Sciences 9 (1948) 1389– 1402.

[9] E. Hopf, The partial differential equation ut + uux = µuxx , Comm. Pure Appl. Maths. 3 (1950) 201–230. [10] J. D. Cole, On a quasilinear parabolic equation occurring in aerodynamics, Quart. Appl. Maths. 9 (1951) 225–236.

M

[11] J. M. Burgers, A mathematical model illustrating the theory of turbulence, Advances in Applied Mechanics 1 (1948) 171–199.

ED

[12] C. A. J. Fletcher, Burgers equation: A model for all reasons, in: J. Noye (Ed.), Numerical Solution of Partial Differential Equations, North-Holland, Amsterdam, 1982, pp. 139–225.

PT

[13] M. M. Cecchi, R. Nociforo, P. Patuzzo, Space-time finite elements numerical solutions of burgers problems, Le Matematiche LI (1996) 43–57.

CE

[14] J. C. Tannehill, D. A. Anderson, R. H. Pletcher, Computational Fluid Mechanics and Heat Transfer, 2nd Edition, Taylor & Francis, Washington, DC, 1997.

AC

¨ [15] S. Kutluay, A. R. Bahadir, A. Ozde¸ s, Numerical solution of onedimensional burgers equation: explicit and exact-explicit finite difference methods, Journal of Computational and Applied Mathematics 103 (2) (1999) 251–261. [16] A. N. Boules, On the least-squares conjugate-gradient solution of the finite element approximation of burgers’ equation, Applied Mathematical Modelling 25 (9) (2001) 731–741. 27

ACCEPTED MANUSCRIPT

[17] I. A. Hassanien, A. A. Salama, H. A. Hosham, Fourth-order finite difference method for solving burgers equation, Applied Mathematics and Computation 170 (2) (2005) 781–800.

CR IP T

[18] S. Dhawan, S. Kapoor, S. Kumar, S. Rawat, Contemporary review of techniques for the solution of nonlinear burgers equation, Journal of Computational Science 3 (5) (2012) 405–419. [19] A. N. Boules, On the existence of the solution of burgers equation for n < 4, Internat. J. Math. & Math. Sci. 13 (4) (1990) 645–650.

AN US

[20] M. M. Cecchi, R. Nociforo, P. Patuzzo, Burgers problems — theoretical results, Rivista di Matematica Pura ed Applicata 21 (1997) 159–174. [21] H.-O. Kreiss, J. Lorenz, Initial-Boundary Value Problems and the Navier-Stokes Equations, Classics in Applied Mathematics, SIAM, 2004.

[22] W. L. Wood, An exact solution for burgers equation, Communications In Numerical Methods in Engineering 22 (7) (2006) 797–798.

ED

M

[23] H. R. Clark, M. A. Rincon, A. Silva, Analysis and numerical simulation of viscous burgers equation, Numerical Functional Analysis and Optimization 32 (7) (2011) 695–716.

PT

[24] Y. Benia, B.-K. Sadallah, Existence of solutions to burgers equations in domains that can be transformed into rectangles, Electronic Journal of Differential Equations 157 (2016) 1–13. [25] O. M. Alifanov, E. A. Artuhin, S. V. Rumyantsev, Extreme methods for the solution of ill-posed problems, Nauka, Moscow, 1988.

CE

[26] O. M. Alifanov, Inverse Heat Transfer Problems, International Series in Heat and Mass Transfer, Springer–Verlag, Berlin, 1994.

AC

[27] H. W. Engl, M. Hanke, O. Scherzer, A convergence analysis of the landweber iteration for nonlinear ill-posed problems, Numer. Math. 72 (1995) 21–37. [28] H. W. Engl, O. Scherzer, A convergent rate result for a steepest descent method and a minimal error method for the solution of nonlinear illposed problems, ZAA 14 (1995) 369–377. 28

ACCEPTED MANUSCRIPT

[29] V. V. Vasin, On the convergence of gradient-type methods for nonlinear equations, Doklady Mathematics 57 (2) (1998) 173–175.

CR IP T

[30] V. V. Vasin, I. I. Eremin, Operators and iterative processes of Fejer type. Theory and applications, 2005.

[31] S. I. Kabanikhin, O. Scherzer, M. A. Shishlenin, Iteration methods for solving a two dimensional inverse problem for a hyperbolic equation, J. Inverse and Ill-Posed Problems 11 (1) (2003) 87–109.

AN US

[32] V. T. Volkov, N. E. Grachev, N. N. Nefedov, A. N. Nikolaev, On the formation of sharp transition layers in two-dimensional reaction-diffusion models, Computational Mathematics and Mathematical Physics 47 (8) (2007) 1301–1309.

[33] V. T. Volkov, N. N. Nefedov, Asymptotic-numerical investigation of generation and motion of fronts in phase transition models, Lecture Notes in Computer Science 8236 (2013) 524–531.

ED

M

[34] A. B. Alshin, E. A. Alshina, N. N. Kalitkin, A. B. Koryagina, Rosenbrock schemes with complex coefficients for stiff and differential algebraic systems, Computational Mathematics and Mathematical Physics 46 (8) (2006) 1320–1340. [35] N. N. Kalitkin, A. B. Alshin, E. A. Alshina, B. V. Rogov, Computations on Quasi-Uniform Grids [in Russian], Fizmatlit, Moscow, 2005.

PT

[36] E. Hairer, G. Wanner, Solving of Ordinary Differential Equations. Stiff and Differential-Algebraic Problems, Springer, 2002.

AC

CE

[37] H. H. Rosenbrock, Some general implicit processes for the numerical solution of differential equations, Comput. J. 5 (4) (1963) 329330.

29