Numerical Recovering of Space-Dependent Sources in Hyperbolic Transmission Problems

Abstract

A body may have a structural, thermal, electromagnetic or optical role. In wave propagation, many models are described for transmission problems, whose solutions are defined in two or more domains. In this paper, we consider an inverse source hyperbolic problem on disconnected intervals, using solution point constraints. Applying a transform method, we reduce the inverse problems to direct ones, which are studied for well-posedness in special weighted Sobolev spaces. This means that the inverse problem is said to be well posed in the sense of Tikhonov (or conditionally well posed). The main aim of this study is to develop a finite difference method for solution of the transformed hyperbolic problems with a non-local differential operator and initial conditions. Numerical test examples are also analyzed.

1. Introduction

Many phenomena in science and engineering are described by hyperbolic partial differential equations (HPDEs). Typical properties of HPDEs are finite propagation and Huygen’s principle. In the context of mathematical modeling, hyperbolic systems describe wave propagation, processes in fluid dynamics and heat conduction.

Layers with material properties that differ substantially from those of the ambient medium occur in a wide range of applications [1,2,3]. The layer may assume structural, thermal, electromagnetic or optical roles. Processes in regions with layers can be described by boundary value problems whose solutions are defined in two (or more) regions. In some instances these domains are disconnected. Such a situation arises, for example, when the solution in the medial domain is known or can be obtained by solving a simpler equation. The influence of the intermediate domain can be accounted for by non-local jump conditions [3,4]. A class of interface finite hyperbolic problems is called transmission problems.

The isotropic elastic wave equation in a bounded domain with a coefficient jump at a nested set of interfaces, satisfying natural transmission conditions, is investigated in [5]. The microlocal behavior of such solutions, like reflection, transmission and mode conversion of sand P waves, evanescent models, and Rayleigh and Stoneley waves, are analyzed. These solutions are of particular interest to scientists and practitioners (see, e.g., [2,3]).

The paper [6] studies the asymptotic behavior of 1D bodies that are composed of two different types of materials—one thermoelastic and the other without thermal effects. A mathematical model is used to represent the hyperbolic transmission problem.

The transmission problem for semilinear parabolic–hyperbolic equations in a multidimensional structure is studied in [7]. Rote’s method is used to prove the existence and uniqueness of the solution.

In [8], the authors investigate an initial boundary value problem for one-dimensional hyperbolic equations on two disconnected intervals. They construct and analyze a finite difference scheme for solving this problem. The corresponding 2D problem is studied in [9]. The existence and uniqueness of the solution is established and a priori estimates for the weak solutions in appropriate Sobolev-like spaces are obtained. Additionally, finite difference schemes are proposed and the convergence is proved.

In many mathematical problems, modeling heat conduction, thermoelasticity, or processes in biology, chemistry, environmental pollution, finance, industry, etc., apart from the solution, one or more of the following are unknown: coefficients, source terms, boundary conditions, boundaries (Cauchy problem) or initial conditions (retrospective inverse problem), since these quantities cannot be directly measured. In addition, some overspecified data are given. Such problems are called inverse problems. Typically they are ill-posed in the sense of Hadamard, i.e., the solution does not exist or it is not unique. In addition, the solution may be unstable [10,11,12,13,14]. The ill-posedness is usually overcome by using regularization methods [10,12,14,15,16,17]. Another often used approach is to transform the inverse problem to a well-posed direct problem [14,18,19,20].

Recently, many investigations have been devoted to inverse problems for PDEs, in particular, inverse problems for hyperbolic-type equations. Inverse problems for HPDEs in the connected domain are well studied in the literature (see, e.g., [13,15,17,21] and references therein).

The uniqueness and stability analysis for an inverse problem for recovering the source function in evolution problems from final time measurements are studied in [13]. One-dimensional and multi-dimensional parabolic, hyperbolic, Euler–Bernoulli beam and Kirchhoff plate equations are investigated. In [22], the author studies theoretically and numerically an inverse source problem for the Euler–Bernoulli beam and Kirchhoff–Love plate. A numerical algorithm, based on the Landweber–Fridman iterative regularization method is developed. A boundary control method is utilized in [21] to solve the inverse coefficient problem for wave equations. The identification of the spatial part of the source term of a degenerate wave equation based on the final observation data is studied in [17]. The inverse problem is formulated as a nonlinear optimization problem. The Tikhonov regularization technique is used to handle the noisy measurements.

An inverse problem for recovering two space-dependent coefficients in an acoustic equation of the hyperbolic type, on the basis of interior measurements of solutions, is considered in [23]. The uniqueness of the solution of the inverse problem is established using Lipschitz stability estimates, on the basis of a Carleman type estimate.

The authors of [15] introduce a two-dimensional wave equation on a double connected domain and study the inverse problem of recovering the interior boundary curve from overspecified data on the exterior boundary. The inverse problem is formulated as a system of boundary integral equations. An iterative numerical approach is constructed, linearizing the equation by the Frechet derivative. The generated ill-conditioned linear system is solved, applying Tikhonov regularization.

Papers [18,19] deal with inverse source problems in 1D and 2D thermoelasticity systems, considering a coupled system of a parabolic (heat) equation and a vectorial hyperbolic equation for the displacement. In [18], a 1D problem is considered for determination of the time-dependent heat source from temperature measurements inside the body. Eliminating the source function by observation, the inverse problem is reduced to a direct problem.

The authors prove the existence and uniqueness of the solution and construct a numerical algorithm for solving the problem. Recovery of a space-dependent vector source under final time measurements in a coupled hyperbolic–parabolic system is studied in [19]. The uniqueness of the solution is established and a stable iterative numerical algorithm, based on succession of well-posed direct problems, is constructed for solving the inverse problem.

In [20], a numerical algorithm is proposed based on a decomposition technique for determination of the right-hand side in a time-fractional parabolic equation, defined on disjoint domains.

A retrospective inverse problem for parabolic equations on disjoint domains is considered in [16]. The inverse problem is discretized by weighted time discretization. Then, an iterative conjugate gradient method is used to solve the obtained ill-posed systems of the difference equations.

To the best of our knowledge, in contrast to the hyperbolic standard transmission problems, results for inverse problems for hyperbolic transmission problems on disconnected domains are missing in the literature. Precisely because of this, the main goal of the present work is to construct efficient numerical methods for identifying a space-dependent source term in a hyperbolic problem on disconnected intervals.

The motivation is explained by the capacity of the hyperbolic equations to describe a large number of problems in elasticity, thermoelasticity and wave propagation (see, e.g., [2,3,4,5,6,12,18,22]). Also, it comes from the theoretical and numerical challenges of the present specific problem (see, e.g., [5,8,9,10,11,12,15,17,18,21,23]).

The remaining part of the paper is organized as follows: In the next section, the direct (forward) problem is formulated and well-posedness is discussed. Section 3 is devoted to the inverse problem for identification of the source in each domain based upon fixed time measurements. We propose a method for reducing the inverse problem to a direct non-local-in-time problem and show that it is well posed. In Section 4, efficient numerical approaches are developed for the solution of the inverse problem. Numerical simulations are discussed in Section 5 and then the paper is completed with some conclusions.

2. The Forward Problem

In this section, we introduce a direct problem–initial boundary value problem for HPDE on disconnected intervals. This is a nonstandard interface problem with jump conditions of Robin’s type. Well-posedness is discussed as well.

2.1. Formulation of the Initial Boundary Value Problem

We consider the following model initial boundary value problem for unknown functions $u_1=u_1(x,t)$, $u_2=u_2(x,t)$, satisfying the system of hyperbolic equations [8]

(1)$\begin{array}{ccc}& & \frac{{\partial }^2{u}_1}{\partial t^2}-\frac{\partial }{\partial x}\left(p_1\left(x\right)\frac{\partial u_1}{\partial x}\right)=f_1\left(x\right)g_1\left(t\right),(x,t)\in Q_{1T},\hfill \end{array}$

(2)$\begin{array}{ccc}& & \frac{{\partial }^2{u}_2}{\partial t^2}-\frac{\partial }{\partial x}\left(p_2\left(x\right)\frac{\partial u_2}{\partial x}\right)=f_2\left(x\right)g_2\left(t\right),(x,t)\in Q_{2T},\hfill \end{array}$

where $Q_{jT}={\Omega }_j\times (0, T)$, ${\Omega }_j=(a_j,b_j)$, $j=1, 2$. The Robin–Dirichlet type internal boundary conditions

(3)$\begin{array}{ccc}& & p_1\left(b_1\right)\frac{\partial u_1}{\partial x}(b_1,t)+{\alpha }_1{u}_1(b_1,t)={\beta }_1{u}_2(a_2,t)+{\gamma }_1\left(t\right),t\in (0, T],\hfill \end{array}$

(4)$\begin{array}{ccc}& & -p_2\left(a_2\right)\frac{\partial u_2}{\partial x}(a_2,t)+{\alpha }_2{u}_2(a_2,t)={\beta }_2{u}_1(b_1,t)+{\gamma }_2\left(t\right),t\in (0, T],\hfill \end{array}$

and external Dirichlet boundary conditions

(5)$u_1(a_1,t)={\mu }_1\left(t\right),u_2(b_2,t)={\mu }_2\left(t\right),t\in (0, T],$

are imposed and initial conditions are given by

(6)$\begin{array}{ccc}& & u_1(x,0)={\phi }_1\left(x\right),\frac{\partial u_1}{\partial t}(x,0)={\psi }_1\left(x\right),x\in {\Omega }_1,\hfill \end{array}$

(7)$\begin{array}{ccc}& & u_2(x,0)={\phi }_2\left(x\right),\frac{\partial u_2}{\partial t}(x,0)={\psi }_2\left(x\right),x\in {\Omega }_2.\hfill \end{array}$

The sense of this simple model, described by the systems of Equations (1)–(7) is the vibration of two standard one-dimensional linear solids rod of length $b_j-a_j$, $j=1,2$, whose displacement at point x and t is modeled by functions $u_j$, $j=1,2$, which are solutions of (1) and (2). The source functions in the Equations (1) and (2) represent external force. At the left end of the first rod and right end of the second rod, a source ${\mu }_j$, $j=1,2$ is applied. Next, there is a displacement exchange between the rods, described by the interface conditions (3) and (4). In general, the problems of type (1)–(7) are very difficult to solve analytically and so a numerical approach is required.

Further, we assume

