Numerical Solution of External Boundary Conditions Inverse Multilayer Diffusion Problems

Abstract

The present study is concerned with the numerical solution of external boundary conditions in inverse problems for one-dimensional multilayer diffusion, using the difference method. First, we formulate multispecies parabolic problems with three types of Dirichlet–Neumann–Robin internal boundary conditions that apply at the interfaces between adjacent layers. Then, using the symmetry of the diffusion operator, we prove the well-posedness of the direct (forward) problem in which the coefficients, the right-hand side, and the initial and boundary conditions are given. In inverse problems, instead of external boundary conditions of the first and the last layers, point observations of the solution within the entire domain are posed. Rothe’s semi-discretization of differential problems combined with a symmetric exponential finite difference solution for elliptic problems on each time layer is proposed to develop an efficient semi-analytical approach. Next, using special solution decomposition techniques, we numerically solve the inverse problems for the identification of external boundary conditions. Numerical test examples are discussed.

1. Introduction

Many problems in industries, environmental science, biology, and applied mechanics are modeled by diffusion processes in layered materials. Examples include the transport of solutes, chemicals, and gases in porous media; heat and mass transfer in composite materials; geological profiles; oceanography; meteorology; drug-delivery effectiveness in living tissues; and the study of biological tissues using infrared light [1,2,3,4,5,6].

Analytical solutions to such problems are of significant interest, although their construction is challenging. Nonetheless, they are highly valuable as they offer insight into solution behavior and serve as a benchmark for numerical methods.

The method of separate solutions in each layer is frequently employed for constructing analytical solutions. This approach often results in finding eigenfunctions that are coupled through internal boundary conditions at the interfaces between adjacent layers. However, it proves to be ineffective for differential equations with variable coefficients, as discussed in [1,2,7].

The study in [5] provides an exact stability analysis of a multilayer diffusion–reaction problem, accommodating a potentially large number of layers. An analytical approach based on Laplace transforms for solving general 2D multilayer diffusion problems with different convection heat transfer coefficients in each layer is presented in [8]. The stability of the 2D multilayer diffusion–reaction system is investigated in [9].

In [10], the authors propose exact solutions for diffusive transport models in domains experiencing spatially nonuniform growth. An embedding method is developed in [11] to analytically solve multilayer diffusion problems with time-dependent inhomogeneous boundary conditions. Applying the Laplace transform, semi-analytical solutions for axisymmetric reaction–diffusion equations with a piecewise constant coefficient and continuous transient boundary conditions have been developed in [12]. Semi-analytical solutions for advection–dispersion equations in multilayer porous media with Robin boundary conditions, as well as continuous and general interface conditions, are presented in [1,2]. The exact representation of the solution on the interface is introduced in [3].

The finite difference method is widely recognized for solving diffusion problems in single-layer systems (for example, [13,14,15,16,17,18]). It has been effectively used to address parabolic problems on disjoint domains (for example, [19,20,21]) and interface time-fractional diffusion equations (see [22]). Additionally, it has been applied to multilayer diffusion problems with both continuous and jump conditions at interfaces [4]. A second-order finite difference scheme for elliptic systems with perfect and imperfect contact at interfaces has been developed in studies such as [23,24].

The finite volume-difference method is developed in [6] for solving multilayer diffusion problems with different interface conditions, namely perfect and imperfect contact, as well as continuous and more general diffusive flux across the interface.

In [25], the authors provide theoretical and numerical investigation for solving a stationary diffusion problem in a domain with two distinct material constituents and significant differences in conductivity.

One effective method for solving interface problems is the immersed interface method (IIM). This approach achieves second-order accuracy and does not necessitate that grid points be aligned with the interfaces [26]. A diffusion equation with nonlinear localized chemical reactions is solved numerically in [27] using IIM.

Exponential difference schemes are highly effective in resolving boundary and interior layers of small diffusion problems [28]. In [29], the authors construct exponential difference schemes for solving convection–diffusion boundary value problems in the multidimensional case. Exponential difference schemes have been developed and analyzed in [30] for addressing nonlinear problems arising in variably saturated flow with discontinuous absolute permeability.

Numerous researchers (see e.g., [14,18,31,32,33]) have investigated the inverse problems related to determining unknown boundary conditions in parabolic equations, because of their variety of applications, especially in heat equations on a single domain. From a mathematical perspective, these inverse problems are classified as ill-posed problems [13,16,31,34,35,36], meaning that small changes in the input data can result in significant errors in the output. Furthermore, deriving analytical solutions for inverse problems involving multiple layers is often infeasible. Therefore, developing effective numerical methods for addressing such inverse problems is crucial.

The existence and uniqueness of the solution of inverse problems for one-layer parabolic and fractional diffusion equations in which a boundary condition function is unknown upon fixed-point measurements were investigated in [20,32,33,37].

In [14], a semigroup method was developed to analyze inverse problems for reconstruction of the right boundary condition in parabolic problems on the base of point observation. Numerical methods for the determination of BCs in parabolic, time-fractional parabolic, and parabolic–hyperbolic problems on disjoint domains were developed and analyzed in [20,21,38].

In [7], the authors solved the inverse heat conduction problem in layered materials by representing the numerical solution as a series of exponential basis functions defined in space and time. A numerical approach utilizing the method of fundamental solutions is presented in [39] for determining a moving boundary from Cauchy data within a one-dimensional heat equation involving a multilayer domain.

In this work, we prove the well-posedness of the direct multilayer interface diffusion problem and propose a numerical method for solving both direct and inverse problems to identify external time-dependent boundary conditions based on additional observations within different layers. The exponential finite difference scheme, extended for layered problems, is employed (see [30]). Boundary conditions are discretized by developing a new, second-order approximation. Additionally, a decomposition technique is applied to transform the ill-posed problem into a series of well-posed problems. Correct formulas are proposed for the computation of the solution to the inverse problem. This summarizes our main scientific contribution and its practical significance.

The structure of the paper is as follows: Section 2 formulates the multilayer interface problem along with the interface conditions. Section 3 addresses the existence and uniqueness of the solution within specific Sobolev spaces. In Section 4, exponential difference schemes for multilayer initial boundary value problems are developed. Section 5 formulates and numerically solves inverse problems for unknown external boundary conditions. Section 6 discusses computational results, and the paper concludes with final remarks and future perspectives.

2. Multilayer Diffusion Problems

This section provides a brief overview of the internal boundary conditions applied at the interfaces between adjacent layers.

We consider a diffusion process on the interval $[l_0,L]$, divided into M subintervals, as seen in the studies such as [4,5,6]. The partitioning is defined by $l_0<l_1<\dots <l_{M-1}$$<l_M=L$. For each subinterval (layer), ${\mathsf{\Omega }}_m=(l_{m-1},l_m)$, $m=1,2,\dots ,M$, the differential equation is given by [2,7]

$${\displaystyle \frac{\partial u_m}{\partial t}}=d_m{\displaystyle \frac{{\partial }^2{u}_m}{\partial x^2}}+f(x,t),x\in (l_{m-1},l_m),0<t<T,$$

(1)

where $u_m(x,t)$ represents the solution (for instance, temperature or concentration [40]) at position x and time t in the m-th layer, and $d_m>0$ is the diffusion coefficient in that layer.

The initial conditions are given as follows:

$$u_m(x,0)=u_m^0\left(x\right),m=1,\dots ,M$$

(2)

and external BCs are defined as follows:

$$\begin{array}{ccc}& & a_0{u}_1(l_0,t)-b_0{\displaystyle \frac{\partial u_1}{\partial x}}(l_0,t)=g_0\left(t\right),0<t<T,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & a_L{u}_M(L,t)+b_L{\displaystyle \frac{\partial u_M}{\partial x}}(L,t)=g_L\left(t\right),0<t<T,\hfill \end{array}$$

(4)

where the coefficients $a_0,b_0,a_L,b_L$ are non-negative constants satisfying $a_0+b_0>0$ and $a_L+b_L>0$, and $g_0\left(t\right)$ and $g_L\left(t\right)$ are specified time-dependent functions.

The internal boundary conditions at the interfaces between adjacent layers depend on the physical processes being examined. A review of various interface conditions is provided in [2]. In this context, we focus on some common conditions typically encountered in physics and mechanics, particularly in heat and mass transfer applications (for example, see [3,41]).

Typically, the continuity of the solution and the diffusive flux is implicitly assumed at each interface, that is,

$$\begin{array}{ccc}\hfill & & {\left[u_m\right]}_{l_m}=u_{m+1}(l_m,t)-u_m(l_m,t)=0,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& & {\left[{\displaystyle \frac{\partial u_m}{\partial x}}\right]}_{l_m}=d_{m+1}{\displaystyle \frac{\partial u_{m+1}}{\partial x}}(l_m,t)-d_m{\displaystyle \frac{\partial u_m}{\partial x}}(l_m,t)=0,\hfill \end{array}$$

(6)

for $t\in (0,T]$ and for all $m=1,\dots ,M-1$.

However, in fields such as mechanics, physics, and biology, many real processes require relations that differ from those in Equations (5) and (6).

Perfect contact conditions are

$$\begin{array}{ccc}& & {\left[u_m\right]}_{l_m}=u_{m+1}(l_m,t)-u_m(l_m,t)=0,0<t\le T,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}& & {\left[{\displaystyle \frac{\partial u_m}{\partial x}}\right]}_{l_m}={\gamma }_{m+1}{\displaystyle \frac{\partial u_{m+1}}{\partial x}}(l_m,t)-{\gamma }_m{\displaystyle \frac{\partial u_m}{\partial x}}(l_m,t)=0,0<t\le T,\hfill \end{array}$$

(8)

for ${\gamma }_m>0$. These interface conditions generalize Equations (5) and (6). For instance, in heat transfer, the thermal diffusivity is defined as $d_m:=k_m/\left({\rho }_m{c}_{pm}\right)$, where $k_m$ is the thermal conductivity and ${\rho }_m{c}_{pm}$ represents the volumetric capacity in layer m.

The concentrated reaction interface conditions, as described in studies such as [27], are utilized in modeling processes in chemistry, mechanics, physics, and other fields. These conditions are expressed as follows:

$$\begin{array}{ccc}& & {\left[u_m\right]}_{l_m}=u_{m+1}(l_m,t)-u_m(l_m,t)=0,\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& & {\left[{\displaystyle \frac{\partial u_m}{\partial x}}\right]}_{l_m}=r_{m+1}\left(t\right){\displaystyle \frac{\partial u_{m+1}}{\partial x}}(l_m,t)-r_m\left(t\right){\displaystyle \frac{\partial u_m}{\partial x}}(l_m,t)=u_m(l_m,t)+g_m\left(t\right),\hfill \end{array}$$

