Numerical Solution of the Retrospective Inverse Parabolic Problem on Disjoint Intervals

: The retrospective inverse problem for evolution equations is formulated as the reconstruction of unknown initial data by a given solution at the ﬁnal time. We consider the inverse retrospective problem for a one-dimensional parabolic equation in two disconnected intervals with weak solutions in weighted Sobolev spaces. The two solutions are connected with nonstandard interface conditions, and thus this problem is solved in the whole spatial region. Such a problem, as with other inverse problems, is ill-posed, and for its numerical solution, speciﬁc techniques have to be used. The direct problem is ﬁrst discretized by a difference scheme which provides a second order of approximation in space. For the resulting ordinary differential equation system, the positive coerciveness is established. Next, we develop an iterative conjugate gradient method to solve the ill-posed systems of the difference equations, which are obtained after weighted time discretization, of the inverse problem. Test examples with noisy input data are discussed.


Introduction
Interface problems have many applications in biology, applied mechanics, heat and mass transfer, etc. [1][2][3][4][5].Regularity of linear inhomogeneous parabolic interface transmission problems was investigated in [1].In [2], the author considered two-phase parabolic free boundary problems and established a monotonicity formula for heat functions in disjoint domains.The existence and uniqueness of strong solutions for linear parabolic partial differential equations in two adjoining domains connected through nonlinear Neumanntype interface conditions was studied in [4].In [5], the advantages of exact representation of the solution in the interface in numerical schemes is discussed.
In this work, we consider a parabolic transmission problem on disjoint domains with exact interface conditions.The derivation of this conditions is discussed in [6,7].
Problems in disjoint domains can be considered a specific instance of interface problems.Mathematically, interface problems result in partial differential equations featuring discontinuities in both the input data and solutions across one or more hypersurfaces with dimensions lower than the domain in which the problem is defined.Different types of conjugation conditions that link domain solutions and their derivatives are known.
The direct parabolic transmission problem in disjoint domains is well studied in the literature.The existence and uniqueness of the weak and strong solutions and the a priori estimate in an appropriate Sobolev-like space were proven in [6,[8][9][10][11].Numerical methods for solving such problems were constructed and analyzed in [8][9][10][11].
Due to many applications in physics, mechanics, biology, finance, etc., a large number of results for inverse problems was obtained in recent decades [12][13][14][15][16][17][18][19][20][21].As examples of inverse problem applications, we can mention the reconstruction of thermal sources, intensity estimation in the heat transfer, as well as the initial condition estimation based on final time temperature measurements.The backward problem of determination of the initial condition of the direct (forward) problem from a known final time value of the parabolic problem solution is often referred to as a retrospective inverse problem (RIP).
Implementing an inverse source problem for a parabolic transmission problem by using a single measurement on a part of a time interval was considered in [22].In [23], the authors constructed and analyzed a numerical method for simultaneously determining the initial value and source in a parabolic transmission problem.In [24,25], an inverse boundary value problem for the heat equation in a two-layered medium with unknown inclusions is solved numerically.
In our previous works, we constructed numerical algorithms for solving the parabolic transmission problem in disjoint domains.The inverse problem for reconstruction of a time-dependent source in classical and time-fractional problems from point or integral observations was solved in [26,27], respectively.The inverse problem for identification of external boundary conditions from point measurements was developed and analyzed in [28].In all these papers, using implicit-explicit time stepping for each time level, the full inverse problem was decoupled into two Dirichlet inverse problems.Then, a decomposition technique or loaded equation method was applied.
The retrospective inverse problem for parabolic equation consists of reconstruction of the unknown initial condition from a given final time observation.It is not well posed, as a small perturbation of the input data can produce large perturbations in the solution (see, for example, [29][30][31][32][33][34][35]).The main difficulty in the construction of a numerical approximation of the solution comes from the strong ill-posedness of the differential problem and the poor conditioning of the corresponding algebraic equations.
To the best of our knowledge, [36] is the first work that studied an RIP numerically.The RIP for heat conduction was solved as an optimal control problem of an object with distributed parameters in [29].In [31,33,34], the authors constructed an iterative numerical method for the solution of an RIP for one-and multi-dimensional parabolic equations.An efficient numerical method for solving the RIP for a time-fractional parabolic equation was developed in [32].The authors applied a conjugate gradient-type regularization method to solve the discrete ill-posed linear system.In [37,38], the existence and uniqueness of a quasi-solution to backward time-fractional and space-time-fractional diffusion equations were studied.The Levenberg-Marquardt regularization method was applied to solve the ill-posed problem.A regularization approach for solving the retrospective diffusion problem was also proposed in [30].In [35], the variation method was utilized for solution of the RIP for a nonlinear, heterogeneous Burgers' equation.Another important method, at least in terms of theoretical aspects, is the method of regularization of a parabolic equation backward in time through nonlocal initial boundary value problems (see [16] (Section 3.3)) [39].A disdvantage of this method is the nonlocal initial conditions which arise that require specific numerical techniques.
However, there are no results for the inverse retrospective problem of our type in the case where the governing equations are situated in disconnected domains.
In this paper, we unfold the finite-difference method proposed in [31,33,34] and solve numerically the nonstandard RIP for recovering the initial condition in a parabolic equation in disjoint domains.
The remaining part of this paper is organized as follows.In the next section, we introduce the physical model problem.Section 3 is devoted to the well-posedness of the direct problem in specific weighted Sobolev spaces.In Section 4, we describe the semidiscretization of the direct problem.In Section 5, we describe the iterative method for solving the inverse problem in Equations ( 1)-( 5) and ( 7)- (9).The stability of the numerical approach is discussed in Section 6.In the next section, we propose our computational results.This paper is finalized in the Conclusions section.

Model Problem
In this section, we briefly describe the direct and inverse problems for our physical model.
The retrospective inverse problem for the system in Equations ( 1)-( 6) consists of the recovery of an a priory unknown initial state (Equation ( 6)) and the solution u(x, t) = {u 1 (x, t), u 2 (x, t)}, 0 ≤ t ≤ T with the known observations at the final time: The sense of the simple model, described by the system in Equations ( 1)-( 9), is heat conduction into interacted one-dimensional rods of lengths b j − a j (j = 1, 2) whose temperature at point x and time t is modeled by the functions u j (x, t) (j = 1, 2), which solves Equation (1).The left end of the first rod and the right end of the second rod are isolated if ϕ 1 (t) = ϕ 2 (t) = 0 in Equatons (4) and ( 5) and retain the temperature otherwise.Next, there is a heat exchange between the rods described by the conditions in Equations ( 2) and ( 3).The initial temperature (Equation ( 6)) of the rods is unknown.It is required to reconstruct the functions u • j (j = 1, 2) from the measured temperatures φ j (x) (j = 1, 2) at t = T.

Well-Posedness of the Direct Problem
In this section, we present the results for the existence and uniqueness of a strong solution to the problem in Equations ( 1)-( 6) as obtained in [8].However, by assuming more smoothness (Equation ( 7)), in comparison with [8], we obtain a smooth solution u = {u 1 , u 2 } to the differential problem in Equations ( 1)- (6).Our main concern is the coerciveness of the bilinear form A. This property will be used in the construction of the iterative algorithm.
We consider the product space equipped with the inner product and norm In addition, we introduce the space equipped with the following inner product and norm: where

Specifically, we assign
We introduce the bilinear form corresponding to the problem in Equations ( 1)-( 6): On the base of the results in [8,9], we establish the following statements: Lemma 1.Under the conditions in Equations( 7) and (8), the bilinear form A, defined by Equation (10), is symmetric and bounded on H 2 × H 2 .Moreover, this form is also coercive on H 2 0 (i.e., there exists a constant c Let u(t) be a function mapping Ω ∈ R n into a Hilbert space H. Furthermore, we define [40] the Sobolev space H k (Ω, H) with an inner product where k in a non-integer.For k = 0, we set L 2 (Ω, H) = H 0 (Ω, H).
Next, we introduce the space Theorem 1.Let the assumptions in Equations ( 7) and (8) hold and u 0 = (u 10 , u 20 ) ∈ H 1 0 .Then, the initial boundary value problem in Equations (1)-( 6) has a unique strong solution u = (u 1 , u 2 ) ∈ H 2,1 , and the following a priory estimate holds true Let us note that the increase in the smoothness of the input data leads to the increase in the solution's smoothness.Furthermore, using this fact, we assume that the sufficient smoothness of the differential problem's solution requires the construction of corresponding difference schemes.

Semidiscrerization of the Direct Problem
First, we a construct finite difference scheme for the direct problem in Equations ( 1)-( 6).We introduce a uniform partition of the intervals Ω j through meshes ω h j (j = 1, 2): We also denote the mesh function y j (x, t) at node (x j,i j , t) by y i,j i (t) = y i,j i .Next, we also introduce the notations h j y j,i j w j,i j , y = y j , y j .
By applying the finite-volume method and eliminating the external boundary grid nodes x 1,0 = a 1 , x 2,N 2 = b 2 , we obtain the spatial semidiscretization of Equations ( 1)-(6): or i j = 2 − j, 3 − j, . . ., N j − j + 1, where L j u j,i j (t) = L j u j,i j = q j,i j u j,i j Furthermore, for convenience, we rearrange the indexes, representing the unknown solution and corresponding functions in the form With these notations, the full discrete scheme in Equation ( 11) is rewritten as follows: where and In this framework, the scalar product and the corresponding norm are defined as follows: Considering the main goal of the present research, namely the numerical solution of the retrospective problem in Equations ( 1)-( 5), ( 7) and ( 9), the following results are of basic importance: smooth in the space solution to the problem in Equations (1)-( 8) and {U i (t)} be the solution to the discrete problem in Equation (12).Then, the operator L is positively definite, and the difference scheme is of the second order of approximation.
Proof.The second order of approximation directly follows from the construction of Equation (11).
Next, we multiply the first N 1 number of equations in Equation ( 12) by β 2 and the remaining N 2 number of equations by β 1 , and in the resulting discrete system, we set Therefore, under the conditions in Equation ( 8), we obtain where C is a positive constant.
Let us note that from the definiteness (i.e., the coerciveness of the operator L), the solution to the energy stability for the ODEs (Equation ( 11)) and its global existence directly follow.

Iterative Solution to the Inverse Problem
Now, we consider the uniform temporal mesh ω τ = {t n : t n = nτ, n = 0, 1, . . ., M, τ = T/M}, such that the computational domain Q j T = Ω j × [0, T] is discretized by a rectangular grid ω h j × ω τ (j = 1, 2).The mesh function y j (x, t) at node (x j,i j , t n ) is denoted by y n i,j i .Furthermore, the σ-weighted time stepping (0 < σ ≤ 1) leads to the following full discrete scheme: where L n+σ j (y) = L j y, σt n+1 + (1 − σ)t n , g n+σ j,i j = g j,i j (t n+1 + (1 − σ)t n ).As before, we rearrange the indexes, where the unknown solution and corresponding functions are With these notations, the full discrete scheme (Equation ( 13)) is rewritten as follows: where In order to solve the inverse problem in Equations ( 1)-( 5) and ( 9) for reconstruction of the initial conditions (Equation ( 6)), we consider the discretization in Equation ( 13) (or Equation ( 14)) associated with Equation (9), namely Let E be the unit matrix of a size N × N and introduce the notations Thus, from Equation ( 14), for the solution vector Consequently, after time integration, from Equation ( 15), we obtain where Therefore, for n = M − 1, from Equation ( 16) we obtain where Finally, the initial condition is determined by solving the linear equation Since the operator S n+σ , n = 0, 1, . . ., M − 1 is self-adjoint, and S n+σ = (S n+σ ) * and is positive, it follows from Lemma 2 that its inverse exists and is also self-adjoint and positive.Consequently, for the solution of the system of linear algebraic Equation ( 18) with a full coefficient matrix, it is appropriate to use the conjugate gradient method (CGM).The steps are described in Algorithm 1.
Algorithm 1 CGM for solving an RIP.
Require: φ j (x), j = 1, 2, initial guess V 0 for U 0 , consistent with external boundary conditions, accuracy ε > 0 Ensure: S k+σ P k by solving the problem The main advantage of the proposed iterative method for solving an ill-posed discrete linear system for the inverse problem is that in contrast to the often-used regularization methods for solving the corresponding minimization problems, the use of uncertainties such as a regularization parameter or quadratic norm is avoided.Instead, in the conjugate gradient iterative procedure, the regularization is performed during the iterations, since the number of iterations serves as a regularization parameter.
Furthermore, the iterations execute quickly, since both the direct problem and the problem for identifying Z k have the same coefficient matrix.

Stability Discussion
Using the standard theory of difference schemes (see, for example, Chapter 7 in [43]), one can prove easily the following assertion: Theorem 2. Let {U n } be the solution to the weighted difference scheme in Equation (13).Then, there exists a constant C > 0 independent of τ, h 1 and h 2 , and τ 0 > 0 such that for all 0 < τ ≤ τ 0 , we have Regarding the direct problem in Equations ( 1)-( 6), solved by the difference scheme in Equation ( 14), from the above considerations, when taking into account that L n+σ i is self-adjoint, positively defined and bounded, in light of Chapter 7 in [43], we deduce that the numerical scheme is stable if and the following a priory estimate holds: Therefore, the solution to the finite difference scheme in Equation ( 14) is stable for σ ≥ 1/2 and convergent with the order O(τ ) otherwise.For σ < 1/2, the numerical discretization in Equation ( 14) is stable for a sufficiently small time step [33].Now, we consider the RIP.

Numerical Tests
In this section, some examples will be presented to illustrate the algorithm performance.Let The right-hand sides f j (x, t), γ j (t) and initial conditions u • j (x) (j = 1, 2) are determined such that u = (u 1 , u 2 ), u 1 (x, t) = e −t/2 sin(πx/2), u 2 (x, t) = e −t/4 cos(πx/2) will be the exact solution to the direct problem in Equations ( 1)-( 6) and (20).We will refer to this solution as the true (exact) solution.All computations were performed for Example 1. (Direct problem) First, we test the accuracy of the numerical solution u n j,i j , j = 1, 2, i j = 0, 1, . . ., N j , n = 0, 1, . . ., M for the direct problem in Equations (1)-( 6), computed by Equation (13) or (14).In Tables 1 and 2, we give the errors (E N j j ) and orders of convergence (CR j ) at the final time , for σ = 1/2 and σ = 1.In the case where σ = 1/2, we fixed the ratio to τ/h = 1, while for σ = 1, we had τ = h 2 .The computations show that the accuracy of the numerical solution of the direct problem computed by Equation (14) was O(τ 2 + h 2 ) for σ = 1/2, and it was O(τ + h 2 ) otherwise.In Figures 1 and 2, we plotted the exact (true) and numerical solutions to the direct problem and the corresponding error u j (x, T) − u T j,i j at the final time for σ = 1/2 and σ = 1, N 1 = 80 and τ = h.Exact solution u 1 (x, T), u 2 (x, T) (solid line) and numerical solution to the direct problem u T 1,i 1 , u T 2,i 2 (line with circles) at the final time (left) and the corresponding error (right), where σ = 1 and τ = h (Example 1).

Example 2. (Inverse problem)
We demonstrate the efficiency of Algorithm 1 in recovering the initial condition u • 1 , u • 2 and solution u 1 , u 2 .We consider perturbed measurements φ j (x j,i j ) = u j (x j,i j , T) + 2ρ j u j (x j,i j , T)( j (x j ) − 0.5), i j = 1, 2, . . ., N j − 1, j = 1, 2, where j (x j ) is a random function uniformly distributed in the interval [0, 1], where ρ j is the amplitude and u j (x j,i j , T) is the true solution.
In Figures 3 and 4, we depict the true function u ) and recovered initial function u 0 j = (u 0 1,i 1 , u 0 2,i 2 ) and the corresponding error u • j (x) − u 0 j,i j for σ = 1/2 and σ = 1, respectively.In Figures 5 and 6, we plot the true solution u(x, T) = (u 1 (x, T), u 2 (x, T)) and recovered solution u M j = (u M 1,i 1 , u M 2,i 2 ) and the corresponding error u T − u T j for σ = 1/2 and σ = 1, respectively.We observed better precision from Algorithm 1 for σ = 1 in comparison with σ = 1/2, especially for the solution at the final time.The precision was almost the same as that for the numerical solution to the direct problem.
In Figures 7 and 8, we plotted the true and recovered initial functions, the solution at the final time, the corresponding errors for σ = 1/2 and the smoothing of the measurements φ j (j = 1, 2) using polynomial curve fitting of the fifth degree.The accuracy improved significantly for both the initial function and solution, but the accuracy of the solution for σ = 1 was still higher.

Conclusions
In this paper, we studied a retrospective inverse problem, consisting of reconstruction of the initial data from the final time observation for a parabolic problem defined on disjoint intervals.The problem was ill-posed, and for its approximate solution, we suggested and validated a second-order accurate difference scheme.We proved the positive definiteness of the basic semidiscrete space difference operator in a Sobolev weighted norm.For solving the resulting linear difference system of equations, we developed an iterative conjugate gradient algorithm.Numerical experiments are provided, demonstrating the efficiency of the proposed approach.Also, the numerical results validate the theoretical statements and show that for σ = 1, the proposed numerical algorithm achieved better precision in recovering the initial condition and solution, in contrast to the smaller weights (σ = 1/2).
A natural extension of the present approach would be concerned with two-or multidimensional parabolic equations with linear or nonlinear interface conditions (see, for example, [4]).
Let us note that in [37], the quasi-solution to problem for identifying the initial data from the final time observations was studied.We plan to combine this method with the one in the present paper for studying the retrospective inverse problem for the time-fractional diffusion system in [27,28].
Finally, it is important to note that the authors of [44] solved the inverse problem of recovering the initial condition of a degenerate parabolic equation on the basis of final time observations by implementing the Landberg iteration method.In future work, we intend to develop this idea to a retrospective inverse problem for a degenerate atmospheric model [45,46].

Figure 1 .
Figure 1.Exact solution u 1 (x, T), u 2 (x, T) (solid line) and numerical solution u T 1,i 1 , u T 2,i 2 (line with circles) of the direct problem at the final time (left) and the corresponding error (right), where σ = 1/2 and τ = h (Example 1).

Figure 2 .
Figure 2.Exact solution u 1 (x, T), u 2 (x, T) (solid line) and numerical solution to the direct problem u T 1,i 1 , u T 2,i 2 (line with circles) at the final time (left) and the corresponding error (right), where σ = 1 and τ = h (Example 1).
2 for a sufficiently small τ.

Table 1 .
Errors and spatial convergence rate of the solution to the direct problem, where σ = 1/2 and τ = h (Example 1).

Table 2 .
Errors and spatial convergence rate of the solution to the direct problem, where (σ = 1 and τ = h 2 ) (Example 1).