(8)$\begin{array}{cc}& {\displaystyle p_j\left(x\right)\in L_{\infty }\left({\Omega }_j\right),p_j\left(x\right)\ge p_{0j}>0,x\in {\Omega }_j,j=1,2}\end{array}$

(9)$\begin{array}{cc}\hfill & {\displaystyle {\alpha }_j>0,{\beta }_j>0,j=1,2,{\beta }_1{\beta }_2\le {\alpha }_1{\alpha }_2.}\end{array}$

The problem of the hyperbolic system (1) and (2), with a known coefficient and right-hand side, obeying the internal boundary conditions (3) and (4), external boundary conditions (5), and initial conditions (6) and (7), we call a direct or forward problem (see, e.g., [10,11,12,14]).

2.2. Well-Posedness of the Direct Problem

In this subsection, following the results of paper [8], we discuss the well-posedness, i.e., existence, uniqueness and continuous dependence on input data, of the problem (1)–(7).

For simplicity of exposition, we assume homogeneous external and internal boundary conditions, i.e., ${\mu }_j\left(t\right)=0$, and ${\gamma }_j\left(t\right)=0,t\in (0, T]$ for $j=1,2$. Otherwise, introducing variables

$\begin{array}{ccc}& & w_1(x,t)=u_1^1(x,t)-A_1\left(t\right)(x-a_1)(x-b_1),A_1\left(t\right)=\frac{{\gamma }_1\left(t\right)}{p_1\left(b_1\right)(b_1-a_1)},\hfill \\ & & w_2(x,t)=u_2^1(x,t)-A_2\left(t\right)(x-a_2)(x-b_2),A_2\left(t\right)=\frac{{\gamma }_2\left(t\right)}{p_2\left(b_2\right)(b_2-a_2),}\hfill \end{array}$

where

$u_1^1(x,t)=u_1(x,t)-{\phi }_1\left(x\right)-{\mu }_1\left(t\right),u_2^1(x,t)=u_2(x,t)-{\phi }_2\left(x\right)-{\mu }_2\left(t\right),$

yields to this homogeneous case.

Our results need the following construction: We introduce the product space endowed with the inner product and associated norm:

$L=L_2\left({\Omega }_1\right)\times L_2\left({\Omega }_2\right)=\{v=(v_1,v_2)|v_j\in L_2\left({\Omega }_i\right)\},{(u_j,v_j)}_{L_2\left({\Omega }_i\right)}={\int }_{{\Omega }_i}{u}_j{v}_{j}dx,j=1,2.$

Then, we define:

${(u,v)}_L={\beta }_2{(u_1,v_1)}_{L_2\left({\Omega }_1\right)}+{\beta }_1{(u_2,v_2)}_{L_2\left({\Omega }_2\right)},{\parallel v\parallel }_L=(v,v).$

We also define the Sobolev spaces:

$H^k=\left\{(v_1,v_2)\right|v_j\in H^k\left({\Omega }_j\right)\},k=1,2,\cdots ,$

where:

$\begin{array}{ccc}& & {(u,v)}_{H^k}={\beta }_2{(u_1,v_1)}_{H^k\left({\Omega }_1\right)}+{\beta }_1{(u_2,v_2)}_{H^k\left({\Omega }_2\right)},{\parallel v\parallel }_{H^k}={(v,v)}^{1/2},\hfill \\ & & {(u_j,v_j)}_{H^k\left({\Omega }_j\right)}=\sum _{j=0}^k{\left(\frac{d^i{u}_j}{d{x}^i},\frac{d^i{v}_j}{d{x}^i}\right)}_{L_2\left({\Omega }_j\right)},i=1,2,k=1,2,\cdots ,\hfill \\ & & H_0^1=\{v=(v_1,v_2)\in H^1|v_1\left(a_1\right)=0,v_2\left(b_2\right)=0\}.\hfill \end{array}$

Now, we can construct the bilinear form:

$\begin{array}{ccc}\hfill A(u,v)& =& {\beta }_2{\int }_{{\Omega }_1}{p}_1\frac{d{u}_1}{dx}\frac{d{v}_1}{dx}+{\beta }_1{\int }_{{\Omega }_2}{p}_2\frac{d{u}_2}{dx}\frac{d{v}_2}{dx}\hfill \\ & & +{\beta }_1{\alpha }_1{u}_1\left(b_1\right)v_1\left(b_1\right)+{\beta }_2{\alpha }_2{u}_2\left(a_2\right)v_2\left(a_2\right)-{\beta }_1{\beta }_2[u_1\left(b_1\right)v_2\left(a_2\right)+u_2\left(a_2\right)v_1\left(b_1\right)].\hfill \end{array}$

The following assertion holds.

Lemma 1 

([8]). Let the conditions (8) and (9) be fulfilled. Then, the bilinear form A is symmetric bounded on $H^1\times H^1$ and coercive on $H_0^1$, i.e., there exists a constant $c>0$, such that

$A(v,v)\ge {c\parallel u\parallel }_{H^1}^2,\forall v\in H_0^1.$

Hence, the basic result follows:

Theorem 1 

([8]). Suppose that the conditions (8) and (9) hold and

$\begin{array}{cc}& {\displaystyle u_0\left(x\right)=({\phi }_1\left(x\right),{\phi }_2\left(x\right))\in H_0^1,u_1\left(x\right)\equiv ({\psi }_1\left(x\right),{\psi }_2\left(x\right))\in L,}\\ & {\displaystyle \tilde{f}=({\tilde{f}}_1,{\tilde{f}}_2)\in L_2((0, T),L),{\tilde{f}}_i=f_i{g}_i,i=1,2.}\end{array}$

Then, the hyperbolic initial boundary value problem (1)–(7) has a unique weak solution $u\in L_2((0,T),H_0^1)$ and it depends continuously on $\tilde{f}$, φ and ψ.

3. Inverse Problem

We formulate the inverse problem for identifying the space-dependent part of the source term in the hyperbolic interface problem (1)–(7). A solution method based on the reduction of the inverse problem to a direct well-posed problem is discussed.

3.1. Formulation of the Inverse Problem

In this work, we consider an inverse problem for identifying the space-dependent parts $f_j\left(x\right)$ of the source terms in (1)–(2) and the solutions $u_j$, $j=1,2$ of (1)–(7), if additional observations of the solution at time $0<t^{*}\le T$ are given, namely,

(10)$u_j(x,t^{*})=m_j\left(x\right),x\in {\Omega }_j,j=1,2.$

In general, the problems of type (1)–(9), as well as the corresponding inverse problem (1)–(10), are very difficult to solve analytically and so a numerical approach is required.

Inverse source problems in wave propagations in a single domain are well studied (see, e.g., [13,15,17,21,22]). In contrast, there are no results in the literature for inverse hyperbolic problems on separate intervals. The sense of the simple model, described by the system of Equations (1)–(7) and data (10), is the vibration of two interacting one-dimensional rods of length $b_j-a_j$, $j=1,2$ under a given vertical displacement at time $t=t^{*}$. In some cases, the external forces cannot be measured. However, having observations of the displacement in a fixed moment $t^{*}$, our aim is to reconstruct approximately the space-dependent parts of the source functions.

3.2. Solving Inverse Problem

In this section, we discuss the solution method for the inverse problem. It extends the idea proposed in [24] for the Dirichlet problem for a parabolic equation on a single domain in order to convert the inverse source hyperbolic problem on disjoint domains to a direct one.

Let

(11)$v_j=\frac{\partial u_j}{\partial t},j=1,2.$

Differentiating the Equations (1)–(4) with respect to t yields

(12)$\begin{array}{ccc}& & \frac{{\partial }^2{v}_1}{\partial t^2}-\frac{\partial }{\partial x}\left(p_1\left(x\right)\frac{\partial v_1}{\partial x}\right)=f_1\left(x\right)\frac{d{g}_1\left(t\right)}{dlt},(x,t)\in Q_{1T},\hfill \end{array}$

(13)$\begin{array}{ccc}& & \frac{{\partial }^2{v}_2}{\partial t^2}-\frac{\partial }{\partial x}\left(p_2\left(x\right)\frac{\partial v_2}{\partial x}\right)=f_2\left(x\right)\frac{d{g}_2\left(t\right)}{dt},(x,t)\in Q_{2T},\hfill \end{array}$

(14)$\begin{array}{ccc}& & p_1\left(b_1\right)\frac{\partial v_1}{\partial x}(b_1,t)+{\alpha }_1{v}_1(b_1,t)={\beta }_1{v}_2(a_2,t)+\frac{d{\gamma }_1\left(t\right)}{dt},t\in (0,T],\hfill \end{array}$

(15)$\begin{array}{ccc}& & -p_2\left(a_2\right)\frac{\partial v_2}{\partial x}(a_2,t)+{\alpha }_2{v}_2(a_2,t)={\beta }_2{v}_1(b_1,t)+\frac{d{\gamma }_2\left(t\right)}{dt},t\in (0,T],\hfill \end{array}$

(16)$\begin{array}{ccc}& & v_1(a_1,t)=\frac{d{\mu }_1\left(t\right)}{dt},v_2(b_2,t)=\frac{d{\mu }_2\left(t\right)}{dt},t\in (0,T],\hfill \end{array}$

From (6) and (7), we obtain

(17)$v_j(x,0)={\psi }_j\left(x\right),x\in {\Omega }_j,j=1,2.$

Next, from (1)–(4) for $t=0$ and (6) and (7), we derive

(18)$\begin{array}{ccc}& & \frac{\partial v_1}{\partial t}(x,0)-\frac{d}{dx}\left(p_1\left(x\right)\frac{d{\phi }_1\left(x\right)}{dx}\right)=f_1\left(x\right)g_1\left(0\right),x\in {\Omega }_1,\hfill \end{array}$

(19)$\begin{array}{ccc}& & \frac{\partial v_2}{\partial t}(x,0)-\frac{d}{dx}\left(p_2\left(x\right)\frac{d{\phi }_2\left(x\right)}{dx}\right)=f_2\left(x\right)g_2\left(0\right),x\in {\Omega }_2,\hfill \end{array}$

(20)$\begin{array}{ccc}& & p_1\left(b_1\right)\frac{d{\phi }_1}{dx}\left(b_1\right)+{\alpha }_1{\phi }_1\left(b_1\right)={\beta }_1{\phi }_2\left(a_2\right)+{\gamma }_1\left(0\right),\hfill \end{array}$