(10)

where $t\in (0,T]$ and $m=1,2,\dots ,M-1$.

It is important to note that the first two interface conditions (5)–(8) are homogeneous. However, in nonhomogeneous cases, a linear transformation can often be found that reduces the nonhomogeneous interface conditions to homogeneous ones. For instance, the nonhomogeneous boundary conditions (3) and (4) in conjunction with (5) and (6) for $M=2$ are addressed in [42]. Our primary focus is on the problem involving the interface conditions (9) and (10). We will demonstrate this reduction for the problem defined by Equations (1)–(4), (9) and (10) in the stationary case.

For simplicity, we take $l_0=-1,l_1=0,l_2=1$ and seek the solution in the form

$$u=(u_1,u_2),v=(v_1,v_2),w=(w_1,w_2),$$

where

$$v\left(x\right)=\left\{\begin{array}{cc}{v}_1=g_0+ax,\hfill & x\in (-1,0),\hfill \\ v_2=g_2+(1-x)b,\hfill & x\in (0,1).\hfill \end{array}\right.$$

Using the interface relations, we find

$$a={\displaystyle \frac{r_2(g_2-g_0)-g_1}{1+r_1+r_2}},b={\displaystyle \frac{(1+r_1)(g_0-g_2)-g_1}{1+r_1+r_2}}.$$

The constructed functions have the following properties:

$$\begin{array}{ccc}& & v_1(-1)=g_0,v_2\left(1\right)=g_2\hfill \\ & & v_2\left(0\right)=v_1\left(0\right)=0,r_2{\displaystyle \frac{d{v}_2}{dx}}\left(0\right)-r_1{\displaystyle \frac{d{v}_1}{dx}}\left(0\right)=u_1\left(0\right)+g_1,\hfill \\ & & -{\displaystyle \frac{d}{dx}}\left(d_1{\displaystyle \frac{d{v}_1}{dx}}\right)+e_1\left(x\right)(g_0+ax)=z_1\left(x\right),x\in (0,1),\hfill \\ & & -{\displaystyle \frac{d}{dx}}\left(d_2{\displaystyle \frac{d{v}_2}{dx}}\right)+e_2\left(x\right)(g_2+(1-x)b)=z_2\left(x\right),x\in (0,1),\hfill \end{array}$$

where

$$\begin{array}{ccc}& & z_1\left(x\right)=-a{\displaystyle \frac{d{d}_1}{dx}}{e}_1\left(x\right)(g_0+ax),\hfill \\ & & z_2\left(x\right)=-b{\displaystyle \frac{d{d}_2}{dx}}{e}_2\left(x\right)(g_2+(1-x)b).\hfill \end{array}$$

We define the functions $w_1$ and $w_2$ with zero boundary conditions and homogeneous interface conditions as follows:

$$\begin{array}{ccc}& & w_1(-1)=0,w_2\left(1\right)=0,\hfill \\ & & w_2\left(0\right)=w_1\left(0\right),r_2{\displaystyle \frac{d{w}_2}{dx}}\left(0\right)-r_1{\displaystyle \frac{d{w}_1}{dx}}\left(0\right)=0,\hfill \\ & & -{\displaystyle \frac{d}{dx}}\left(d_1{\displaystyle \frac{d{w}_1}{dx}}\right)+e_1\left(x\right)w_1=f_1-z_1,\hfill \\ & & -{\displaystyle \frac{d}{dx}}\left(d_2{\displaystyle \frac{d{w}_2}{dx}}\right)+e_2\left(x\right)w_2=f_2-z_2.\hfill \end{array}$$

3. Existence and Uniqueness of Solutions to Forward Problems

We now analyze the existence and uniqueness of weak solutions to forward problems defined in the previous section.

For definiteness, we shall consider the case of $M=2$ and the interface problem with relations (9) and (10), using the Dirac delta function as follows:

$${\displaystyle \frac{\partial u}{\partial t}}-{\displaystyle \frac{\partial }{\partial x}}\left(d\left(x\right){\displaystyle \frac{\partial u}{\partial x}}\right)+\left(e\left(x\right)+P\delta \left(x\right)\right)u=g(x,t),(x,t)\in Q_T,$$

(11)

where $P>0$ is a real constant, $Q_T=\mathsf{\Omega }\times (0,T)$, $\mathsf{\Omega }={\mathsf{\Omega }}_1\cup {\mathsf{\Omega }}_2$, ${\mathsf{\Omega }}_1=(-1,0)$, ${\mathsf{\Omega }}_2=(0,1)$ and the initial boundary conditions are

$$u(x,0)=u^0\left(x\right),u(-1,t)=u_l\left(t\right),u(1,t)=u_r\left(t\right).$$

(12)

A weak formulation is obtained by multiplying (11) by test functions $v\in L^2(0,T;H^1\left(\mathsf{\Omega }\right))$, integrating over $\mathsf{\Omega }$ and using integration by parts.

In this formulation, we seek $u\in L^2(0,T;H^1\left(\mathsf{\Omega }\right))$ with ${\displaystyle \frac{\partial u}{\partial t}\in L^2(0,T;H^{-1}\left(\mathsf{\Omega }\right))}$ and $g\in L^2\left(Q_T\right)$, such that

$$\begin{array}{ccc}\hfill & & \underset{0}{\overset{T}{\int }}\left[-\underset{-1}{\overset{1}{\int }}u{\displaystyle \frac{\partial v}{\partial t}}dx+{\int }_{-1}^1\left(d{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial v}{\partial x}}+e\left(x\right)uv\right)dx+Pu(0,t)v(0,t)-\underset{-1}{\overset{1}{\int }}gvdx\right]dt\hfill \\ & & =\underset{-1}{\overset{1}{\int }}{u}^0\left(x\right)v(x,0)dx.\hfill \end{array}$$

Suppose that such a weak solution exists and that ${\displaystyle \frac{\partial u}{\partial t}\in L_2((0,T);\mathsf{\Omega })}$. Taking $v(x,t)=w\left(x\right)\phi \left(t\right)$, where $w\in H^1\left(\mathsf{\Omega }\right)$, $w(-1)=u_l$, $u\left(1\right)=u_r$, ${\displaystyle \frac{d\phi }{dt}\in L_2(0,T),\phi \left(T\right)=0}$ and performing integration by parts, we obtain

$$\begin{array}{ccc}& & \underset{0}{\overset{T}{\int }}\left[\underset{-1}{\overset{1}{\int }}{\displaystyle \frac{\partial u}{\partial t}}wdx+\underset{-1}{\overset{1}{\int }}\left(d{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial w}{\partial x}}+e\left(x\right)uw\right)dx+Pu\left(0\right)w\left(0\right)-{\int }_{-1}^{1}gwdx\right]dt\hfill \\ & & +\phi \left(0\right){\int }_{-1}^1(u(x,0)-u^0\left(x\right))wdx=0.\hfill \end{array}$$

Hence, taking into account the arbitrariness of the function $\phi \left(t\right)$, it follows that

$${\int }_{-1}^1{\displaystyle \frac{\partial u}{\partial t}}wdx+{\int }_{-1}^1\left(d{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial w}{\partial x}}+e\left(x\right)uw\right)dx+Pu\left(0\right)w\left(0\right)={\int }_{-1}^{1}gwdx,$$

for almost $t\in (0,T)$ and

$$\underset{-1}{\overset{1}{\int }}u(x,0)w\left(x\right)dx=\underset{-1}{\overset{1}{\int }}{u}_0\left(x\right)w\left(x\right)dx.$$

Further, we shall work with the weak solution of the problem in (11) and (12) satisfying the last two identities.

Theorem 1.

Suppose that $a,e\in L_{\infty }\left(\mathsf{\Omega }\right)$ and the conditions

$$d_1\ge d\left(x\right)\ge d_0>0,e_1\ge e\left(x\right)\ge e_0>0forx\in \overline{\mathsf{\Omega }},$$

(13)

are satisfied and $g\in L^2\left(Q_T\right)$, $u_0\in H^1\left(\mathsf{\Omega }\right)$. Then, the problem in (11) and (12) admits a weak solution $u\in H^1\left(Q_T\right)$.

Proof. 

We follow the scheme of [43] (Section 3.2).

• Step 1. Construction of the difference scheme with respect to the time variable t to obtain the approximating solutions:

For any positive integer k and a function $w(x,t)$, denote $w^{k,j}\left(x\right)=w(x,jh),0,1,\dots ,k,$ where $\tau =T/m$.

Consider the approximation of Equation (11):

$${\displaystyle \frac{u^{m,j}-u^{m,j-1}}{\tau }}-{\displaystyle \frac{\partial }{\partial x}}\left(d\left(x\right){\displaystyle \frac{\partial u^{m,j}}{\partial x}}\right)+(e\left(x\right)+P\delta \left(x\right))u^{m,j}=g^{m,j}\left(x\right).$$

(14)

Taking into account that $u^{m,0}=u^0\in H^1\left(\mathsf{\Omega }\right)$ and the symmetry of the elliptic operator, we can suppose that $u^{m,j-1}\in H^1\left(\mathsf{\Omega }\right)$ is known.

We will prove that Equation (14) admits a weak solution $u^{m,j}\in H^1\left(\mathsf{\Omega }\right)$ in Step 2. Thus, by induction, we will obtain $u^{m,1},u^{m,2},\dots ,u^{m,m}$ in $H^1\left(\mathsf{\Omega }\right)$, for which we have a uniform estimate.

Then, by employing a standard linear procedure to progress to the new time level and applying estimates to the approximating solutions similar to the approach in [43] (Section 3.2), we establish the proof.

• Step 2. An auxiliary stationary problem with a prototype to Equation (14) is studied.

In more detail, we shall consider the symmetrized self-adjoint problem

$${\left(d\left(x\right)u^{\prime }\right)}^{\prime }+(e\left(x\right)+P\delta \left(x\right))u=g\left(x\right)+Q\delta \left(x\right),x\in \mathsf{\Omega },P=\mathrm{const}.>0,$$

(15)

where ${\mathsf{\Omega }}_1=(-1,0)$, ${\mathsf{\Omega }}_2=(0,1)$ and

$$u(-1)=u_l,u\left(1\right)=u_r.$$

(16)

In $H^1\left(\mathsf{\Omega }\right)$, we introduce a scalar product equivalent to the standard symmetric one, $(u,v)={\int }_{\mathsf{\Omega }}(u^{\prime }{v}^{\prime }+uv)dx$, and define the new scalar product as follows:

$${(u,v)}_{H_{+}^1}\left(\mathsf{\Omega }\right)={\int }_{\mathsf{\Omega }}(d{u}^{\prime }{v}^{\prime }+euv)dx+Pu\left(0\right)v\left(0\right).$$

□

Lemma 1.

Let the conditions in (13) be satisfied. Then, for each $g\in L_2\left(\mathsf{\Omega }\right)$, there exists a unique weak solution $u\in H^1\left(\mathsf{\Omega }\right)$ to the problem (15) and (16), and uniform estimates hold

$${\parallel u\parallel }_{H_{+}^1\left(\mathsf{\Omega }\right)}\le {C(\parallel g\parallel }_{L_2\left(\mathsf{\Omega }\right)}+|u_l|+|u_r\left|\right),$$

(17)

where the constant C is independent of $g,u_l$ and $u_r$. Also, the interface jump conditions are satisfied

$${\left[u\right]}_{x=0}=u(0+)-u(0-)=0,{\left[d{u}^{\prime }\right]}_{x=0}=Pu\left(0\right).$$

(18)

Moreover, let $d,e\in L_{\infty }\left(\mathsf{\Omega }\right)\cap C^1\left({\mathsf{\Omega }}_1\right)\cap C^1\left({\mathsf{\Omega }}_2\right)$, $g\in L_{\infty }\left(\mathsf{\Omega }\right)\cap C\left({\mathsf{\Omega }}_1\right)\cap C\left({\mathsf{\Omega }}_2\right)$. Then $u\in C\left(\mathsf{\Omega }\right)\cap C^2\left({\mathsf{\Omega }}_1\right)\cap C^2\left({\mathsf{\Omega }}_2\right)$, i.e., u is a classical solution of (15) and (16) and the estimate holds:

$${\parallel u\parallel }_{\infty }\le max\left(|u_{-1}|+|u_1|,{\displaystyle \frac{1}{e_0}}{\parallel g\parallel }_{\infty }\right).$$

Proof. 

The generalized solution of (15) and (16) is defined by the identity as follows:

$${(u,w)}_{H_{+}^1\left(\mathsf{\Omega }\right)}\le {(g,w)}_{L_2\left(\mathsf{\Omega }\right)}+\left(uw\right)(-1)+\left(uw\right)\left(1\right)\forall w\in H^1\left(\mathsf{\Omega }\right).$$

For fixed $g\in L_2\left(\mathsf{\Omega }\right)$, ${(g,w)}_{L_2\left(\mathsf{\Omega }\right)}$ is a linear functional, defined on $H^1\left(\mathsf{\Omega }\right),w\in H^1\left(\mathsf{\Omega }\right)$.

Applying the Lax theorem [12,44], we obtain the inequality (17).

The left equality in (18) follows from the trace theorem [44].

With integration by parts of (15) and having in mind that $w\left(x\right)$ is a continuous function, we find

$$[d(0+)u^{\prime }(0+)-d(0-)u^{\prime }(0-)-Pu\left(0\right)]w\left(0\right)=0.$$

Since in general $w\left(0\right)\ne 0$, we obtain the right jump interface conditions in (13). □

4. Exponential Finite Difference Schemes

In this section, we consider the model multilayer parabolic initial boundary value problem with a self-adjoint (symmetric) elliptic part (1)–(4), (9) and (10). Introducing uniform temporal mesh

$${\overline{\omega }}_{\tau }=\{t_n=n\tau ,n=0,1,\dots ,N,\tau =T/N\},$$

using the notation $u_m^n=u_m^n\left(x\right)=u_m(x,t_n)$, and applying the Rothe method to the problem (1)–(4), (9) and (10) yields

$$\begin{array}{ccc}& & {\displaystyle \frac{\partial }{\partial x}}\left(d_m\left(x\right){\displaystyle \frac{\partial u_m^{n+1}}{\partial x^2}}\right)-{\displaystyle \frac{u_m^{n+1}}{\tau }}=-{\displaystyle \frac{u_m^n}{\tau }}-f_m^{n+1}\left(x\right),x\in (l_{m-1},l_m),\hfill \end{array}$$

(19)

$$\begin{array}{ccc}& & {\left[u_m^{n+1}\right]}_{l_m}=u_{m+1}^{n+1}\left(l_m\right)-u_m^{n+1}\left(l_m\right)=0,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& & {\left[{\displaystyle \frac{\partial u_m^{n+1}}{\partial x}}\right]}_{l_m}=r_{m+1}^{n+1}{\displaystyle \frac{\partial u_{m+1}^{n+1}}{\partial x}}\left(l_m\right)-r_m^{n+1}{\displaystyle \frac{\partial u_m^{n+1}}{\partial x}}\left(l_m\right)=u_m^{n+1}\left(l_m\right)+g_m^{n+1},\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & a_0{u}_1^{n+1}\left(0\right)-b_0{\displaystyle \frac{\partial u_1^{n+1}}{\partial x}}\left(0\right)=g_0^{n+1},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}& & a_L{u}_M^{n+1}\left(L\right)+b_L{\displaystyle \frac{\partial u_M^{n+1}}{\partial x}}\left(L\right)=g_L^{n+1},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}& & u_m(x,0)=u_m^0\left(x\right).\hfill \end{array}$$

(24)

We utilize and extend the approximation approach, proposed in [30] for a single layer, in order to obtain the full discretization of the multilayer problem (1)–(4), (9) and (10). Moreover, we develop different second-order approximations for Robin’s type boundary conditions (3) and (4).

Let us consider the nonuniform spatial mesh in the whole interval $[l_0,L]$

$${\overline{\omega }}_h=\{x_i:x_i=x_{i-1}+h_i,i=1,2,\dots ,I,x_0=l_0,x_I=L\}$$

and denote $h=\underset{1\le i{\le }_m}{max}{h}_i$, ${\hslash }_i=0.5(h_i+h_{i+1})$, $i=0,1,\dots ,I-1$. The numerical solution at grid node $(x_i,t_n)$ is denoted by $u_{m,i}^n=u_m^n\left(x_i\right)$, $i=0,1,\dots ,I$.

Suppose that $l_m$, $m=0,1,\dots ,M$ are grid nodes of the spatial mesh ${\overline{\omega }}_h$.

Let

$$\begin{array}{cc}\hfill & u_i^n=u_{m,i}^n,\mathrm{if}{x}_i\in [l_{m-1},l_m],\hfill \\ \hfill & u_{m,i}^n=u_{m+1,i}^n\mathrm{for}{x}_i=l_m,\mathrm{taking}\mathrm{into}\mathrm{account}\left(9\right),\hfill \\ \hfill & u=(u_1,u_2,\dots ,u_M).\hfill \end{array}$$

(25)

We introduce a discontinuous with respect to x function $v(x,t)$, with a finite number of discontinuities in $[l_0,L]$ and present the Steklov-type averaging

$$v_i^{-}\left(t\right)={\displaystyle \frac{1}{h_i}}\underset{x_{i-1}}{\overset{x_i}{\int }}v(x,t)dx,v_i^{+}\left(t\right)={\displaystyle \frac{1}{h_{i+1}}}\underset{x_i}{\overset{x_{i+1}}{\int }}v(x,t)dx.$$

(26)

From (19), in view of (20) and convention (25), we introduce the following problems:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle \frac{d}{dx}}\left(d_i^{-}{\displaystyle \frac{d{u}^{n+1}}{dx}}\right)-{\displaystyle \frac{u^{n+1}}{\tau }}=-{\left({\tilde{f}}_i^{n+1}\right)}^{-},x\in (x_{i-1},x_i)\in (l_{m-1},l_{m+1}),\hfill \\ & u^{n+1}\left(x_{i-1}\right)=u_{i-1},u^{n+1}\left(x_i\right)=u_i^{n+1},\hfill \end{array}\end{array}$$

(27)

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle \frac{d}{dx}}\left(d_i^{+}{\displaystyle \frac{d{u}^{n+1}}{dx}}\right)-{\displaystyle \frac{u^{n+1}}{\tau }}=-{\left({\tilde{f}}_i^{n+1}\right)}^{+},x\in (x_i,x_{i+1})\in (l_{m-1},l_{m+1}).\hfill \\ & u^{n+1}\left(x_i\right)=u_i^{n+1},u^{n+1}\left(x_{i+1}\right)=u_{i+1}^{n+1},\hfill \end{array}\end{array}$$

(28)

where ${d_i}^{\pm }$ is determined as follows:

$$\begin{array}{ccc}& d_i^{-}={d_i}^{+}=d_{m,i},& \mathrm{if}(x_{i-1},x_{i+1})\in (l_{m-1},l_m),\hfill \\ & d_i^{-}=d_{m,i},d_i^{+}=d_{m+1,i},& \mathrm{if}{x}_i=l_m,\hfill \\ & d_i^{-}=d_i^{+}=d_{m+1,i},& \mathrm{if}(x_{i-1},x_{i+1})\in (l_m,l_{m+1}).\hfill \end{array}$$

and similarly for functions ${\displaystyle {\left({\tilde{f}}_i^{n+1}\right)}^{\pm }={\left(\frac{u_i^n}{\tau }+f_i^{n+1}\right)}^{\pm }}$.

The coefficients $d_i^{\pm }$ and the right-hand sides ${\left({\tilde{f}}_i^{n+1}\right)}^{\pm }$ can be derived by the precise solving of the corresponding integrals (26) or by numerical integration.

We start with solving the problem (27), which can be written as a boundary value problem for a second-order ODE with constant coefficients as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}& d_i^{-}{\displaystyle \frac{d^2{u}^{n+1}}{d{x}^2}}-{\displaystyle \frac{u^{n+1}}{\tau }}=-{\left({\tilde{f}}_i^{n+1}\right)}^{-},x\in (x_{i-1},x_i),\hfill \\ & u^{n+1}\left(x_{i-1}\right)=u_{i-1}^{n+1},u^{n+1}\left(x_i\right)=u_i^{n+1}.\hfill \end{array}\end{array}$$

(29)