(21)$\begin{array}{ccc}& & -p_2\left(a_2\right)\frac{d{\phi }_2}{dx}\left(a_2\right)+{\alpha }_2{\phi }_2\left(a_2\right)={\beta }_2{\phi }_1\left(b_1\right)+{\gamma }_2\left(0\right),\hfill \end{array}$

Finally, we use (1) and (2) at $t=t^{*}$ and the overspecified data (10) to determine

(22)$\begin{array}{ccc}& & {\displaystyle f_1\left(x\right)=\frac{\frac{\partial v_1}{\partial t}(x,t^{*})-\frac{d}{dx}\left(p_1\left(x\right)\frac{d{m}_1\left(x\right)}{dx}\right)}{g_1\left(t^{*}\right)},x\in {\Omega }_1,}\hfill \end{array}$

(23)$\begin{array}{ccc}& & {\displaystyle f_2\left(x\right)=\frac{\frac{\partial v_2}{\partial t}(x,t^{*})-\frac{d}{dx}\left(p_2\left(x\right)\frac{d{m}_2\left(x\right)}{dx}\right)}{g_2\left(t^{*}\right)},x\in {\Omega }_2.}\hfill \end{array}$

Therefore, to recover the source functions $f_1\left(x\right)$, $f_2\left(x\right)$, we need to solve the forward non-local problem (12)–(23). We have two hyperbolic equations with time-loaded terms (see, e.g., [14]) ${\displaystyle \frac{\partial v_j}{\partial t}(x,t^{*})}$, $j=1,2$, subject with boundary conditions (14)–(16) and initial conditions (17)–(19). Let us note that (18)–(19) are non-local boundary conditions and (20)–(21) are only compatibility conditions concerning ${\phi }_j$ and ${\gamma }_j$, $j=1,2$.

3.3. Well-Posedness of the Problem (12)–(19)

We introduce the notation $\dot{v}=\frac{dv}{dt}$ and substitute $f_j$, $j=1,2$ from (22) and (23) in Equations (12)–(19) to obtain the equivalent problem for $v_j$, $j=1,2$

(24)$\begin{array}{ccc}& & \frac{{\partial }^2{v}_1}{\partial t^2}-\frac{\partial }{\partial x}\left(p_1\left(x\right)\frac{\partial v_1}{\partial x}\right)-\frac{{\dot{g}}_1\left(t\right)}{g_1\left(t^{*}\right)}\frac{\partial v_1}{\partial t}(x,t^{*})=F_1(x,t),\hfill \end{array}$

(25)$\begin{array}{ccc}& & \frac{{\partial }^2{v}_2}{\partial t^2}-\frac{\partial }{\partial x}\left(p_2\left(x\right)\frac{\partial v_2}{\partial x}\right)-\frac{{\dot{g}}_2\left(t\right)}{g_2\left(t^{*}\right)}\frac{\partial v_2}{\partial t}(x,t^{*})=F_2(x,t),\hfill \end{array}$

where ${\displaystyle F_j(x,t)=-\frac{{\dot{g}}_j\left(t\right)}{g_j\left(t^{*}\right)}\frac{d}{dx}\left(p_j\left(x\right)\frac{d{m}_j\left(x\right)}{dx}\right)}$, with initial conditions (17) and

(26)$\begin{array}{ccc}& & \frac{\partial v_1}{\partial t}(x,0)-\frac{g_1\left(0\right)}{g_1\left(t^{*}\right)}\frac{\partial v_1}{\partial t}(x,t^{*})=\frac{d}{dx}\left(p_1\left(x\right)\frac{d{\phi }_1\left(x\right)}{dx}\right)+F_1(x,0),\hfill \end{array}$

(27)$\begin{array}{ccc}& & \frac{\partial v_2}{\partial t}(x,0)-\frac{g_2\left(0\right)}{g_2\left(t^{*}\right)}\frac{\partial v_2}{\partial t}(x,t^{*})=\frac{d}{dx}\left(p_2\left(x\right)\frac{d{\phi }_2\left(x\right)}{dx}\right)+F_2(x,0)\hfill \end{array}$

and boundary conditions (14)–(16).

Theorem 2. 

Let the assumptions (8) and (9) hold and the internal and external boundary conditions are homogeneous. Suppose that $\{{\psi }_1,{\psi }_2\}\in L$. Then, the initial boundary value problem (24)–(27) has a unique solution $v=(v_1,v_2)\in L_2\left((0,T),H_0^1\right)$ and depends continuously on the input data

Proof. 

The proof is performed in two steps.

Step 1. Let

(28)$\frac{\partial v_j}{\partial t}(x,t)=w_j(x,t)⟹v_j(x,t)={\int }_0^t{w}_j(x,s)ds+v_j(x,0)=I_{j}w+{\psi }_j\left(x\right),j=1,2.$

Thus, the problem (24)–(27) with (17)–(19) is equivalent to the following one: find the functions $w_j(x,t)$, $j=1,2$, satisfying

(29)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{\partial w_j}{\partial t}-\frac{\partial }{\partial x}\left(p_j\left(x\right)\frac{\partial }{\partial x}{I}_{j}w\right)-\frac{{\dot{g}}_j\left(t\right)}{g_j\left(t^{*}\right)}{w}_j(x,t^{*})\hfill \\ & =F_j(x,t)+\frac{\partial }{\partial x}\left(p_j\left(x\right)\frac{d{\psi }_j\left(x\right)}{dx}\right)-\frac{{\dot{g}}_j\left(t\right)}{g_j\left(t^{*}\right)}{\psi }_j\left(x\right),j=1,2,\hfill \end{array}\end{array}$

with initial conditions

(30)$\frac{\partial w_j}{\partial t}(x,0)-\frac{g_j\left(0\right)}{g_j\left(t^{*}\right)}{w}_j(x,t^{*})=\frac{d}{dx}\left(p_j\left(x\right)\frac{d{\phi }_j\left(x\right)}{dx}\right)+F_j(x,0),j=1,2,$

internal boundary conditions

(31)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & p_1\left(b_1\right)\frac{\partial w_1}{\partial x}(b_1,t)+{\alpha }_1{w}_1(b_1,t)={\beta }_1{w}_2(a_2,t)+\frac{d^2{\gamma }_1\left(t\right)}{d^{2}t},t\in (0,T],\hfill \\ & -p_2\left(a_2\right)\frac{\partial w_2}{\partial x}(a_2,t)+{\alpha }_2{w}_2(a_2,t)={\beta }_2{w}_1(b_1,t)+\frac{d{\gamma }_2^2\left(t\right)}{d{t}^2},t\in (0,T],\hfill \end{array}\end{array}$

and external boundary conditions

(32)$w_1(a_1,t)=\frac{d^2{\mu }_1\left(t\right)}{d{t}^2},w_2(b_2,t)=\frac{d^2{\mu }_2\left(t\right)}{d{t}^2},t\in (0,T].$

Step 2. For any positive integer $M^{*}$ and function $w(x,t)$, denote

$w^{M^{*},m}\left(x\right)=w(x,m\tau ),m=0,1,\cdots ,M^{*},$

where $\tau =t^{*}/M^{*}$.

We consider the approximating equations of (30)

(33)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{w_j^{M^{*},m}\left(x\right)-w_j^{M^{*},m-1}\left(x\right)}{\tau }-\frac{\partial }{\partial x}\left(p_1\left(x\right)\tau \sum _{s=1}^m{w}_j^{M^{*},m}\right)-\frac{{\dot{g}}_j\left(t^m\right)}{g_j\left(t^{*}\right)}{w}_j(x,t^{*})\hfill \\ & =-\frac{{\dot{g}}_j\left(t^m\right)}{g_j\left(t^{*}\right)}\frac{d}{dx}\left(p_r\left(x\right)\frac{d{m}_j\left(x\right)}{dx}\right),w_j(x,t^{*})=w_j^{M^{*},M^{*}},j=1,2,\hfill \end{array}\end{array}$

and (30)–(32) in a similar way.

For $m=1$, we take

$\frac{d{v}_j^{M,0}\left(x\right)}{dx}=\frac{d{\psi }_j\left(x\right)}{dx},j=1,2.$

Thus, from (33), we obtain

$w_j^{M,1}\left(x\right)=\frac{{\dot{g}}_j\left(t^1\right)}{g_j\left(t^{*}\right)}\frac{\partial v_j}{\partial t}(x,t^{*})+\tau \frac{\partial }{\partial x}\left(p_j\left(x\right)\frac{v_j^{M,0}}{\partial x}\right)+\tau F(x,t^m).$

We define an approximate (Rothe) solution for $t\in (0,\tau )$ by interpolation

$w_j^{M,1}(x,t)=\frac{(\tau -t)}{\tau }{w}^{M,0}\left(x\right)+\frac{t}{\tau }{w}_j^{M,1}\left(x\right),j=1,2.$

Then, with the corresponding approximation of $v_i(x,t)$, we have

$v_j^1(x,t)\approx \left(-t+\frac{t^2}{2\tau }\right)w^{M,0}\left(x\right)+\frac{t^2}{2\tau }{w}_j^{M,1}\left(x\right),j=1,2,$

for $0\le t<\tau .$

Further, using the trace inequality (see, e.g., [25])

$z^2\left(d\right)\le \frac{c}{{(d-c)}^2}{\parallel z\parallel }_{H^1(c,d)}^2,c=0,d=t^{*},z=w_j(x,t^{*}),j=1,2,$

we write a weak formulation similar to that of the direct problem, Section 2.2, and then obtain a priori estimates of the solution $\{w_1,w_2\}$.Then, on this base, following the Rothe’s methodology, see e.g., [7,18,19,26], we prove it’s unique convergence as $M^{*}\to \infty$ to the unique solution of (29)–(32) on the interval $[0,t^{*}]$.    □

4. Numerical Method

In this section, we construct numerical methods for solving the problem (12)–(13).

We develop iterative finite difference methods based on concepts from papers [27,28,29,30,31], which address the solution of linear parabolic problems with non-local terms only in the initial condition, and [32], where a non-local term occurs both in the the differential operator and the initial condition, but again in the linear parabolic problem on a single connected domain.

Let us define a uniform mesh in time and space in each domain ${\Omega }_j$, $j=1,2$,

$\begin{array}{ccc}& & {\overline{\omega }}_{\tau }=\{t_n=n\tau ,n=0,1,\cdots ,M,M\tau =T\},\hfill \\ & & {\overline{\omega }}_{hj}=\{x_{j,i_j}=i_j{h}_j,i_j=1,\cdots ,N_j,x_{j,0}=a_j,x_{j,N_j}=b_j\},j=1,2\hfill \end{array}$