The solution of the homogeneous ODE in (29) is $u_i^h=C_{i,1}^{-}{e}^{{\lambda }_{i,1}^{-}x}+C_{i,1}^{-}{e}^{{\lambda }_{i,2}^{-}x}$, where $C_{i,1}^{-}$ and $C_{i,1}^{-}$ are constants, and ${\lambda }_{i,1}^{-}$, ${\lambda }_{i,2}^{-}$ are roots of the characteristic equation

$${\lambda }_{i,1}^{-}={\displaystyle \frac{\sqrt{D_i^{-}}}{2d_i^{-}}},{\lambda }_{i,2}^{-}={\displaystyle \frac{-\sqrt{D_i^{-}}}{2d_i^{-}}},{\lambda }_{i,1}^{-}>0,{\lambda }_{i,2}^{-}<0,D_i^{-}={\displaystyle \frac{4}{\tau }}{d}_i^{-}>0.$$

In fact, $C_{i,j}^{-}={\left(C_{i,j}^{n+1}\right)}^{-}$, $j=1,2$, but for simplicity, we omit the superscript $n+1$. Taking into account that ${\left({\tilde{f}}_i^{n+1}\right)}^{-}$ is a constant, we obtain the particular solution of the nonhomogeneous ODE in (29), namely $({\eta }_i^{n+1}{)}^{-}=-{\left({\tilde{f}}_i^{n+1}\right)}^{-}{\left(d_i^{-}{\lambda }_{i,1}^{-}{\lambda }_{i,2}^{-}\right)}^{-1}=\tau {\left({\tilde{f}}_i^{n+1}\right)}^{-}$. Therefore, the general solution is

$$\begin{array}{c}\hfill u^{n+1}\left(x\right)=C_{i,1}^{-}{e}^{{\lambda }_{i,1}^{-}x}+C_{i,1}^{-}{e}^{{\lambda }_{i,2}^{-}x}+\tau {\left({\tilde{f}}_i^{n+1}\right)}^{-},x\in [x_{i-1},x_i].\end{array}$$

(30)

The constants $C_{i,1}^{-}$ and $C_{i,2}^{-}$ are determined using BCs in (29) as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & C_{i,1}^{-}={\displaystyle \frac{(u_{i-1}^{n+1}-{\eta }_i^{-})e^{{\lambda }_{i,2}^{-}{x}_i}-(u_i^{n+1}-{\eta }_i^{-})e^{{\lambda }_{i,2}^{-}{x}_{i-1}}}{e^{{\lambda }_{i,1}^{-}{x}_{i-1}+{\lambda }_{i,2}^{-}{x}_i}-e^{{\lambda }_{i,1}^{-}{x}_i+{\lambda }_{i,2}^{-}{x}_{i-1}}}},\hfill \\ & C_{i,2}^{-}={\displaystyle \frac{(u_i^{n+1}-{\eta }_i^{-})e^{{\lambda }_{i,1}^{-}{x}_{i-1}}-(u_{i-1}^{n+1}-{\eta }_i^{-})e^{{\lambda }_{i,1}^{-}{x}_i}}{e^{{\lambda }_{i,1}^{-}{x}_{i-1}+{\lambda }_{i,2}^{-}{x}_i}-e^{{\lambda }_{i,1}^{-}{x}_i+{\lambda }_{i,2}^{-}{x}_{i-1}}}}.\hfill \end{array}\end{array}$$

(31)

Next, by substituting (31) in (30), differentiating $u^{n+1}\left(x\right)$, taking the limit $x\to x_i-0$, and applying some evaluation and reductions, we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{d{u}^{n+1}}{dx}}{|}_{x=x_i-0}=& -{\displaystyle \frac{({\lambda }_{i,1}^{-}-{\lambda }_{i,2}^{-})e^{{\lambda }_{i,2}^{-}{h}_i}}{1-e^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}}{u}_{i-1}^{n+1}+{\displaystyle \frac{{\lambda }_{i,1}^{-}-{\lambda }_{i,2}^{-}{e}^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}{1-e^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}}{u}_i^{n+1}\hfill \\ & -{\left({\eta }_i^{n+1}\right)}^{-}{\displaystyle \frac{{\lambda }_{i,1}^{-}(1-e^{{\lambda }_{i,2}^{-}{h}_i})+{\lambda }_{i,2}^{-}{e}^{{\lambda }_{i,2}^{-}{h}_i}(1-e^{-{\lambda }_{i,1}^{-}{h}_i})}{1-e^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}},i=1,2,\dots ,I-1.\hfill \end{array}\end{array}$$

(32)

In the same way, we solve the problem (28), represented in an equivalent form as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}& d_i^{+}{\displaystyle \frac{d^2{u}^{n+1}}{d{x}^2}}-{\displaystyle \frac{u^{n+1}}{\tau }}=-{\left({\tilde{f}}_i^{n+1}\right)}^{+},x\in (x_i,x_{i+1}),\hfill \\ & u^{n+1}\left(x_i\right)=u_i^{n+1},u^{n+1}\left(x_{i+1}\right)=u_{i+1}^{n+1}.\hfill \end{array}\end{array}$$

The corresponding solution is

$$u^{n+1}\left(x\right)=C_{i,1}^{+}{e}^{{\lambda }_{i,1}^{+}x}+C_{i,1}^{+}{e}^{{\lambda }_{i,2}^{+}x}+{\eta }_i^{+},{\eta }_i^{+}=\tau {\left({\tilde{f}}_i^{n+1}\right)}^{+},x\in [x_i,x_{i+1}],$$

(33)

where

$$\begin{array}{ccc}& & C_{i,1}^{+}={\displaystyle \frac{(u_i-{\eta }_i^{+})e^{{\lambda }_{i,2}^{-}{x}_{i+1}}-(u_{i+1}-{\eta }_i^{+})e^{{\lambda }_{i,2}^{+}{x}_i}}{e^{{\lambda }_{i,1}^{+}{x}_i+{\lambda }_{i,2}^{+}{x}_{i+1}}-e^{{\lambda }_{i,1}^{+}{x}_{i+1}+{\lambda }_{i,2}^{+}{x}_i}}},C_{i,2}^{+}={\displaystyle \frac{(u_{i+1}-{\eta }_i^{+})e^{{\lambda }_{i,1}^{+}{x}_i}-(u_i-{\eta }_i^{-})e^{{\lambda }_{i,1}^{+}{x}_{i+1}}}{e^{{\lambda }_{i,1}^{+}{x}_i+{\lambda }_{i,2}^{+}{x}_{i+1}}-e^{{\lambda }_{i,1}^{+}{x}_{i+1}+{\lambda }_{i,2}^{+}{x}_i}}},\hfill \\ & & {\lambda }_{i,1}^{+}={\displaystyle \frac{\sqrt{D_i^{+}}}{2{\left(d_i^{n+1}\right)}^{+}}},{\lambda }_{i,2}^{+}={\displaystyle \frac{-\sqrt{D_i^{-}}}{2{\left(d_i^{n+1}\right)}^{+}}},{\lambda }_{i,1}^{+}>0,{\lambda }_{i,2}^{+}<0,D_i^{+}={\displaystyle \frac{4}{\tau }}{d}_i^{+}>0.\hfill \end{array}$$

Differentiating $u^{n+1}\left(x\right)$ in (33) and taking the limit $x\to x_i+0$, we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{d{u}^{n+1}}{dx}}{|}_{x=x_i+0}=& {\displaystyle \frac{({\lambda }_{i,1}^{+}-{\lambda }_{i,2}^{+})e^{-{\lambda }_{i,1}^{+}{h}_{i+1}}}{1-e^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}}}{u}_{i+1}^{n+1}-{\displaystyle \frac{{\lambda }_{i,1}^{+}{e}^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}-{\lambda }_{i,2}^{+}}{1-e^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}}}{u}_i^{n+1}\hfill \\ & -{\left({\eta }_i^{n+1}\right)}^{+}{\displaystyle \frac{{\lambda }_{i,1}^{+}{e}^{-{\lambda }_{i,1}^{+}{h}_{i+1}}(1-e^{{\lambda }_{i,2}^{+}{h}_{i+1}})+{\lambda }_{i,2}^{+}(1-e^{-{\lambda }_{i,1}^{+}{h}_{i+1}})}{1-e^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}}}.\hfill \end{array}\end{array}$$

(34)

At regular points $x_i\ne l_m$, $m=1,2,\dots ,M-1$, we have ${\displaystyle d\left(x\right)\frac{d{u}^{n+1}\left(x\right)}{dx}{|}_{x=x_i-0}=d\left(x\right)\frac{d{u}^{n+1}\left(x\right)}{dx}{|}_{x=x_i+0}}$, while at interface nodes $x_i=l_m$, $m=1,2,\dots ,M-1$ the relation (10) holds. Therefore, from (21), (32) and (34), we obtain a symmetric system of algebraic equations

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & -{\chi }_{m+1}^{n+1}{\displaystyle \frac{({\lambda }_{i,1}^{+}-{\lambda }_{i,2}^{+})e^{-{\lambda }_{i,1}^{+}{h}_{i+1}}}{1-e^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}}}{u}_{i+1}^{n+1}\hfill \\ & +\left({\chi }_{m+1}^{n+1}{\displaystyle \frac{{\lambda }_{i,1}^{+}{e}^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}-{\lambda }_{i,2}^{+}}{1-e^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}}}+{\chi }_m^{n+1}{\displaystyle \frac{{\lambda }_{i,1}^{-}-{\lambda }_{i,2}^{-}{e}^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}{1-e^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}}+{\tilde{\chi }}_m\right)u_i^{n+1}\hfill \\ & -{\chi }_m^{n+1}{\displaystyle \frac{({\lambda }_{i,1}^{-}-{\lambda }_{i,2}^{-})e^{{\lambda }_{i,2}^{-}{h}_i}}{1-e^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}}{u}_{i-1}^{n+1}={\overline{F}}_i^{n+1},\hfill \\ & {\overline{F}}_i^{n+1}=-{\chi }_{m+1}^{n+1}{\left({\eta }_i^{n+1}\right)}^{+}{\displaystyle \frac{{\lambda }_{i,1}^{+}{e}^{-{\lambda }_{i,1}^{+}{h}_{i+1}}(1-e^{{\lambda }_{i,2}^{+}{h}_{i+1}})+{\lambda }_{i,2}^{+}(1-e^{-{\lambda }_{i,1}^{+}{h}_{i+1}})}{1-e^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}}}\hfill \\ & +{\chi }_m^{n+1}{\left({\eta }_i^{n+1}\right)}^{-}{\displaystyle \frac{{\lambda }_{i,1}^{-}(1-e^{{\lambda }_{i,2}^{-}{h}_i})+{\lambda }_{i,2}^{-}{e}^{{\lambda }_{i,2}^{-}{h}_i}(1-e^{-{\lambda }_{i,1}^{-}{h}_i})}{1-e^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}}-{\tilde{\chi }}_m{g}_m^{n+1},i=1,\dots ,I-1,\hfill \end{array}\end{array}$$

(35)

where

$${\chi }_m^{n+1}=\left\{\begin{array}{cc}{r}_m^{n+1},\hfill & x_i=l_m,m=1,2.\dots ,M,\hfill \\ d_m,\hfill & x_i\ne l_m,m=1,2.\dots ,M\hfill \end{array}\right.$$

and

$${\tilde{\chi }}_m=\left\{\begin{array}{cc}1,\hfill & x_i=l_m,m=1,2.\dots ,M,\hfill \\ 0,\hfill & x_i\ne l_m,m=1,2,\dots ,M.\hfill \end{array}\right.$$

Next, we approximate the boundary conditions. For $b_0=0$ and $b_L=0$, we obtain Dirichlet BCs, and the discretization is trivial. Let $b_0\ne 0$ and $b_L\ne 0$. To derive the discretization for (22), we proceed as follows: At each time layer, we consider the following boundary value ODE problem:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle \frac{\partial u^{n+1}}{\partial x}}\left(x\right)-{\displaystyle \frac{a_0}{b_0}}{u}^{n+1}\left(x\right)=-{\displaystyle \frac{g_0^{n+1}}{b_0}},x\in (0,h),\hfill \\ & u^{n+1}\left(x_0\right)=u_0^{n+1},u^{n+1}\left(x_1\right)=u_1^{n+1}.\hfill \end{array}\end{array}$$

(36)

The general solution is

$${\displaystyle u^{n+1}\left(x\right)=C{e}^{\frac{a_0}{b_0}x}+\frac{g_0^{n+1}}{a_0}.}$$

Using boundary conditions in (36), we obtain

$$e^{{\displaystyle \frac{a_0}{b_0}}h}{u}_0^{n+1}-u_1^{n+1}={\displaystyle \frac{g_0^{n+1}}{a_0}}\left(e^{{\displaystyle \frac{a_0}{b_0}}h}-1\right).$$

(37)

We discretize the right boundary condition (23) in the same manner. Consider the following boundary value ODE problem:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle \frac{\partial u^{n+1}}{\partial x}}\left(x\right)+{\displaystyle \frac{a_L}{b_L}}{u}^{n+1}\left(x\right)={\displaystyle \frac{g_L^{n+1}}{b_L}},x\in (L-h,L),\hfill \\ & u^{n+1}\left(x_{I-1}\right)=u_{I-1}^{n+1},u^{n+1}\left(x_I\right)=u_I^{n+1}.\hfill \end{array}\end{array}$$

(38)

The general solution is

$${\displaystyle u^{n+1}\left(x\right)=C{e}^{-\frac{a_L}{b_L}x}+\frac{g_L^{n+1}}{a_L}.}$$

Using boundary conditions in (38), we obtain

$$e^{{\displaystyle \frac{a_L}{b_L}}h}{u}_I^{n+1}-u_{I-1}^{n+1}={\displaystyle \frac{g_L^{n+1}}{a_L}}\left(e^{{\displaystyle \frac{a_L}{b_L}}h}-1\right).$$

(39)

The scheme is completed by the following initial condition:

$$u(x,0)=u^0\left(x\right).$$

(40)

Let us rewrite the discretization (35), (37), (39) and (40) in a more compact symmetric form:

$$\begin{array}{ccc}& & D_0^{n+1}{u}_0^{n+1}+B_0^{n+1}{u}_1^{n+1}=F_0^{n+1},\hfill \\ & & A_i^{n+1}{u}_{i-1}^{n+1}+D_i^{n+1}{u}_i^{n+1}+B_i^{n+1}{u}_{i+1}^{n+1}={\overline{F}}_i^{n+1},i=1,2,\dots ,I-1,\hfill \\ & & A_I^{n+1}{u}_{I-1}^{n+1}+D_I^{n+1}{u}_I^{n+1}=F_I^{n+1},\hfill \\ & & u_i^0=u^0\left(x_i\right),i=0,1,\dots ,I.\hfill \end{array}$$

(41)

where

$$\begin{array}{ccc}& & {\displaystyle D_0^{n+1}=e^{\frac{a_0}{b_0}h},B_0^{n+1}=-1,F_0^{n+1}=\frac{g_0^{n+1}}{a_0}\left(e^{\frac{a_0}{b_0}h}-1\right),}\hfill \\ & & A_i^{n+1}=-{\chi }_m^{n+1}{\displaystyle \frac{({\lambda }_{i,1}^{-}-{\lambda }_{i,2}^{-})e^{{\lambda }_{i,2}^{-}{h}_i}}{1-e^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}},\hfill \\ & & D_i^{n+1}={\chi }_{m+1}^{n+1}{\displaystyle \frac{{\lambda }_{i,1}^{+}{e}^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}-{\lambda }_{i,2}^{+}}{1-e^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}}}+{\chi }_m^{n+1}{\displaystyle \frac{{\lambda }_{i,1}^{-}-{\lambda }_{i,2}^{-}{e}^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}{1-e^{({\lambda }_{i,2}^{-}-{\lambda }_{i,1}^{-})h_i}}}+{\tilde{\chi }}_m,\hfill \\ & & B_i^{n+1}=-\left({\chi }_{m+1}^{n+1}\right){\displaystyle \frac{({\lambda }_{i,1}^{+}-{\lambda }_{i,2}^{+})e^{-{\lambda }_{i,1}^{+}{h}_{i+1}}}{1-e^{({\lambda }_{i,2}^{+}-{\lambda }_{i,1}^{+})h_{i+1}}}},\hfill \\ & & A_I^{n+1}=-1,D_I^{n+1}=e^{{\displaystyle \frac{a_L}{b_L}}h},F_I^{n+1}={\displaystyle \frac{g_L^{n+1}}{a_L}}\left(e^{{\displaystyle \frac{a_L}{b_L}}h}-1\right).\hfill \end{array}$$

Similar to [30], from the discrete maximum principle [45], it follows that, for $u(x,0)\ge 0$, $g_0\left(x\right)\ge 0$, $g_L\left(x\right)\ge 0$, and $g_m\left(x\right)\le 0$, the numerical scheme (41) is positivity preserving, i.e., $u_i^n\ge 0$, $i=0,1,\dots ,I$, $n=1,2,\dots ,N$. Moreover, the positivity of the solution is preserved, and the scheme remains second-order accurate, even when there is a convection term in (1).

5. Inverse Problems with Unknown External Boundary Conditions

In this section, we formulate and solve numerically inverse problems for recovering functions $g_0$ and/or $g_0$ for given measurements in different layers.

5.1. Formulation of the Inverse Problems

The inverse problems are formulated considering the functions $g_0\left(t\right)$ and/or $g_L\left(t\right)$ to be unknown and must be identified using some overspecified data.

• Depending on the observed data, the inverse problems IP for the determination of the boundary functions $g_0\left(t\right)$ and/or $g_L\left(t\right)$ in (3) and (4) can be formulated as follows:

• IP1: For given observation $u_1(x_1^{*},t)={\psi }_1\left(t\right)$, $x_1^{*}\in {\mathsf{\Omega }}_1$ or $u_M(x_M^{*},t)={\psi }_M\left(t\right),x_M^{*}\in {\mathsf{\Omega }}_M$, find $\left(u_1,\dots ,u_M,g_0\left(t\right)\right)$ or $\left(u_1,\dots ,u_M,g_L\left(t\right)\right)$, respectively.
• IP2: For given observation $u_s(x_s^{*},t)={\psi }_s\left(t\right)$, $x_s^{*}\in {\mathsf{\Omega }}_s$, $s\ne 1$ and $s\ne M$, find $\left(u_1\dots ,u_M,g_L\left(t\right)\right)$ or $\left(u_1\dots ,u_M,g_0\left(t\right)\right)$.
• IP3: For given observations $u_1(x_1^{*},t)={\psi }_1\left(t\right)$, $x_1^{*}\in {\mathsf{\Omega }}_1$ and $u_M(x_M^{*},t)={\psi }_M\left(t\right)$, $x_M^{*}\in {\mathsf{\Omega }}_M$, find $\left(u_1,\dots ,u_M,g_0\left(t\right),g_L\left(t\right)\right)$.
• IP4: For given observations $u_1(x_1^{*},t)={\psi }_1\left(t\right)$, $x_1^{*}\in {\mathsf{\Omega }}_1$ and $u_i(x_s^{*},t)={\psi }_s\left(t\right)$, $x_s^{*}\in {\mathsf{\Omega }}_s$, $s\ne M$, find $\left(u_1,\dots ,u_M,g_0\left(t\right),g_L\left(t\right)\right)$.
• IP5: For given observations $u_s(x_s^{*},t)={\psi }_s\left(t\right)$, $x_s^{*}\in {\mathsf{\Omega }}_s$ and $u_M(x_s^{*},t)={\psi }_M\left(t\right)$, $x_M^{*}\in {\mathsf{\Omega }}_M$, $snot=1$, find $\left(u_1,\dots ,u_M,g_0\left(t\right),g_L\left(t\right)\right)$.
• IP6: For given observations $u_{s_1}(x_{s_1}^{*},t)={\psi }_{s_1}\left(t\right)$, $x_{s_1}^{*}\in {\mathsf{\Omega }}_{s_1}$ and $u_{s_2}(x_{s_2}^{*},t)={\psi }_{s_2}\left(t\right)$, $x_{s_2}^{*}\in {\mathsf{\Omega }}_{s_2}$, $s_1,s_2\ne \{1,M\}$, find $\left(u_1,\dots ,u_M,g_0\left(t\right),g_L\left(t\right)\right)$.

5.2. Numerical Solution

In this section, we discuss numerically solving IP1–IP6. Because of the unified discretization of the direct problem, which is valid at each space grid node, independent of the region ${\mathsf{\Omega }}_m$, we can classify these IPs into two groups.