and denote by $V_{j,i_j}^n=v_j(x_{j,i_j},t_n)$ the approximate value of the function $v_j$ at grid node $(x_{j,i_j},t_n)\in Q_{jT}$, $j=1,2$. We consider the case when $t^{*}$ is a grid node and $t^{*}=n^{*}\tau$. We introduce also the notations [33]

$\begin{array}{cc}& {\displaystyle p_{j,i_j\pm 1/2}=p_j\left(x_{j,i_j\pm 1/2}\right),V_{j,x_{i_j}}=\frac{V_{j,i_j+1}-V_{j,i_j}}{h_j},V_{j,{\overline{x}}_{i_j}}=\frac{V_{j,i_j}-V_{j,i_j-1}}{h_j},}\\ & {\displaystyle V_{j,t_{i_j}}^n=\frac{V_{j,i_j}^n-V_{j,i_j}^{n-1}}{\tau },V_{j,\overline{t}{t}_{i_j}}^{n+1}=\frac{V_{j,i_j}^{n+1}-2V_{j,i_j}^n+V_{j,i_j}^{n-1}}{{\tau }^2},}\\ & {\displaystyle {\Lambda }_{i_j}\left(V_j^n\right)=\frac{1}{h_j}\left(p_{j,i_j+1/2}{V}_{j,x_{i_j}}^n-p_{j,i_j-1/2}{V}_{j,{\overline{x}}_{i_j}}^n\right),P_j\left(x\right)=\frac{d}{dx}\left(p_j\left(x\right)\frac{d{\phi }_1\left(x\right)}{dx}\right)}\\ & {\displaystyle G_j=\frac{d{g}_j}{dt},{\Phi }_j=\frac{d{\phi }_j}{dx},{\Psi }_j=\frac{d{\psi }_j}{dx},{\Gamma }_j=\frac{d{\gamma }_j}{dt},{\mathcal{M}}_j=\frac{d{\mu }_j}{dt}.}\end{array}$

4.1. Simple Discretization

We apply the finite difference method, implicit time-stepping for the differential equations at the inner domain, and implicit-explicit time-stepping for the interface conditions. Using second-order central finite difference approximations for the temporal and spatial derivatives in (12) and (13), we obtain

(34)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & V_{1,\overline{t}{t}_{i_1}}^{n+1}-{\Lambda }_{i_1}\left(V_1^{n+1}\right)=f_{1,i_1}{G}_1^{n+1},n=1,2\cdots ,M-1,i_1=1,2,\cdots ,N_1-1,\hfill \\ & V_{2,\overline{t}{t}_{i_2}}^{n+1}-{\Lambda }_{i_2}\left(V_2^{n+1}\right)=f_{2,i_2}{G}_2^{n+1},n=1,2\cdots ,M-1,i_2=1,2,\cdots ,N_2-1.\hfill \end{array}\end{array}$

For $n=0$, the Equation (34) involves the solutions $V_j$ at the time layer $n=-1$, namely,

$\begin{array}{c}\hfill \begin{array}{cc}& \frac{V_{j,i_1}^1}{{\tau }^2}-{\Lambda }_{i_1}\left(V_1^1\right)=\frac{2{\psi }_{j,i_1}-V_{j,i_1}^{-1}}{{\tau }^2}+f_{1,i_1}{G}_1^1,i_1=1,2,\cdots ,N_1-1,\hfill \\ & \frac{V_{j,i_2}^1}{{\tau }^2}-{\Lambda }_{i_2}\left(V_2^1\right)=\frac{2{\psi }_{j,i_2}-V_{j,i_2}^{-1}}{{\tau }^2}+f_{2,i_2}{G}_2^1,i_2=1,2,\cdots ,N_2-1.\hfill \end{array}\end{array}$

Thus, we apply a second-order central finite difference approximation for the initial condition (18)–(19) to obtain

(35)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{2V_{j,i_1}^1}{{\tau }^2}-{\Lambda }_{i_1}\left(V_1^1\right)=\frac{2{\psi }_{j,i_1}}{{\tau }^2}+\frac{2P_{1,i_1}}{\tau }+\left(\frac{2g_1^0}{\tau }+G_1^1\right)f_{1,i_1},i_1=1,2,\cdots ,N_1-1,\hfill \\ & \frac{2V_{j,i_2}^1}{{\tau }^2}-{\Lambda }_{i_2}\left(V_2^1\right)=\frac{2{\psi }_{j,i_2}}{{\tau }^2}+\frac{2P_{2,i_2}}{\tau }+\left(\frac{2g_2^0}{\tau }+G_2^1\right)f_{2,i_2},i_2=1,2,\cdots ,N_2-1.\hfill \end{array}\end{array}$

The boundary conditions (16), (20) and (21) are approximated using first-order finite differences [33]

(36)$\begin{array}{c}\hfill \begin{array}{c}{V}_{1,0}^{n+1}={\mathcal{M}}_1^{n+1},n=0,1,\cdots ,n^{*}-1,\hfill \\ p_{1,N_1}\frac{V_{1,N_1}^{n+1}-V_{1,N_1-1}^{n+1}}{h_1}+{\alpha }_1{V}_{1,N_1}^{n+1}={\beta }_1{V}_{2,0}^n+{\Gamma }_1^{n+1},\hfill \\ p_{2,0}\frac{V_{2,1}^{n+1}-V_{2,0}^{n+1}}{h_2}+{\alpha }_2{V}_{2,0}^{n+1}={\beta }_2{V}_{1,N_1}^n+{\Gamma }_2^{n+1},\hfill \\ V_{2,N_2}^{n+1}={\mathcal{M}}_2^{n+1},n=0,1,\cdots ,n^{*}-1.\hfill \end{array}\end{array}$

Since the interface conditions are discretized by first-order approximation, we cannot expect the spatial order of accuracy of the numerical scheme (34)–(36) to be more than one.

The right-hand side is determined from the approximation of (22) and (23)

(37)$f_{j,i_j}=\frac{V_{j,t_{i_j}}^{n^{*}}-{\Lambda }_{i_j}\left(m_j\right)}{g_j^{n^{*}}},i_j=1,2,\cdots ,N_j-1,j=1,2,$

where $V_j^{n^{*}}$ is the solution $V_j$ at time layer $t^{*}$. Note that to determine ${\Lambda }_{i_j}\left(m_j\right)$ at $j=1$, $i_1=1$ and $j=2$, $i_2=N_2-1$, we use external boundary conditions (5), i.e., $m_{1,0}={\mu }_1^{n^{*}}$ and $m_{2,N_2}={\mu }_2^{n^{*}}$. The values ${\Lambda }_{N_1-1}\left(m_1\right)$ and ${\Lambda }_1\left(m_2\right)$ are computed as follows: first, we represent the second term in the nominator of (22) and (23) in a nondivergent form, and then we replace the corresponding derivatives by one-sided finite differences, namely,

$\begin{array}{cc}\hfill & {\displaystyle {\Lambda }_{N_1-1}\left(m_1\right)={\tilde{P}}_1\left(x_{1,N_1-1}\right),{\Lambda }_1\left(m_2\right)={\tilde{P}}_2\left(x_{2,1}\right),}\\ & {\displaystyle {\tilde{P}}_j\left(x\right)=\frac{d{p}_j\left(x\right)}{dx}\frac{d{m}_j\left(x\right)}{dx}+p_j\left(x\right)\frac{d^2{m}_j\left(x\right)}{d{x}^2},}\end{array}$

where the first and second derivatives of $m_j$ are approximated by left and right second-order discretizations

(38)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{d{m}_1}{dx}\left(x_{1,N_1-1}\right)=m_{1,{\stackrel{<}{x}}_{N_1-1}}+O\left(h_1^2\right),\frac{d{m}_2}{dx}\left(x_{1,1}\right)=m_{2,{\stackrel{>}{x}}_1}+O\left(h_2^2\right),\hfill \\ & \frac{d^2{m}_1}{d^{2}x}\left(x_{1,N_1-1}\right)=m_{1,{\stackrel{<}{xx}}_{N_1-1}}+O\left(h_1^2\right),\frac{d^2{m}_2}{d{x}^2}\left(x_{1,1}\right)=m_{2,{\stackrel{>}{xx}}_1}+O\left(h_2^2\right),\hfill \end{array}\end{array}$

using the notations

$\begin{array}{ccc}& & m_{j,{\stackrel{<}{x}}_{i_j}}=\frac{3m_{j,i_j}-4m_{j,i_j-1}+m_{j,i_j-2}}{h_1},\hfill \\ & & m_{j,{\stackrel{>}{x}}_{i_j}}=\frac{-3m_{j,i_j}+4m_{j,i_j+1}-m_{j,i_j+2}}{h_2},\hfill \\ & & m_{j,{\stackrel{<}{xx}}_{i_j}}=\frac{2m_{j,i_j}-5m_{1,i_j-1}+4m_{1,i_j-2}-m_{1,i_j-3}}{h_1^2},\hfill \\ & & m_{j,{\stackrel{>}{xx}}_{i_j}}=\frac{2m_{j,i_j}-5m_{j,i_j+1}+4m_{j,i_j+2}-m_{j,i_j+3}}{h_2^2}.\hfill \end{array}$

The numerical recovering is performed by an iteration process for $k=0,1,\cdots$. Denote by $V_j^{\left(k\right),n}$ the solution $V_j$ at the k-th iteration, at time layer $t_n$. We start with an initial guess for $f_j$, $j=1,2$ and compute the problem (34)–(36) in the whole domain ${\overline{\omega }}_{\tau }\times {\overline{\omega }}_{h1}\times {\overline{\omega }}_{h2}$ to find $V_{j,i_j}^{\left(1\right),n}$, $i_j=0,1,\cdots ,N_j$, $n=1,2,\cdots ,M$, $j=1,2$. Then, we update the function $f_{j,i_j}$ by (37) and compute again the scheme (34)–(36) to find the solution at the next iteration. This process continues up to reaching the desired accuracy $ϵ$, i.e., $max\{\parallel V_1^{(k+1),n+1}-V_1^{\left(k\right),n+1}\parallel ,\parallel V_2^{(k+1),n+1}-V_2^{\left(k\right),n+1}\parallel \}\le ϵ$, where $\parallel v_j\parallel =\underset{0\le i_j\le N_j}{max}\left|v_{i_j}\right|$.

Once the functions $f_j$ are restored, we may compute the solution $U_j$, $j=1,2$.