• IPa: Find $\left(u_1,\dots ,u_M,g_0\left(t\right)\right)$ or $\left(u_1,\dots ,u_M,g_L\left(t\right)\right)$, if measurements at one space point $x_s^{*}\in {\overline{\omega }}_h$, $x_s^{*}\ne \{x_0,l_m,x_L\}$, $m=1,2,\dots ,M$ is given, i.e., $u(x_s^{*},t)={\psi }_s\left(t\right)$, $x_s^{*}\in {\mathsf{\Omega }}_s$.
• IPb: Find $\left(u_1,\dots ,u_M,g_0\left(t\right),g_L\left(t\right)\right)$, if measurements at two space points $x_{s_1}^{*}\in {\overline{\omega }}_h$, $x_{s_2}^{*}\in {\overline{\omega }}_h$, $x_{s_1}^{*},x_{s_2}^{*}\ne \{x_0,l_m,x_L\}$, $m=1,2,\dots ,M$ are given, i.e., $u(x_{s_1}^{*},t)={\psi }_{s_1}\left(t\right)$, $x_{s_1}^{*}\in {\mathsf{\Omega }}_{s_1}$ and $u(x_{s_2}^{*},t)={\psi }_{s_2}\left(t\right)$, $x_{s_2}^{*}\in {\mathsf{\Omega }}_{s_2}$.

Since we deal with nonuniform spatial mesh, the measurement can always be chosen to be grid nodes.

Consider IPa. In order to find $\left(u_1,\dots ,u_M,g_0\left(t\right)\right)$, we introduce the symmetric with respect to $g_0\left(t\right)$ decomposition

$$u_i^{n+1}=U_i^{n+1}+g_0\left(t\right)V_i^{n+1}.$$

(42)

Substituting in (41), we obtain the following problems:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & D_0^{n+1}{U}_0^{n+1}+B_0^{n+1}{U}_1^{n+1}={\tilde{F}}_0^{n+1},\hfill \\ & A_i^{n+1}{U}_{i-1}^{n+1}+D_i^{n+1}{U}_i^{n+1}+B_i^{n+1}{U}_{i+1}^{n+1}={\overline{F}}_i^{n+1},i-1,2,\dots ,I-1,\hfill \\ & A_I^{n+1}{U}_{I-1}^{n+1}+D_I^{n+1}{U}_I^{n+1}=F_I^{n+1},\hfill \end{array}\end{array}$$

(43)

where

$${\tilde{F}}_0^{n+1}={\displaystyle \frac{1}{a_0}}\left(e^{{\displaystyle \frac{a_0}{b_0}}h}-1\right),$$

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & D_0^{n+1}{V}_0^{n+1}+B_0^{n+1}{V}_1^{n+1}=1,\hfill \\ & A_i^{n+1}{V}_{i-1}^{n+1}+D_i^{n+1}{V}_i^{n+1}+B_i^{n+1}{V}_{i+1}^{n+1}=0,i-1,2,\dots ,I-1,\hfill \\ & A_I^{n+1}{V}_{I-1}^{n+1}+D_I^{n+1}{V}_I^{n+1}=0.\hfill \end{array}\end{array}$$

(44)

After solving the problems in (43) and (44), we use the measurement and (42) to obtain

$$g_0^{n+1}={\displaystyle \frac{{\psi }_s^{n+1}-U^{n+1}\left(x_s^{*}\right)}{V{\left(x_s^{*}\right)}^{n+1}}},n=0,1,\dots ,N-1.$$

(45)

In the same manner, we find $\left(u_1,\dots ,u_M,g_L\left(t\right)\right)$. Consider the symmetric form with respect to $g_L\left(t\right)$; decomposition is

$$u_i^{n+1}=U_i^{n+1}+g_L\left(t\right)V_i^{n+1}.$$

(46)

Substituting in (41), we obtain the following problems:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & D_0^{n+1}{U}_0^{n+1}+B_0^{n+1}{U}_1^{n+1}=F_0^{n+1},\hfill \\ & A_i^{n+1}{U}_{i-1}^{n+1}+D_i^{n+1}{U}_i^{n+1}+B_i^{n+1}{U}_{i+1}^{n+1}={\overline{F}}_i^{n+1},i-1,2,\dots ,I-1,\hfill \\ & A_I^{n+1}{U}_{I-1}^{n+1}+D_I^{n+1}{U}_I^{n+1}={\tilde{F}}_I^{n+1},\hfill \end{array}\end{array}$$

(47)

where

$${\tilde{F}}_I^{n+1}={\displaystyle \frac{1}{a_L}}\left(e^{{\displaystyle \frac{a_L}{b_L}}h}-1\right),$$

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & D_0^{n+1}{V}_0^{n+1}+B_0^{n+1}{V}_1^{n+1}=0,\hfill \\ & A_i^{n+1}{V}_{i-1}^{n+1}+D_i^{n+1}{V}_i^{n+1}+B_i^{n+1}{V}_{i+1}^{n+1}=0,i-1,2,\dots ,I-1,\hfill \\ & A_I^{n+1}{V}_{I-1}^{n+1}+D_I^{n+1}{V}_I^{n+1}=1.\hfill \end{array}\end{array}$$

(48)

Using the solutions of (47) and (48), the measurement, and (46), we find

$$g_L^{n+1}={\displaystyle \frac{{\psi }_s^{n+1}-U^{n+1}\left(x_s^{*}\right)}{V{\left(x_s^{*}\right)}^{n+1}}},n=0,1,\dots ,N-1.$$

Further, we will focus on IPb. In order to find $\left(u_1,\dots ,u_M,g_0\left(t\right),g_L\left(t\right)\right)$, we consider the decomposition

$$u_i^{n+1}=U_i^{n+1}+g_0\left(t\right)V_i^{n+1}+g_L\left(t\right)W_i^{n+1}.$$

(49)

Substituting in (41), we obtain three problems as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & D_0^{n+1}{U}_0^{n+1}+B_0^{n+1}{U}_1^{n+1}={\tilde{F}}_0^{n+1},\hfill \\ & A_i^{n+1}{U}_{i-1}^{n+1}+D_i^{n+1}{U}_i^{n+1}+B_i^{n+1}{U}_{i+1}^{n+1}={\overline{F}}_i^{n+1},i-1,2,\dots ,I-1,\hfill \\ & A_I^{n+1}{U}_{I-1}^{n+1}+D_I^{n+1}{U}_I^{n+1}={\tilde{F}}_I^{n+1},\hfill \end{array}\end{array}$$

(50)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & D_0^{n+1}{V}_0^{n+1}+B_0^{n+1}{V}_1^{n+1}=1,\hfill \\ & A_i^{n+1}{V}_{i-1}^{n+1}+D_i^{n+1}{V}_i^{n+1}+B_i^{n+1}{V}_{i+1}^{n+1}=0,i-1,2,\dots ,I-1,\hfill \\ & A_I^{n+1}{V}_{I-1}^{n+1}+D_I^{n+1}{V}_I^{n+1}=0.\hfill \end{array}\end{array}$$

(51)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & D_0^{n+1}{W}_0^{n+1}+B_0^{n+1}{W}_1^{n+1}=0,\hfill \\ & A_i^{n+1}{W}_{i-1}^{n+1}+D_i^{n+1}{W}_i^{n+1}+B_i^{n+1}{W}_{i+1}^{n+1}=0,i-1,2,\dots ,I-1,\hfill \\ & A_I^{n+1}{W}_{I-1}^{n+1}+D_I^{n+1}{W}_I^{n+1}=1.\hfill \end{array}\end{array}$$

(52)

Once we find the solutions of (50)–(52), we use the measurements and decomposition (49) to determine $g_0^{n+1}$ and $g_L^{n+1}$, $n=0,1,\dots ,N-1$, solving the system as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\psi }_{s_1}^{n+1}=U^{n+1}\left(x_{s_1}^{*}\right)+g_0^{n+1}{V}^{n+1}\left(x_{s_1}^{*}\right)+g_L\left(t\right)W^{n+1}\left(x_{s_1}^{*}\right),\hfill \\ & {\psi }_{s_2}^{n+1}=U^{n+1}\left(x_{s_1}^{*}\right)+g_0^{n+1}{V}^{n+1}\left(x_{s_2}^{*}\right)+g_L\left(t\right)W^{n+1}\left(x_{s_2}^{*}\right).\hfill \end{array}\end{array}$$

Therefore, we derive

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & g_0^{n+1}={\displaystyle \frac{W^{n+1}\left(x_{s_2}^{*}\right)\left({\psi }_{s_1}^{n+1}-U{\left(x_{s_1}^{*}\right)}^{n+1}\right)-W^{n+1}\left(x_{s_1}^{*}\right)\left({\psi }_{s_2}^{n+1}-U_{x_{s_2}^{*}}^{n+1}\right)}{V^{n+1}\left(x_{s_1}^{*}\right)W^{n+1}\left(x_{s_2}^{*}\right)-V^{n+1}\left(x_{s_2}^{*}\right)W^{n+1}\left(x_{s_1}^{*}\right)}},\hfill \\ & g_L^{n+1}={\displaystyle \frac{V^{n+1}\left(x_{s_1}^{*}\right)\left({\psi }_{s_2}^{n+1}-U{\left(x_{s_2}^{*}\right)}^{n+1}\right)-V^{n+1}\left(x_{s_2}^{*}\right)\left({\psi }_{s_1}^{n+1}-U{\left(x_{s_1}^{*}\right)}^{n+1}\right)}{V^{n+1}\left(x_{s_1}^{*}\right)W^{n+1}\left(x_{s_2}^{*}\right)-V^{n+1}\left(x_{s_2}^{*}\right)W^{n+1}\left(x_{s_1}^{*}\right)}}.\hfill \end{array}\end{array}$$

(53)

6. Numerical Results

In this section, we illustrate the accuracy of the proposed numerical method for solving direct and inverse problems IPa and IPb. Additionally, we show the influence of the location of the point(s) (in which region(s) ${\mathsf{\Omega }}_1$, ${\mathsf{\Omega }}_2$, or ${\mathsf{\Omega }}_3$ belong) of measurement(s) on the precision of the solution, i.e., we study IP1–IP6. We take a test example with Robin’s BC at $x=l_0$ and Dirichlet BC at $x=L$, namely $a_0=b_0=a_L=1$, $b_L=0$. We let $M=3$ and $l_0=0$, $l_1=1$, $l_2=2$, $L=l_3=3$. All computations are performed on uniform spatial mesh with step size h.

Example 1

(Direct problem: convergence test). First, we examine the order of convergence of the difference scheme (41) for solving the direct problem. We consider a test example with the following coefficients, right-hand sides, and initial conditions:

$$\begin{array}{ccc}{\displaystyle d_1\left(x\right)=3x+1,}\hfill & f_1(x,t)=-e^{-t}-38x^3-12x^2,\hfill & x\in [l_0,l_1],\hfill \\ {\displaystyle d_2\left(x\right)=5,}\hfill & {\displaystyle f_2(x,t)=-e^{-t}-\frac{{\pi }^2}{16}cos\frac{\pi x}{2},}\hfill & x\in [l_1,l_2],\hfill \\ {\displaystyle d_3\left(x\right)=x^2+2,}\hfill & {\displaystyle f_3(x,t)=-e^{-t}-\frac{x}{2},}\hfill & x\in [l_2,l_3],\hfill \\ g_0\left(t\right)=e^{-t},\hfill & {\displaystyle g_L\left(t\right)=\frac{3}{4}+1+e^{-t},}\hfill \\ {\displaystyle r_1=\frac{1}{4},r_2=\frac{4t}{3},}\hfill & g_1\left(t\right)=-e^{-t}+t-2,\hfill \\ {\displaystyle r_3=\frac{7t}{12},}\hfill & {\displaystyle g_2\left(t\right)=\frac{\pi t}{6}-\frac{5}{2}-e^{-t},}\hfill \end{array}$$

$$u^0\left(x\right)=\left\{\begin{array}{cc}{\displaystyle u_1^0\left(x\right)=x^4+1,}\hfill & x\in [l_0,l_1],\hfill \\ {\displaystyle u_2^0\left(x\right)=\frac{1}{4}sin\frac{\pi x}{2}+\frac{3x}{4}+1,}\hfill & x\in [l_1,l_2],\hfill \\ {\displaystyle u_3^0\left(x\right)=\frac{x}{4}+2,}\hfill & x\in [l_2,l_3].\hfill \end{array}\right.$$

Therefore, the solution of the problem (1)–(4), (9) and (10) is

$$u(x,t)=\left\{\begin{array}{cc}{\displaystyle u_1(x,t)=x^4+e^{-t},}\hfill & x\in [l_0,l_1],\hfill \\ {\displaystyle u_2(x,t)=\frac{1}{4}sin\frac{\pi x}{2}+\frac{3x}{4}+e^{-t},}\hfill & x\in [l_1,l_2],\hfill \\ {\displaystyle u_3(x,t)=\frac{x}{4}+1+e^{-t},}\hfill & x\in [l_2,l_3].\hfill \end{array}\right.$$

Denote by ${\tilde{u}}_i^n$ the numerical solution, obtained by (41) at time layer $t_n$ and spatial grid node $x_i$ and let $E_i^n={\tilde{u}}_i^n-u(x_i,t_n)$. The computational results are estimated by [15,46]

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle {\mathcal{E}}_{\infty }={\mathcal{E}}_{\infty }^I=\underset{1\le n\le N}{max}\underset{0\le i\le I}{max}\left|E_i^n\right|,{\mathcal{E}}_2={\mathcal{E}}_2^I=\sqrt{h\tau \sum _{i=0}^I\sum _{n=1}^N{\left(E_i^n\right)}^2},}\hfill \\ & {\displaystyle C{R}_{\infty }={log}_2\frac{{\mathcal{E}}_{\infty }^I}{{\mathcal{E}}_{\infty }^{2I}},C{R}_2={log}_2\frac{{\mathcal{E}}_2^I}{{\mathcal{E}}_2^{2I}}.}\hfill \end{array}\end{array}$$

(54)

In

Table 1. Errors and convergence rates of the numerical solution of the direct problem, $\tau \approx h^2$; Example 1.

I	${\mathcal{E}}_{\mathsf{\infty }}$	${\mathit{CR}}_{\mathsf{\infty }}$	${\mathcal{E}}_2$	${\mathit{CR}}_2$
120	2.1992$\times {10}^{-2}$		1.0855$\times {10}^{-2}$
240	1.5012$\times {10}^{-2}$	0.5509	7.3008$\times {10}^{-3}$	0.5723
480	7.2866$\times {10}^{-3}$	1.0428	3.5299$\times {10}^{-3}$	1.0484
960	2.5956$\times {10}^{-3}$	1.4892	1.2671$\times {10}^{-3}$	1.4781
1920	6.7216$\times {10}^{-4}$	1.9492	3.5650$\times {10}^{-4}$	1.8295

, we present the errors and convergence rates of the numerical solution of the direct problem, computed for fixed ratio $\tau \approx h^2$, while the results in

Table 2. Errors and convergence rates of the numerical solution of the direct problem, $\tau =h$; Example 1.

I	${\mathcal{E}}_{\mathsf{\infty }}$	${\mathit{CR}}_{\mathsf{\infty }}$	${\mathcal{E}}_2$	${\mathit{CR}}_2$
120	7.7520$\times {10}^{-3}$		3.9829$\times {10}^{-3}$
240	3.7614$\times {10}^{-3}$	1.0433	1.8154$\times {10}^{-3}$	1.1335
480	1.8516$\times {10}^{-3}$	1.0225	8.6736$\times {10}^{-4}$	1.0656
960	9.1846$\times {10}^{-4}$	1.0115	4.2413$\times {10}^{-4}$	1.0321
1920	4.5739$\times {10}^{-4}$	1.0058	2.0975$\times {10}^{-4}$	1.0158

are obtained for $\tau =h$. The final time is $T=0.5$. We observe that the order of convergence in space is second, while in time, it is first. Therefore, the accuracy of the discrete scheme (41) is $O(\tau +h^2)$.

For comparison, only in this test, we consider layers with different lengths, i.e., $l_0=0$, $l_1=1.5$, $l_2=2$, $l_3=3$. For the computations, we use the same coefficients $d_m$ and $r_m$ as before. The initial functions; right-hand sides; functions $g_0$, $g_2$, and $g_L$; and solutions in $[l_1,l_2]$ and $[l_2,l_3]$ are the same as in the previous example. The only modification is the initial solution at the first layer, which reflects the solution, the function $g_1$, and the corresponding right-hand side $f_1$. Now, we take $u^0\left(x\right)={\displaystyle \frac{2}{9}}{x}^4+1+{\displaystyle \frac{\sqrt{2}}{8}}$, $u(x,t)=e^{-t}+{\displaystyle \frac{2}{9}}{x}^4+{\displaystyle \frac{\sqrt{2}}{8}}$, $x\in [l_0,l_1]$. The results from computations for $\tau \approx h^2$, $T=0.5$ are given in

Table 3. Errors and convergence rates of the numerical solution of the direct problem, $\tau \approx h^2$; Example 1.

I	${\mathcal{E}}_{\mathsf{\infty }}$	${\mathit{CR}}_{\mathsf{\infty }}$	${\mathcal{E}}_2$	${\mathit{CR}}_2$
120	1.8200$\times {10}^{-2}$		1.0099$\times {10}^{-2}$
240	1.1761$\times {10}^{-2}$	0.6299	6.4996$\times {10}^{-3}$	0.6358
480	5.6813$\times {10}^{-3}$	1.0497	3.1291$\times {10}^{-3}$	1.0546
960	2.0976$\times {10}^{-3}$	1.4375	1.1520$\times {10}^{-3}$	1.4416
1920	6.0960$\times {10}^{-4}$	1.7828	3.3755$\times {10}^{-4}$	1.7710

. We observe that the results in Table 1 and Table 3 are close, and the spatial order of convergence is one and the same. The different lengths of the layers (position of the interface points) do not significantly influence the accuracy.

Note that the discretization (41) is derived semi-analytically, and in the stationary case for piecewise constant diffusion, it is exact; see [30]. Therefore, if each diffusion coefficient, $d_m$, $m=1,2,\dots ,M$ is a constant, potentially different for each layer, and the numerical scheme (41) is very precise. Moreover, it is exact and meshless in space and only involves error from the time discretization. We will compare the precision of the numerical scheme (41) with the second-order problem in space, i.e., the implicit finite difference scheme proposed in [6], in the case of spatially dependent coefficients $d_m$. The results from computations with the method proposed in [6] for $T=0.5$ and $\tau =h$ are presented in

Table 4. Errors and convergence rates of the numerical solution of the direct problem, [6], $\tau =h$; Example 1.

I	${\mathcal{E}}_{\mathsf{\infty }}$	${\mathit{CR}}_{\mathsf{\infty }}$	${\mathcal{E}}_2$	${\mathit{CR}}_2$
120	8.0572$\times {10}^{-3}$		3.9049$\times {10}^{-3}$
240	3.9901$\times {10}^{-3}$	1.0139	1.9440$\times {10}^{-3}$	1.0063
480	1.9855$\times {10}^{-3}$	1.0069	9.7031$\times {10}^{-4}$	1.0025
960	9.9040$\times {10}^{-4}$	1.0034	4.8479$\times {10}^{-4}$	1.0011
1920	4.9461$\times {10}^{-4}$	1.0017	2.4231$\times {10}^{-4}$	1.0005

. We observe that for space-dependent diffusion, the scheme in (41) is a little bit more precise in comparison with the one in [6]; see Table 2 and Table 4. However, the scheme in (41) is not only second-order accurate but also preserves the non-negativity of the solution, regardless of whether a convection term is present in (1) or not [30].

Example 2

(Inverse problems: discrete measurements). We consider the same test example as in Example 1 and recover the functions $g_0$ or/and $g_L$. We take measurements from the numerical solution of the direct problem at each time layer of the mesh and at fixed spatial grid node(s). Since in our test example, the exact functions $g_0$ and $g_L$ are known, we may compute the error between the exact and recovered ones, computed by (42)–(45) or (49)–(53). Denote by $g_0^r$ and $g_L^r$ the numerically recovered functions $g_0$ and $g_L$, respectively, and introduce the errors $e_0^n=g_0\left(t_n\right)-g_0^r\left(t_n\right)$, $e_L^n=g_L\left(t_n\right)-g_L^r\left(t_n\right)$. To verify the efficiency of the developed numerical methods for solving IP1–IP6, we define

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle {ϵ}_{\infty }\left(g_0\right)=\underset{1\le n\le N}{max}\left|e_0^n\right|,{ϵ}_2\left(g_0\right)=\sqrt{\tau \sum _{n=1}^N{\left(e_0^n\right)}^2},}\hfill \\ & {\displaystyle {ϵ}_{\infty }\left(g_L\right)=\underset{1\le n\le N}{max}\left|e_L^n\right|,{ϵ}_2\left(g_L\right)=\sqrt{\tau \sum _{n=1}^N{\left(e_L^n\right)}^2}.}\hfill \end{array}\end{array}$$