If $t^{*}<T$, from (11), (6) and (7), we find $U_{j,i_j}^n$, $n=1,2,\cdots ,n^{*}$, $i_j=0,1,\cdots ,N_j$, $j=1,2$, as follows

(39)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & U_{j,i_j}^1=U_{j,i_j}^{-1}+2\tau V_{j,i_j}^1=U_{j,i_j}^1-2\tau ({\psi }_{j,i_j}-V_{j,i_j}^1),\hfill \\ & U_{j,i_j}^{n+1}=U_{j,i_j}^{n-1}+2\tau V_{j,i_j}^{n+1},n=1,2,\cdots ,n^{*}-1,\hfill \\ & U_{1,0}^{n+1}={\mu }_1^{n+1},U_{2,N_2}^{n+1}={\mu }_2^{n+1}.\hfill \end{array}\end{array}$

Further, to obtain the solution $U_{j,i_j}^{n+1}$, $n=n^{*},n^{*}+1,\cdots ,M-1$, $i_j=0,1,\cdots ,N_1$, $j=1,2$, we approximate the direct problem (1)–(7), namely,

(40)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & U_{1,0}^{n+1}={\mu }_1^{n+1},\hfill \\ & U_{1,\overline{t}{t}_{i_1}}^{n+1}-{\Lambda }_{i_1}\left(U_1^{n+1}\right)=f_{1,i_1}{g}_1^{n+1},i_1=1,2,\cdots ,N_1-1,\hfill \\ & p_{1,N_1}\frac{U_{1,N_1}^{n+1}-U_{1,N_1-1}^{n+1}}{h_1}+{\alpha }_1{U}_{1,N_1}^{n+1}={\beta }_1{U}_{2,0}^n+{\gamma }_1^{n+1},\hfill \\ & p_{2,0}\frac{U_{2,1}^{n+1}-U_{2,0}^{n+1}}{h_2}+{\alpha }_2{U}_{2,0}^{n+1}={\beta }_2{U}_{1,N_1}^n+{\gamma }_2^{n+1},\hfill \\ & U_{2,\overline{t}{t}_{i_2}}^{n+1}-{\Lambda }_{i_2}\left(U_2^{n+1}\right)=f_{2,i_2}{g}_2^{n+1},i_2=1,2,\cdots ,N_2-1,\hfill \\ & U_{2,N_2}^{n+1}={\mu }_2^{n+1}.\hfill \end{array}\end{array}$

Note that now in the term $U_{j,\overline{t}{t}_{i_j}}^{n+1}$, the values $U_{j,i_j}^n$ and $U_{j,i_j}^{n-1}$ are known from (39).

If $t^{*}=T$, then the solution $U_{j,i_j}^n$, $n=1,2,\cdots ,M$, $i_j=0,1,\cdots ,N_1$, $j=1,2$ is computed only from (39), $n^{*}=M$.

The above steps are described in Algorithm 1.

Algorithm 1: Recovering by $O\left(h\right)$ discretization

Require: $a_j$, $b_j$, T, $0<t^{*}\le T$, $m_j\left(x\right)$, ${\phi }_j\left(x\right)$, ${\psi }_j\left(x\right)$, ${\mu }_j\left(x\right)$, $p_j\left(x\right)$, $x\in {\Omega }_j$, ${\alpha }_j$, ${\beta }_j$,

${\gamma }_j\left(t\right)$, $g_j\left(t\right)$, $t\in [0,t^{*}]$, $ϵ$, $N_j$, M, $j=1,2$;

Ensure: $f_{j,i_j}$, $i_j=1,2,\cdots ,N_j-1$, $U_{j,i_j}^n$, $i_j=0,1,\cdots ,N_j$, $n=1,2,\cdots ,M$, $j=1,2$;

Part 1: Recovering of $f_{j,i_j}$, $i_j=1,2,\cdots ,N_j-1,j=1,2.$

$k\leftarrow 0$, $\delta \left(0\right)\leftarrow ϵ+1$;

while $\delta \left(k\right)>ϵ$do