(55)

To estimate the numerical solution, we use the same formulas as in (55), but now $E_i^n={\mathrm{u}}_i^n-{\tilde{u}}_i^n$, where ${\mathrm{u}}_i^n$ is the solution of the inverse problem (42)–(45) or (49)–(53) at grid node $(x_i,t_n)$.

First, we consider IP1 and IP2 for recovering the function $g_0$ in the left (Robin’s) BC. The computations are performed for $I=480$, $\tau =2h$, and $T=1$. In

Table 5. Errors of the numerical solution of the IP1, IP2; Example 2.

${\mathit{x}}^{*}$	${\mathcal{E}}_{\mathsf{\infty }}$	${\mathcal{E}}_2$	${\mathit{ϵ}}_{\mathsf{\infty }}\left({\mathit{g}}_0\right)$	${\mathit{ϵ}}_2\left({\mathit{g}}_0\right)$
0.50	1.4040$\times {10}^{-12}$	1.4735$\times {10}^{-13}$	1.5858$\times {10}^{-10}$	5.1343$\times {10}^{-11}$
1.25	1.2267$\times {10}^{-9}$	1.0814$\times {10}^{-10}$	1.5917$\times {10}^{-8}$	6.9252$\times {10}^{-9}$
2.50	4.7723$\times {10}^{-8}$	4.4621$\times {10}^{-9}$	6.5239$\times {10}^{-7}$	2.8830$\times {10}^{-7}$

, we present errors of the function determined by (42)–(45), function $g_0$, and solution $\mathrm{u}$, using the measurements at spatial point $x^{*}$ located in some of the regions ${\mathsf{\Omega }}_s$, $s=1,2,3$. The results show that the solution of the numerical approach (42)–(45) recovers the function $g_0$ and the numerical solution of the direct problem almost exactly. Better precision is obtained in the case in which the measurements are in the adjacent region ${\mathsf{\Omega }}_1$ to the left BC.

Next, we examine IP3–IP6. We take measurements at two spatial grid nodes $x_1^{*}$, $x_2^{*}$. In

Table 6. Errors of the numerical solutions of the IP3–IP6; Example 2.

${\mathit{x}}_1^{*}$	${\mathit{x}}_2^{*}$	${\mathcal{E}}_{\mathsf{\infty }}$	${\mathcal{E}}_2$	${\mathit{ϵ}}_{\mathsf{\infty }}\left({\mathit{g}}_0\right)$	${\mathit{ϵ}}_2\left({\mathit{g}}_0\right)$	${\mathit{ϵ}}_{\mathsf{\infty }}\left({\mathit{g}}_{\mathit{L}}\right)$	${\mathit{ϵ}}_2\left({\mathit{g}}_{\mathit{L}}\right)$
0.50	1.75	1.5197$\times {10}^{-12}$	2.9015$\times {10}^{-13}$	1.9315$\times {10}^{-10}$	5.7772$\times {10}^{-11}$	1.5197$\times {10}^{-12}$	6.6679$\times {10}^{-13}$
1.25	2.50	1.3183$\times {10}^{-9}$	1.1228$\times {10}^{-10}$	1.8044$\times {10}^{-8}$	7.2189$\times {10}^{-9}$	1.5099$\times {10}^{-14}$	5.4342$\times {10}^{-15}$
1.25	1.75	1.4452$\times {10}^{-9}$	1.1143$\times {10}^{-10}$	1.9445$\times {10}^{-8}$	7.0999$\times {10}^{-9}$	1.4988$\times {10}^{-12}$	5.4161$\times {10}^{-13}$
0.50	2.50	1.4908$\times {10}^{-12}$	1.7105$\times {10}^{-13}$	1.6941$\times {10}^{-10}$	6.5879$\times {10}^{-11}$	3.1086$\times {10}^{-14}$	1.1447$\times {10}^{-14}$
0.50	0.80	9.9528$\times {10}^{-10}$	1.6202$\times {10}^{-10}$	1.8100$\times {10}^{-10}$	6.0038$\times {10}^{-11}$	9.9528$\times {10}^{-10}$	4.4503$\times {10}^{-10}$
2.15	2.80	1.7992$\times {10}^{-8}$	1.7575$\times {10}^{-9}$	2.4744$\times {10}^{-7}$	1.1331$\times {10}^{-7}$	3.5527$\times {10}^{-15}$	1.4183$\times {10}^{-15}$

, we present the errors of the recovered function by (49)–(53), functions $g_0$, $g_L$, and solution $\mathrm{u}$ for different locations of the points of measurement. We observe that the numerical method (49)–(53) recovers the function $g_0$, $g_L$ and the numerical solution of the direct problem almost exactly. Better precision of the solution is obtained when one of the measurement points is in the region ${\mathsf{\Omega }}_1$ and the other is in ${\mathsf{\Omega }}_3$.

Therefore, the order of convergence of the solution of the numerical methods (42)–(45) and (49)–(53) is the same as the one of the numerical solution of the direct problem. Note that, taking the measurements close to the interface nodes, the recovery is not successful. The closer the measurement point(s) are to the boundary or boundaries where the boundary condition(s) are being recovered, the higher the accuracy of the restored boundary function(s) and solution.

Example 3

(Inverse problem: noisy measurements). We compute the solution of the inverse problems for the same test example as in Example 2 (left Robin’s BC and right Dirichlet BC) to restore the functions $g_0$ and $g_L$, but now we add noise to the discrete observations as follows:

$$\begin{array}{c}\hfill {\psi }_0^n=u^n(x_1^{*},t_n)+2{\rho }_1({\sigma }_1^n-0.5)u^n(x_1^{*},t_n),\\ \hfill {\psi }_M^n=u^n(x_2^{*},t_n)+2{\rho }_2({\sigma }_2^n-0.5)u^n(x_2^{*},t_n),\end{array}$$

where $u^n(x_i^{*},t_n)$, $i=1,2$ are numerical solutions of the direct problem at time layers $t_n$ and at spatial nodes $x_i^{*}$, ${\rho }_i$ are noise levels and ${\sigma }_i^n$ are random values, uniformly distributed on the interval $[0,1]$. We use polynomial curve fitting of degree 5 to smooth the data. The computations are performed for $I=480$ and $\tau =2h$ for $T=1$.

In Figure 1, we present the exact and recovered functions $g_0$, $g_L$ for $x_1^{*}=0.5$ and $x_2^{*}=2.5$ and ${\rho }_1={\rho }_2=0.01$. In this case, the errors of the solution obtained by numerically solving IP3, compared with the discrete solution of the direct problem, are ${\mathcal{E}}_{\infty }=2.7983\times {10}^{-2}$ and ${\mathcal{E}}_2=2.6706\times {10}^{-3}$. In Figure 2 we present determined boundary functions for ${\rho }_1=0.03$ and ${\rho }_2=0.05$. Now, the corresponding errors are ${\mathcal{E}}_{\infty }=8.3858\times {10}^{-2}$ and ${\mathcal{E}}_2=1.0642\times {10}^{-2}$

In Figure 3, we depict the exact and recovered functions $g_0$, $g_L$ for $x_1^{*}=0.5$, $x_2^{*}=1.9$ and ${\rho }_1=0.03$, ${\rho }_2=0.05$, obtained by numerically solving IP4. The corresponding errors of the solution are ${\mathcal{E}}_{\infty }=1.3969\times {10}^{-1}$ and ${\mathcal{E}}_2=1.4772\times {10}^{-2}$.

Next, we examine the efficiency of the numerical method for recovering nonsmooth function $g_0\left(t\right)$ and discontinuous function $g_L\left(t\right)$. In the above test problem, we take

$$g_0\left(t\right)=\left\{\begin{array}{cc}3t+1,\hfill & 0\le t\le 0.5,\hfill \\ 5(T-t),\hfill & t\le 0.5\le 1,\hfill \end{array}\right.g_L\left(t\right)=\left\{\begin{array}{cc}{\displaystyle 1,}\hfill & 0\le t<0.25,\hfill \\ 5t,\hfill & 0.25<t<0.75,\hfill \\ 4,\hfill & 0.75<t\le 1.\hfill \end{array}\right.$$

In Figure 4, we provide t the exact and recovered functions $g_0$, $g_L$ for $x_1^{*}=0.5$ and $x_2^{*}=1.9$ and ${\rho }_1=0.03$, ${\rho }_2=0.05$, obtained by numerically solving IP4. In Figure 5, we present the solution and error of the solution. In this case, ${\mathcal{E}}_{\infty }=5.3997\times {10}^{-1}$ and ${\mathcal{E}}_2=5.8574\times {10}^{-2}$.

The experiments with noisy measurements illustrate that the numerical method for solving inverse problems recovers the smooth, nonsmooth, and discontinuous functions $g_0$, $g_L$, and the corresponding solution u with optimal accuracy, but the Dirichlet BC (namely the function $g_L$) is recovered more precisely in comparison with Robin’s BC (namely the function $g_0$).

7. Conclusions and Perspectives

In this study, we investigated heat conduction problems involving layered materials. We formulated multilayer parabolic problems incorporating typical interface relations from physics, mechanics, and heat-mass transfer. We established well-posedness in appropriate Sobolev spaces for direct (forward) symmetric (self-adjoint) initial boundary value problems. Additionally, we presented several inverse problems related to unknown external boundary conditions based on different point observations. We constructed a second-order exponential difference scheme for the direct problem and employed special solution decomposition techniques to numerically address the inverse problems. Simulations of test examples demonstrate that the proposed method achieves a second-order convergence in space and a first-order convergence in time for the direct problem. The numerical method for solving the inverse problem effectively recovers functions $g_0$ and $g_L$ and provides high-accuracy numerical solutions for the direct problem, nearly exact in the cases of exact measurements, located at points close to the corresponding boundary. The position of measurement points significantly affects the precision of the method, with better results obtained when measurements are taken in domains adjacent to the boundary. Simulations with perturbed measurements show that the solution achieves optimal accuracy, particularly in recovering Dirichlet boundary conditions. The limitation of the method is related to the complications that arise in the extension of the method to problems with curvilinear interfaces, especially in higher-dimensional cases. In future work, we aim to develop approaches for solving higher-dimensional multilayer nonsymmetric diffusion–convection problems with nonlinear reaction terms. Additionally, we plan to extend the inverse problem solution methods to incorporate other types of measurements, such as integral observations.