$f_{j,i_j}=\left\{\begin{array}{c}{\displaystyle \frac{-{\Lambda }_{i_j}\left(m_j\right)}{g_j^{n^{*}}},\hspace{1em}\hspace{1em}\hspace{1em}k=0, i_j=1,2,\cdots ,N_j-1,}\hfill \\ {\displaystyle \frac{V_{j,t_{i_j}}^{n^{*}}-{\Lambda }_{i_j}\left(m_j\right)}{g_j^{n^{*}}},\hspace{1em}k>0i_j=1,2,\cdots ,N_j-1,;}\hfill \end{array}\right.$

Find $V_{j,i_j}^{(k+1),n+1}$, $i_j=0,1,\cdots ,N_j$, $n=0,1,\cdots ,n^{*}-1$, $j=1,2$, solving (34)–(36);

Determine $V_{j,t_{i_j}}^{n^{*}}$;

$\delta \left(k\right)\leftarrow max\{\parallel V_1^{(k+1),n+1}-V_1^{\left(k\right),n+1}\parallel ,\parallel V_2^{(k+1),n+1}-V_2^{\left(k\right),n+1}\parallel \}$;

$k\leftarrow k+1$.

end while

$V_{j,i_j}^{n+1}\leftarrow V_{j,i_j}^{(k+1),n+1}$, $i_j=0,1,\cdots ,N_j$, $n=0,1,\cdots ,n^{*}-1$, $j=1,2$;

Part 2: Determination of $U_{j,i_j}^n$, $i_j=0,1,\cdots ,N_j$, $j=1,2$, $n=1,2,\cdots ,M$.

if $t^{*}<T$ then,

solve (39) to find $U_{j,i_j}^n$, $n=1,2,\cdots ,n^{*}$, $i_j=0,1,\cdots ,N_j$, $j=1,2$;

solve (40) to find $U_{j,i_j}^n$, $n=n^{*}+1,n^{*}+2,\cdots ,M$, $i_j=0,1,\cdots ,N_j$, $j=1,2$;

else

if $t^{*}=T$ then

solve (39) to find $U_{j,i_j}^n$, $n=1,2,\cdots ,n^{*}=M$, $i_j=0,1,\cdots ,N_j$, $j=1,2$;

end if

end if

4.2. High-Order Numerical Schemes

Now, we construct high-order discretizations. To this end, first, we will clarify some issues.

(1) We will recover the functions $f_j$, $j=1,2$ also at the interface boundaries, since the values $f_1\left(b_1\right)$ and $f_2\left(a_2\right)$ are necessary for the numerical scheme.

(2) We need the values of ${\phi }_j$, $j=1,2$ at the interface boundaries. We suppose that the functions ${\phi }_j$ are given explicitly and ${\phi }_1\left(b_1\right)$, ${\phi }_2\left(a_2\right)$ can be determined directly. But, for example, in the case of table functions ${\phi }_j(x_{j,i_j})$, $i_j=1,2,\cdots ,N_j-1$, we can use (20) and (21) (compatibility conditions) to determine ${\phi }_1\left(b_1\right)$ and ${\phi }_2\left(a_2\right)$.

We implement the finite volume-difference method, implicit-explicit time-stepping in the domain ${\overline{\omega }}_{\tau }\times {\overline{\omega }}_{h1}$ and implicit time-stepping in ${\overline{\omega }}_{\tau }\times {\overline{\omega }}_{h2}$.

The discretization of (12) and (13) is obtained in the same manner as in the previous section. At the inner grid nodes, the numerical scheme is (34) and (35).

For the boundary conditions, we construct a second-order in space approximation. Moreover, in the second domain, we utilize a fully implicit scheme. Using (18)–(21), the boundary conditions (14)–(16) are discretized as follows:

(41)$\begin{array}{c}{V}_{1,0}^{n+1}={\mathcal{M}}_1^{n+1},n=0,1,\cdots ,n^{*}-1,\hfill \\ \left(\frac{1}{{\tau }^2}+\frac{2{\alpha }_1}{h_1}\right)V_{1,N_1}^{n+1}+\frac{2}{h_1}{p}_{1,N_1-1/2}{V}_{1,{\overline{x}}_{N_1}}^{n+1}=\frac{2{\beta }_1}{h_1}{V}_{2,0}^n+\frac{2}{h_1}{\Gamma }_1^{n+1}\hfill \\ +f_{1,N_1}{G}_1^{n+1}+\frac{2V_{1,N_1}^n-V_{1,N_1}^{n-1}}{{\tau }^2},n=1,2,\cdots ,n^{*}-1,\hfill \\ \left(\frac{2}{{\tau }^2}+\frac{2{\alpha }_1}{h_1}\right)V_{1,N_1}^1+\frac{2}{h_1}{p}_{1,N_1-1/2}{V}_{1,{\overline{x}}_{N_1}}^1=\frac{2{\beta }_1}{h_1}{\psi }_{2,0}^0+\frac{2}{h_1}{\Gamma }_1^1\hfill \\ +f_{1,N_1}{G}_1^1+\frac{2{\psi }_{1,N_1}}{{\tau }^2}-\frac{2(P_{1,N_1}+f_{1,N_1}{g}_1^0)}{\tau },\hfill \\ \left(\frac{1}{{\tau }^2}+\frac{2{\alpha }_2}{h_2}\right)V_{2,0}^{n+1}-\frac{2}{h_2}{p}_{2,1/2}{V}_{2,x_0}^{n+1}=\frac{2{\beta }_2}{h_2}{V}_{1,N_1}^n+\frac{2}{h_2}{\Gamma }_2^{n+1}+f_{2,0}^{n+1}{G}_2^{n+1}\hfill \\ +\frac{2V_{2,0}^n-V_{2,0}^{n-1}}{{\tau }^2},n=1,2,\cdots ,n^{*}-1,\hfill \\ \left(\frac{2}{{\tau }^2}+\frac{2{\alpha }_2}{h_2}\right)V_{2,0}^1-\frac{2}{h_2}{p}_{2,1/2}{V}_{2,x_0}^1=\frac{2{\beta }_2}{h_2}{\psi }_{1,N_1}+\frac{2}{h_2}{\Gamma }_2^1+f_{2,0}{G}_2^1\hfill \\ +\frac{2{\psi }_{2,0}}{{\tau }^2}-\frac{2(P_{2,0}+f_{2,0}{g}_2^0)}{\tau },\hfill \\ V_{2,N_2}^{n+1}={\mathcal{M}}_2^{n+1},n=0,1,\cdots ,n^{*}-1.\hfill \end{array}$

(42)$f_{j,i_j}=\frac{V_{j,t_{i_j}}^{n^{*}}-{\Lambda }_{i_j}\left(m_j\right)}{g_j^{n^{*}}},i_j=2-j,3-j,\cdots ,N_j+1-j,j=1,2,$

The values of ${\Lambda }_{N_1-1}\left(m_1\right)$ and ${\Lambda }_1\left(m_2\right)$ are obtained by (38). For the values ${\Lambda }_{N_1}\left(m_1\right)$ and ${\Lambda }_0\left(m_2\right)$, we use the interface conditions (3) and (4) to derive

$\begin{array}{c}{\Lambda }_{N_1}\left(m_1\right)=\frac{2}{h_1}\left(-{\alpha }_1{U}_{1,N_1}^{n^{*}}+{\beta }_1{U}_{2,0}^{n^{*}}+{\gamma }_1^{n^{*}}-p_{1,N_1-1/2}{m}_{1,{\overline{x}}_{N_1}}^n\right),\hfill \\ {\Lambda }_0\left(m_2\right)=\frac{2}{h_2}\left(p_{2,1/2}{m}_{2,x_0}+{\alpha }_2{U}_{2,0}^{n^{*}}-{\beta }_2{U}_{1,N_1}^{n^{*}}-{\gamma }_2^{n^{*}}\right).\hfill \end{array}$

As before, initiating an iteration process for (35), (35), (41) and (42), we compute $f_{j,i_j}$, $i_j=2-j,3-j,\cdots ,N_j+1-j$ and $U_{j,i_j}^n$, $n=1,2,\cdots ,n^{*}$, $i_j=0,1,\cdots ,N_j$, $j=1,2$, using (39). If $T>t^{*}$, then the solution $U_{j,i_j}^{n+1}$, $n=n^{*},n^{*}+1,\cdots ,M-1$, $i_j=0,1,\cdots ,N_1$, $j=1,2$ is obtained by the following numerical scheme:

(43)$\begin{array}{c}\hfill \begin{array}{c}{U}_{1,0}^{n+1}={\mu }_1^{n+1},\hfill \\ U_{1,\overline{t}{t}_{i_1}}^{n+1}-{\Lambda }_{i_1}\left(U_1^{n+1}\right)=f_{1,i_1}{g}_1^{n+1},i_1=1,2,\cdots ,N_1-1,\hfill \\ \left(\frac{1}{{\tau }^2}+\frac{2{\alpha }_1}{h_1}\right)V_{1,N_1}^{n+1}+\frac{2}{h_1}{p}_{1,N_1-1/2}{U}_{1,{\overline{x}}_{N_1}}^{n+1}=\frac{2{\beta }_1}{h_1}{U}_{2,0}^n+\frac{2}{h_1}{\gamma }_1^{n+1}\hfill \\ +f_{1,N_1}{g}_1^{n+1}+\frac{2U_{1,N_1}^n-U_{1,N_1}^{n-1}}{{\tau }^2},\hfill \\ \left(\frac{1}{{\tau }^2}+\frac{2{\alpha }_2}{h_2}\right)U_{2,0}^{n+1}-\frac{2}{h_2}{p}_{2,1/2}{U}_{2,x_0}^{n+1}=\frac{2{\beta }_2}{h_2}{U}_{1,N_1}^n+\frac{2}{h_2}{\gamma }_2^{n+1}+f_{2,0}^{n+1}{g}_2^{n+1},\hfill \\ +\frac{2V_{2,0}^n-V_{2,0}^{n-1}}{{\tau }^2},\hfill \\ U_{2,\overline{t}{t}_{i_2}}^{n+1}-{\Lambda }_{i_2}\left(U_2^{n+1}\right)=f_{2,i_2}{g}_2^{n+1},i_2=1,2,\cdots ,N_2-1,\hfill \\ U_{2,N_2}^{n+1}={\mu }_2^{n+1}.\hfill \end{array}\end{array}$

Algorithm 2 describes the main steps of the numerical recovering.

Algorithm 2: Recovering by $O\left(h^2\right)$ discretization

Require: $a_j$, $b_j$, T, $0<t^{*}\le T$, $m_j\left(x\right)$, ${\phi }_j\left(x\right)$, ${\psi }_j\left(x\right)$, ${\mu }_j\left(x\right)$, $p_j\left(x\right)$, $x\in {\Omega }_j$, ${\alpha }_j$, ${\beta }_j$,

${\gamma }_j\left(t\right)$, $g_j\left(t\right)$, $t\in [0,t^{*}]$, $ϵ$, $N_j$, M, $j=1,2$;

Ensure: $f_{j,i_j}$, $i_j=2-j,3-j,\cdots ,N_j+1-j$, $U_{j,i_j}^n$, $i_j=0,1,\cdots ,N_j$, $n=1,2,\cdots ,M$, $j=1,2$;

Part 1: Recovering of $f_{j,i_j}$, $i_j=2-j,3-j,\cdots ,N_j+1-j$, $j=1,2$.

$k\leftarrow 0$, $\delta \left(0\right)\leftarrow ϵ+1$;

while $\delta \left(k\right)>ϵ$ do

$f_{j,i_j}=\left\{\begin{array}{c}{\displaystyle \frac{-{\Lambda }_{i_j}\left(m_j\right)}{g_j^{n^{*}}},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}k=0, i_j=2-j,3-j,\cdots ,N_j+1-j,}\hfill \\ {\displaystyle \frac{V_{j,t_{i_j}}^{n^{*}}-{\Lambda }_{i_j}\left(m_j\right)}{g_j^{n^{*}}},\hspace{1em}\hspace{1em}k>0, i_j=2-j,3-j,\cdots ,N_j+1-j;}\hfill \end{array}\right.$

Find $V_{j,i_j}^{(k+1),n+1}$, $i_j=0,1,\cdots ,N_j$, $n=0,1,\cdots ,n^{*}-1$, $j=1,2$, solving (34), (35), (41);

Determine $V_{j,t_{i_j}}^{n^{*}}$;

$\delta \left(k\right)\leftarrow max\{\parallel V_1^{(k+1),n+1}-V_1^{\left(k\right),n+1}\parallel ,\parallel V_2^{(k+1),n+1}-V_2^{\left(k\right),n+1}\parallel \}$;

Solve (39) to find $U_{j,i_j}^{(k+1),n+1}$, $n=1,2,\cdots ,n^{*}-1$, $i_j=0,1,\cdots ,N_j$, $j=1,2$;

$k\leftarrow k+1$.

end while

$V_{j,i_j}^{n+1}\leftarrow V_{j,i_j}^{(k+1),n+1}$, $i_j=0,1,\cdots ,N_j$, $n=0,1,\cdots ,n^{*}-1$, $j=1,2$;

$U_{j,i_j}^{n+1}\leftarrow U_{j,i_j}^{(k+1),n+1}$, $i_j=0,1,\cdots ,N_j$, $n=0,1,\cdots ,n^{*}-1$, $j=1,2$;

Part 2: Determination of $U_{j,i_j}^n$, $i_j=0,1,\cdots ,N_j$, $j=1,2$, $n=n^{*}+1,n^{*}+2,\cdots ,M$.

if $t^{*}<T$ then,

solve (43) to find $U_{j,i_j}^n$, $n=n^{*}+1,n^{*}+2,\cdots ,M$, $i_j=0,1,\cdots ,N_j$, $j=1,2$;

end if

5. Numerical Simulations

In this section, we verify the efficiency of the developed Algorithms 1 and 2 for recovering the source terms $f_j$, $j=1,2$. We consider the problem (1)–(7) and (10) for

$\begin{array}{c}{\displaystyle a_1=1,b_1=2,a_2=3,b_2=5,T=1,t^{*}=\frac{1}{2},}\\ {\displaystyle {\alpha }_1=3,{\alpha }_2=1,{\beta }_1=2,{\beta }_2=0.5,g_1\left(t\right)=e^t,g_2\left(t\right)=e^{t/2}.}\end{array}$

To check the accuracy, we deal with the exact solution

$u_1(x,t)=e^{t}cos\frac{\pi x}{4},u_2(x,t)=e^{t/2}cos\frac{\pi x}{2}$

and, in this case, the right-hand sides $f_j$, ${\gamma }_j$, $j=1,2$ become

(44)$\begin{array}{c}\begin{array}{c}{f}_1\left(x\right)=cos(\pi x/4)+\frac{\pi }{4}\frac{d{p}_1\left(x\right)}{dx}sin\frac{\pi x}{4}+\frac{{\pi }^2}{16}{p}_1\left(x\right)cos\frac{\pi x}{4},\hfill \\ f_1\left(x\right)=\frac{1}{4}cos(\pi x/2)+\frac{\pi }{2}\frac{d{p}_2\left(x\right)}{dx}sin\frac{\pi x}{2}+\frac{{\pi }^2}{4}{p}_2\left(x\right)cos\frac{\pi x}{2}.\hfill \end{array}\hfill \\ {\gamma }_1\left(t\right)=-\frac{\pi }{4}{p}_1\left(b_1\right)e^{t}sin\frac{\pi b_1}{4}+{\alpha }_1{e}^{t}cos\frac{\pi b_1}{4}-{\beta }_1{e}^{t/2}cos\frac{\pi a_2}{2},\hfill \\ {\gamma }_2\left(t\right)=\frac{\pi }{2}{p}_2\left(a_2\right)e^{t/2}sin\frac{\pi a_2}{2}+{\alpha }_2{e}^{t/2}cos\frac{\pi a_2}{2}-{\beta }_2{e}^{t}cos\frac{\pi b_1}{4}.\hfill \end{array}$

We give errors in the maximal discrete norm (${\epsilon }_j$, ${\mathcal{E}}_j$) and order of convergence ($c{r}_j$, $C{R}_j$) of the determined functions $f_j$ and solutions $u_j$, obtained by the inverse problem, in comparison with the exact ones $f_j\left(x\right)$ and $u_j(x,t)$, $i=1,2$, respectively

$\begin{array}{c}{\displaystyle {\epsilon }_j={\epsilon }_j\left(N_j\right)\underset{0\le i_j\le N_j}{max}|f_j\left(t_n\right)-f_j^n|,{\mathcal{E}}_j={\mathcal{E}}_j\left(N\right)=\underset{0\le n\le M}{max}\underset{0\le i_j\le N_j}{max}|u_j(x_{i_j},t_n)-U_{j,i_j}^n|,}\\ {\displaystyle c{r}_j={log}_2\frac{{\epsilon }_i\left(N_j\right)}{{\epsilon }_i\left(2N_j\right)},C{R}_j={log}_2\frac{{\mathcal{E}}_j\left(N_j\right)}{{\mathcal{E}}_j\left(2N_j\right)}.}\end{array}$

We consider two test problems:

TP1: $p_1=2x^2+3$, $p_2=sin\left(\pi x\right)+3$,

TP2: $p_1=\left\{\begin{array}{cc}2,\hfill & x\le 1.5,\hfill \\ 3.5-x,\hfill & x>1.5,\hfill \end{array}\right.p_2=\left\{\begin{array}{cc}x+2,\hfill & x\le 4,\hfill \\ 6-x,\hfill & x>4.\hfill \end{array}\right.$

Let $N=N_1+N_2$. The computations are performed for $h=h_1=h_2$ ($N_2=2N_1$), and $ϵ={10}^{-6}$.

Example 1 

(Inverse problem: exact measurements). For the computations, we take measurements $m_j$ in (10) from the exact solution at $t=t^{*}$. In Table 1 and Table 2, we present the errors and order of convergence of the recovered functions $f_j$ and the numerical solution $u_j$, $j=1,2$, computed by Algorithm 1, $\tau =h$, for TP1 and TP2, respectively. Also, the number of iterations k is given. We observe that for both test examples, the convergence rate of the numerical $f_j$ and $u_j$, $j=1,2$ is first order. The number of iterations increases slightly as the mesh becomes finer, but this is not unexpected bearing in mind that we solve an inverse problem. The number of iterations required to determine $f_j$, $j=1,2$ in TP2 is two times greater than in TP1, due to the fact that, in contrast to TP1, in TP2, $f_1$ is a discontinuous function and $f_2$ is not smooth.

In Table 3 and Table 4, we give the results from the computations with Algorithm 2, $\tau =h^2$, for TP1 and TP2, respectively. The results illustrate that the convergence rate of the recovered functions $f_j$ and the numerical solution $u_j$, $j=1,2$ for TP1 is second, while for TP2, it is first. The number of iterations required to reach the desired accuracy is almost the same as for Algorithm 1.

In Figure 1, we plot the quantity $\delta \left(k\right)$ at each iteration for TP1 and TP2, computed by Algorithm 1, $N=240$. We observe that the process is convergent both for the smooth functions $f_j$, $j=1,2$ (TP1) and the non-smooth and discontinuous functions $f_j$, $j=1,2$, at the expense of a larger number of iterations in comparison with TP1.

Example 2 

(Inverse problem: noisy measurements). Now, we test the efficiency of the proposed methods for noisy observations

$m_{j,i_j}^{ϵ}=m_{j,i_j}+{\rho }_j{\sigma }_{j,i_j}{m}_{j,i_j},j=1,2,$

where ${\rho }_j$ is the noise level and ${\sigma }_{j,i_j}$, $j=1,2$ is a random function, uniformly distributed on the interval $[0,1]$. We smooth the measured data, applying polynomial curve fitting of degree 7.

The computations are performed for one and the same mesh parameters: $N=240$, $M=80$ ($\tau =h$). In Figure 2, Figure 3, Figure 4 and Figure 5, we plot the exact and recovered functions $f_j$, $j=1,2$, computed by Algorithms 1 and 2 for TP1 and TP2 and different levels of noise. We find that the recovering is successful with almost the same precision for both algorithms, with slightly better fitting observed for Algorithm 2. This reflects the numerical solution U, as shown on Figure 6 and Figure 7. The number of iterations is similar as in Example 1, i.e., 32 for TP1 and 62–63 for TP2 for both algorithms. This is illustrated in Figure 8.

6. Conclusions

In this work, we constructed a numerical method for the solution of an initial boundary value problem for a hyperbolic equation defined on disjoint intervals. It is a specific interface problem, in which the jump conditions are of the Robin type, while the external boundary conditions are of the Dirichlet type. However, the problem with Neumann external boundary conditions can be treated similarly. We present a method for reduction of the inverse problem for recovering the space right-hand sides to a direct one with a non-local-in-time differential operator and non-local initial conditions. Next, on this basis, we develop efficient iterative numerical approaches to solve the inverse problem.

The numerical results illustrate that for exact measurements, the spatial order of convergence is first or second depending on the method and the smoothness of the recovered functions. The convergence is attained at a moderate number of iterations. The experiments with perturbed data show that the right-hand-side functions are restored with optimal precision. The developed algorithms successfully recover both the smooth unknown functions, as well as the non-smooth and discontinuous right-hand side.

Although in the present paper, our attention is concentrated on the 1D case, it is interesting to address these issues to the corresponding 2D problem. This will be the subject of a future investigation.

The model hyperbolic problem (as well as the corresponding parabolic ones, see, e.g., [16]), in fact is a system of PDEs, coupled at the boundary. The method developed in the paper for studying the inverse problem could also be applied to systems obtained by splitting single equations, as for example in [34,35], where the splitting method is used for solving a Schrodinger-type equation.

Author Contributions

Conceptualization, M.N.K. and L.G.V.; methodology, M.N.K. and L.G.V.; investigation, M.N.K. and L.G.V.; resources, M.N.K. and L.G.V.; writing—original draft preparation, M.N.K. and L.G.V.; writing—review and editing, L.G.V.; validation, M.N.K. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Acknowledgments

The authors are very grateful to the anonymous reviewers, whose valuable comments and suggestions improved the quality of this paper.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The values of [Forumla omitted. See PDF.] at each iteration, Algorithm 1, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], TP1 (left) and TP2 (right), Example 1.

Figure 2. Exact (solid red line) and recovered [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Algorithm 1 (dash-dotted blue line) and Algorithm 2 (dashed green line), [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], TP1, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Example 2.

Figure 3. Exact (solid red line) and recovered [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Algorithm 1 (dash-dotted blue line) and Algorithm 2 (dashed green line), TP1, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Example 2.

Figure 4. Exact (solid red line) and recovered [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Algorithm 1 (dash-dotted blue line) and Algorithm 2 (dashed green line), TP2, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Example 2.

Figure 5. Exact (solid red line) and recovered [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Algorithm 1 (dash-dotted blue line) and Algorithm 2 (dashed green line), TP2, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Example 2.

Figure 6. Error [Forumla omitted. See PDF.]. in the whole space-time domain for TP1, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Algorithm 1 (left) and Algorithm 2 (right), [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Example 2.

Figure 7. Error [Forumla omitted. See PDF.] in the whole space-time domain for TP1, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Algorithm 1 (left) and Algorithm 2 (right), [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Example 2.

Figure 8. The values of [Forumla omitted. See PDF.] at each iteration, TP1, Algorithm 2, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], (left) and TP2, Algorithm 2, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], (right), [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], Example 2.

Table 1

Errors and convergence rate of $f_j$ and $u_j$, $j=1,2$, Algorithm 1, $\tau =h$, TP1, Example 1.

${\mathit{N}}_1$	${\mathit{\epsilon }}_1$	${\mathit{cr}}_1$	${\mathit{\epsilon }}_2$	${\mathit{cr}}_2$	${\mathcal{E}}_1$	${\mathit{CR}}_1$	${\mathcal{E}}_2$	${\mathit{CR}}_2$	k
60	8.766 × ${10}^{-3}$		3.659 × ${10}^{-3}$		2.555 × ${10}^{-4}$		2.312 × ${10}^{-3}$		30
120	4.576 × ${10}^{-3}$	0.938	1.602 × ${10}^{-3}$	1.191	2.394 × ${10}^{-4}$	0.094	6.053 × ${10}^{-4}$	1.933	31
240	2.542 × ${10}^{-3}$	0.848	7.529 × ${10}^{-4}$	1.089	1.644 × ${10}^{-4}$	0.542	1.674 × ${10}^{-4}$	1.854	32
480	1.394 × ${10}^{-3}$	0.867	3.679 × ${10}^{-4}$	1.033	9.638 × ${10}^{-5}$	0.771	5.121 × ${10}^{-5}$	1.709	32
960	7.379 × ${10}^{-4}$	0.918	1.838 × ${10}^{-4}$	1.001	5.206 × ${10}^{-5}$	0.889	1.840 × ${10}^{-5}$	1.477	32

Table 2

Errors and convergence rate of $f_j$ and $u_j$, $j=1,2$, Algorithm 1, $\tau =h$, TP2, Example 1.

${\mathit{N}}_1$	${\mathit{\epsilon }}_1$	${\mathit{cr}}_1$	${\mathit{\epsilon }}_2$	${\mathit{cr}}_2$	${\mathcal{E}}_1$	${\mathit{CR}}_1$	${\mathcal{E}}_2$	${\mathit{CR}}_2$	k
60	1.739 × ${10}^{-2}$		8.355 × ${10}^{-3}$		1.002 × ${10}^{-2}$		3.858 × ${10}^{-3}$		53
120	8.829 × ${10}^{-3}$	0.978	4.904 × ${10}^{-3}$	0.769	5.123 × ${10}^{-3}$	0.968	2.106 × ${10}^{-3}$	0.873	57
240	4.531 × ${10}^{-3}$	0.963	2.675 × ${10}^{-3}$	0.874	2.595 × ${10}^{-3}$	0.982	1.110 × ${10}^{-3}$	0.924	60
480	2.337 × ${10}^{-3}$	0.955	1.404 × ${10}^{-3}$	0.930	1.308 × ${10}^{-3}$	0.988	5.726 × ${10}^{-4}$	0.955	61
960	1.197 × ${10}^{-3}$	0.965	7.226 × ${10}^{-4}$	0.958	6.572 × ${10}^{-4}$	0.993	2.925 × ${10}^{-4}$	0.969	62

Table 3

Errors and convergence rate of $f_j$ and $u_j$, $j=1,2$, Algorithm 2, $\tau =h^2$, TP1, Example 1.

N	${\mathit{\epsilon }}_1$	${\mathit{cr}}_1$	${\mathit{\epsilon }}_2$	${\mathit{cr}}_2$	${\mathcal{E}}_1$	${\mathit{CR}}_1$	${\mathcal{E}}_2$	${\mathit{CR}}_2$	k
15	2.164 × ${10}^{-2}$		3.225 × ${10}^{-2}$		1.337 × ${10}^{-2}$		4.889 × ${10}^{-2}$		29
30	4.259 × ${10}^{-3}$	2.345	5.563 × ${10}^{-3}$	2.536	2.484 × ${10}^{-3}$	2.428	8.482 × ${10}^{-3}$	2.527	32
60	1.285 × ${10}^{-3}$	1.729	1.507 × ${10}^{-3}$	1.884	6.317 × ${10}^{-4}$	1.975	2.137 × ${10}^{-3}$	1.989	32
120	3.411 × ${10}^{-4}$	1.914	3.875 × ${10}^{-4}$	1.960	1.587 × ${10}^{-4}$	1.993	5.360 × ${10}^{-4}$	1.995	33
240	8.783 × ${10}^{-5}$	1.957	9.934 × ${10}^{-5}$	1.964	3.999 × ${10}^{-5}$	1.989	1.351 × ${10}^{-4}$	1.988	33

Table 4

Errors and convergence rate of $f_j$ and $u_j$, $j=1,2$, Algorithm 2, $\tau =h^2$, TP2, Example 1.

N	${\mathit{\epsilon }}_1$	${\mathit{cr}}_1$	${\mathit{\epsilon }}_2$	${\mathit{cr}}_2$	${\mathcal{E}}_1$	${\mathit{CR}}_1$	${\mathcal{E}}_2$	${\mathit{CR}}_2$	k
15	5.576 × ${10}^{-2}$		3.942 × ${10}^{-2}$		6.642 × ${10}^{-2}$		6.515 × ${10}^{-2}$		53
30	2.346 × ${10}^{-2}$	1.249	5.624 × ${10}^{-3}$	2.809	2.783 × ${10}^{-2}$	1.255	9.944 × ${10}^{-3}$	2.712	61
60	1.262 × ${10}^{-2}$	0.895	2.308 × ${10}^{-3}$	1.285	1.386 × ${10}^{-2}$	1.006	2.689 × ${10}^{-3}$	1.887	62
120	6.729 × ${10}^{-3}$	0.907	1.374 × ${10}^{-3}$	0.749	6.901 × ${10}^{-3}$	1.006	1.259 × ${10}^{-3}$	1.095	62
240	3.543 × ${10}^{-3}$	0.925	7.326 × ${10}^{-4}$	0.907	3.442 × ${10}^{-3}$	1.003	6.083 × ${10}^{-4}$	1.049	62

References

1. Caffarelli, L. A monotonicity formula for heat functions in disjoint domains. Bound. Value Probl. Partial. Differ. Equ. Appl.; 1993; 29, pp. 53-60.

2. Datta, A.K. Biological and Bioenvironmental Heat and Mass Transfer; 1st ed. Marcel Dekker: New York, NY, USA, 2002; 424p.

3. Govoli, D. Exact representations on artificial interfaces and applications in mechanics. Appl. Mech. Rev.; 1999; 52, pp. 333-349. [DOI: https://dx.doi.org/10.1115/1.3098920]

4. Govoli, D. Finite element modeling of thin layers. Comput. Model. Eng. Sci.; 2004; 5, pp. 497-514.

5. Stefanov, P.; Uhlmann, G.; Vasy, A. The transmission problem in linear isotropic elasticity. Pure Appl. Anal.; 2021; 3, pp. 109-161. [DOI: https://dx.doi.org/10.2140/paa.2021.3.109]

6. Fatori, L.H.; Lueders, E.; Munoz Rivera, J.E. Tramsmission problem for hyperbolic thermoelastic systems. J. Therm. Stress.; 2003; 26, pp. 739-763. [DOI: https://dx.doi.org/10.1080/713855994]

7. Bouziani, A. Solution of a transmission problem for semilinear parabolic-hyperbolic equations by the time-discretization. Int. J. Stoch. Anal.; 2006; 2006, 061439. [DOI: https://dx.doi.org/10.1155/JAMSA/2006/61439]

8. Jovanović, B.S.; Vulkov, L.G. Numerical solution of a hyperbolic transmission problem. Comput. Methods Appl. Math.; 2008; 8, pp. 374-385. [DOI: https://dx.doi.org/10.2478/cmam-2008-0027]

9. Jovanović, B.S.; Vulkov, L.G. Numerical solution of a two-dimensional hyperbolic transmission problem. J. Comput. Appl. Math.; 2010; 235, pp. 519-534. [DOI: https://dx.doi.org/10.1016/j.cam.2010.05.045]

10. Hasanoglu, A.; Romanov, V.G. Introduction to Inverse Problems for Differential Equations; 1st ed. Springer: Cham, Switzerland, 2017; 261p.

11. Kabanikhin, S.I. Inverse and Ill-Posed Problems; DeGruyer: Berlin, Germany, 2011.

12. Lesnic, D. Inverse Problems with Applications in Science and Engineering; CRC Press: Abingdon, UK, 2021; 349.

13. Romanov, V.; Hasanov, A. Uniqueness and stability analysis of final data inverse source problems for evolution equations. J. Inverse Ill-Posed Probl.; 2022; 30, pp. 425-446. [DOI: https://dx.doi.org/10.1515/jiip-2021-0072]

14. Samarskii, A.A.; Vabishchevich, P.N. Numerical Methods for Solving Inverse Problems of Mathematical Physics; de Gruyter: Berlin, Germany, 2007; 452p.

15. Chapko, R.; Mindrinos, L. On the numerical solution of a hyperbolic inverse boundary value problem in bounded domains. Mathematics; 2022; 10, 750. [DOI: https://dx.doi.org/10.3390/math10050750]

16. Koleva, M.N.; Vulkov, L.G. Numerical solution of the retrospective inverse parabolic problem on disjoint intervals. Computation; 2023; 11, 204. [DOI: https://dx.doi.org/10.3390/computation11100204]

17. Sidi, H.O.; Zaky, M.A.; Qiu, W.; Hendy, A.S. Identification of an unknown spatial source function in a multidimensional hyperbolic partial differential equation with interior degeneracy. Appl. Numer. Math.; 2023; 192, pp. 1-18. [DOI: https://dx.doi.org/10.1016/j.apnum.2023.05.021]

18. Van Bockstal, K.; Slodička, M. Recovery of a time-dependent heat source in one-dimensional thermoelasticity of type-III. Inverse Probl. Sci. Eng.; 2017; 25, pp. 749-770. [DOI: https://dx.doi.org/10.1080/17415977.2016.1199696]

19. Van Bockstal, K.; Slodička, M. Recovery of a space-dependent vector source in thermoelastic systems. Inverse Probl. Sci. Eng.; 2015; 23, pp. 956-968. [DOI: https://dx.doi.org/10.1080/17415977.2014.959008]

20. Koleva, M.N.; Vulkov, L.G. Numerical identification of external boundary conditions for time fractional parabolic equations on disjoint domains. Fractal Fract.; 2023; 7, 326. [DOI: https://dx.doi.org/10.3390/fractalfract7040326]

21. Nursultanov, M.; Oksanen, L. Introduction to inverse problem for hyperbolic PDEs. arXiv; 2023; arXiv: 2311.06020

22. Van Bockstal, K. Identification of an unknown spatial load distribution in a vibrating beam or plate from the final state. J. Inverse Ill-Posed Probl.; 2019; 27, pp. 623-642. [DOI: https://dx.doi.org/10.1515/jiip-2018-0068]

23. Beilina, L.; Cristofol, M.; Li, S.; Yamamoto, M. Lipschitz stability for an inverse hyperbolic problem of determining two coefficients by a finite number of observations. Inverse Probl.; 2018; 34, 015001. [DOI: https://dx.doi.org/10.1088/1361-6420/aa941d]

24. Chadam, J.M.; Yin, H.-M. Determination of an unknown function in a parabolic equation with an overspecified condition. Numer. Methods Appl. Sci.; 1990; 13, pp. 421-430. [DOI: https://dx.doi.org/10.1002/mma.1670130506]

25. Wloka, J. Partial Differential Equations; Cambridge University Press: Cambridge, UK, 1987.

26. Bouziani, A. Application of Rothe’s method to a semilinear hyperbolic equation. Georgian Math. J.; 2010; 17, pp. 437-458. [DOI: https://dx.doi.org/10.1515/gmj.2010.020]

27. Dehghan, M. Numerical schemes for one-dimensional parabolic equations with nonstandard initial condition. Appl. Math. Comput.; 2004; 147, pp. 321-331. [DOI: https://dx.doi.org/10.1016/S0096-3003(02)00669-0]

28. Dehghan, M. Implicit collocation technique for heat equation with non-classic initial condition. Int. J. Nonlinear Sci. Numer. Simul.; 2006; 7, pp. 461-466. [DOI: https://dx.doi.org/10.1515/IJNSNS.2006.7.4.461]

29. Lin, Y. Analytical and numerical solutions for a class nonlocal nonlinear parabolic differential equations. SIAM J. Math. Anal.; 1994; 25, pp. 1577-1594.

30. Lin, Y. Finite-difference solutions for parabolic equations with time weighting initial conditions. Appl. Math. Comput.; 1994; 65, pp. 49-61. [DOI: https://dx.doi.org/10.1016/0096-3003(94)90165-1]

31. Martin-Vaquero, J.; Sajvicius, S. The two level finite difference scheme for the heat equation with nonlocal initial condition. Appl. Math. Comput.; 2019; 342, pp. 160-177. [DOI: https://dx.doi.org/10.1016/j.amc.2018.09.025]

32. Atanasov, A.Z.; Koleva, M.N.; Vulkov, L.G. Inverse problem numerical analysis of forager bee losses in spatial environment without contamination. Symmetry; 2023; 15, 2099. [DOI: https://dx.doi.org/10.3390/sym15122099]

33. Samarskii, A.A. The Theory of Difference Schemes; Marcel Dekker: New York, NY, USA, 2001; 761p.

34. Castro López, R.; Sun, G.-H.; Camacho-Nieto, O.; Yáñez-Márquez, C.; Dong, S.-H. Analytical traveling-wave solutions to a generalized Gross–Pitaevskii equation with some new time and space varying nonlinearity coefficients and external fields. Phys. Lett. A; 2017; 381, pp. 2978-2985. [DOI: https://dx.doi.org/10.1016/j.physleta.2017.07.012]

35. Guo, Y.-S.; Li, W.; Dong, S.-H. Gaussian solitary solution for a class of logarithmic nonlinear Schrödinger equation in (1 + n) dimensions. Results Phys.; 2023; 44, 106187. [DOI: https://dx.doi.org/10.1016/j.rinp.2022.106187]