Open AccessArticle

Inverse Boundary Conditions Interface Problems for the Heat Equation with Cylindrical Symmetry

Miglena N. Koleva

^1,*

and

Lubin G. Vulkov

^2,*

Department of Mathematics, Faculty of Natural Sciences and Education, “Angel Kanchev” University of Ruse, 8 Studentska Str., 7017 Ruse, Bulgaria

Department of Applied Mathematics and Statistics, Faculty of Natural Sciences and Education, “Angel Kanchev” University of Ruse, 8 Studentska Str., 7017 Ruse, Bulgaria

Authors to whom correspondence should be addressed.

Symmetry 2024, 16(8), 1065; https://doi.org/10.3390/sym16081065

Submission received: 12 July 2024 / Revised: 10 August 2024 / Accepted: 16 August 2024 / Published: 18 August 2024

(This article belongs to the Special Issue Symmetry in Mathematical Theory and Simulation Methods for Backward Problems)

Download

Browse Figures

Versions Notes

Abstract

In this paper, we study inverse interface problems with unknown boundary conditions, using point observations for parabolic equations with cylindrical symmetry. In the one-dimensional, two-layer interface problem, the left interval

0 < r < l_{1}

, i.e., the zero degeneracy, causes serious solution difficulty. For this, we investigate the well-posedness of the direct (forward) problem. Next, we formulate and solve five inverse boundary condition problems for the interface heat equation with cylindrical symmetry from internal measurements. The finite volume difference method is developed to construct second-order schemes for direct and inverse problems. The correctness of the proposed numerical solution decomposition algorithms for the inverse problems is discussed. Several numerical examples are presented to illustrate the efficiency of the approach.

Keywords:

parabolic equations with cylindrical symmetry; inverse boundary condition problems; well-posedness; solution decomposition; finite difference schemes

MSC:

65M08; 65M32

1. Introduction

Inverse boundary problems for parabolic problems determine the unknown boundary conditions of the equations from either final time measurements (backward problems) or local space or time measurements of the solutions. Such problems play a crucial role in applied mathematics and physics. For example, inverse and backward source and boundary problems are used for the identification of coefficient and boundary sources of pollution sources affecting water bodies such as rivers and lakes, heat-mass transfer, elasticity, heat conduction systems, medicine, etc. (see, e.g., [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]).

It is challenging to solve inverse problems, even numerically, because, in general, they are ill-posed [2,16,17,18,19,20]. In [9], the authors provide theoretical reasons why the aquifer identification problem can be expected to be numerically unstable. Many papers have studied the inverse problem of identifying an unknown boundary in a parabolic equation. A class of such problems plays a key role in applied mathematics, physics, mechanics, and biology (see, e.g., [2,8,12,15,19,21,22,23]).

The uniqueness of the solution of inverse boundary problems in parabolic and time-fractional diffusion equations by observation of the temperature at a fixed point on the boundary is investigated in [15,24]. In [1], a meshless method is constructed for solving inverse boundary source problems for time-fractional diffusion equations. A semigroup approach is applied in [8,25] for studying inverse problems of recovering the right boundary condition in a linear and quasi-linear parabolic equation under the overspecified data of Dirichlet and Neumann type on the left boundary. A numerical method, based on a decomposition technique for solving the inverse problem for identifying the Dirichlet boundary condition in the heat conduction equation is proposed in [26].

Inverse and backward source and boundary problems are often used to model heat-mass transfer problems in composite materials. For example, in [23], the authors construct a numerical method, based on fundamental solutions and the Tikhonov regularization approach, for recovering a moving boundary in a one-dimensional heat equation on a multilayer domain.

In our previous works [27,28,29], we developed numerical methods, based on the decomposition technique or reduction method combined with the Saul’yev scheme, for recovering external boundary conditions in parabolic, time-fractional parabolic, and parabolic-hyperbolic problems, defined on disjoint domains. In order to solve inverse problems, we usually have to solve well-posed direct (forward) ones as well. There are many studies in this direction. We will mention some of them.

A global well-posedness result for 3D micropolar fluid equations is obtained in [30]. Lyapunov stability over a finite time interval provides well-posedness of the considered differential problem (see [31]). Stronger local and global finite-time stability for time-varying nonlinear systems are investigated. In [32], using comparison principles and the quasilinearization approach, the authors obtain monotone iterative sequences that uniformly and monotonically converge to the unique solution of integro-differential equations with nonlinear boundary conditions.

The work [33] provides an exact analysis of the stability of a multilayer diffusion-reaction problem, where the number of layers can be extensive. In [34], the authors construct a second-order numerical scheme for an elliptic interface system in the case of perfect and imperfect contact. The authors of [35] present a theoretical and computational investigation of an elliptic equation on a domain that consists of two sub-domains with substantially different thermal conductivities and fiber structures. A second-order accurate

L_{2}

norm difference scheme for the axisymmetric Poisson equation is constructed and analyzed in [36]. Semi-analytical solutions of axisymmetric reaction–diffusion equations, describing heat conduction in composite media, are developed in [37].

In particular, inverse and backward source and boundary problems describe mass transfer processes in columns with cylindrical symmetry, taking into account the longitudinal mixing (see, e.g., [10]). Inverse source problems also have applications in multi-wave imaging, geophysics, [3,10,38] and identification of the rate of the outflow in oil reservoirs with walls [4]. Additionally, they have applications in numerous cylindrical constructs, such as pipelines for transporting substantial quantities of natural resources like gas and oil [39].

Recently, many models of biomechanics use cylindrical regions, such as bone and blood vessels. A scheme for determining the inhomogeneous elastic properties of an isotropic cylinder is developed in [5]. An inverse technique for estimating the inlet conditions of fluid in a thick pipe is adapted in [22].

The backward space-dependent source problem from noisy final data for the helium production-diffusion equation is studied in [38]. A coefficient inverse problem for the two-dimensional transient heat equation in a cylindrical coordinate system using the finite difference method is solved in [21]. An inverse problem for reconstructing a space-dependent source in the heat equation of a polar coordinate system from finite temperature measurements is studied in [40]. A modified quasi-boundary value method is developed in [7] to solve an inverse heat conduction problem with spherical symmetry for recovering the internal surface temperature distribution of a hollow sphere from observations at a fixed position inside it. A higher-order difference scheme is studied for the problem, formulated in cylindrical and spherical coordinate systems in [41].

In this article, we develop efficient numerical methods for solving inverse problems for identifying external Dirichlet boundary conditions in parabolic interface problems with cylindrical symmetry.

The main difficulties arise from the degeneration at

r = 0

, the multi-layered structure of the domain, and the discontinuity of the coefficients (interface problem) across the different layers. They are overcome by constructing appropriate discretizations and decomposition techniques. Moreover, the formulated five inverse problems are reduced, as the method is the same regardless of whether the measurements are taken in the adjacent layer or further from the unknown boundary.

The rest of the paper is organized as follows: In Section 2, the direct (forward) and inverse problems are formulated. The well-posedness of the direct problem is proved in Section 3. Finite volume second-order difference schemes for the direct problem are derived in Section 4. A solution decomposition method for numerical solution of the inverse problems is proposed and studied for correctness in Section 5. The results from numerical experiments are presented and discussed in Section 6. The paper is finalized with some conclusions.

2. Direct and Inverse Problems

In this section, we formulate the basic interface direct (forward) problem for the heat equation with cylindrical symmetry.

2.1. Direct (Forward) Problem

We consider the interface initial-boundary value problem (IBVP):

\begin{matrix} \frac{\partial u_{1}}{\partial t} = d_{1} (\frac{\partial^{2} u_{1}}{\partial r^{2}} + \frac{1}{r} \frac{\partial u_{1}}{\partial r}) + f_{1} (r, t), (r, t) \in Q_{1 T} = Ω_{1} \times (0, T), \end{matrix}

(1)

\begin{matrix} \frac{\partial u_{2}}{\partial t} = d_{2} (\frac{\partial^{2} u_{2}}{\partial r^{2}} + \frac{1}{r} \frac{\partial u_{2}}{\partial r}) + f_{2} (r, t), (r, t) \in Q_{2 T} = Ω_{2} \times (0, T), \end{matrix}

(2)

where

u_{i} (r, t), i = 1, 2

is the solution (temperature or concentration, etc.) at position

r \in Ω_{i}

Ω_{1} = (0, l_{1})

r_{2} = (l_{1}, l_{2})

l_{2} = L

and time t,

d_{i}

is the diffusion coefficient in the i-th layer,

i = 1, 2

The initial conditions are given by

u_{i} (r, 0) = u_{i 0} (r), i = 1, 2

(3)

and the external boundary conditions (BCs) are

\begin{matrix} u_{1} (0, t) = g_{0} (t), 0 < t < T, \end{matrix}

(4)

\begin{matrix} u_{2} (L, t) = g_{L} (t), 0 < t < T . \end{matrix}

(5)

Here,

g_{0} (t)

g_{L} (t)

are specified time-dependent continuous functions.

The typical assumption is that the solution and the diffusive flux are continuous at the interface, meaning

{[u]}_{r = l_{1}} = u_{2} (l_{1}, t) - u_{1} (l_{1}, t) = 0,

(6)

{[d \frac{\partial u}{\partial r}]}_{r = l_{1}} = d_{2} \frac{\partial u_{r}}{\partial r} (l_{1}, t) - d_{1} \frac{\partial u_{r}}{\partial r} (l_{1}, t) = 0,

(7)

for

t \in (0, T)

Let us note that, in mechanics, physics, biology, etc., many real processes require different relations than those in (6) and (7) (see e.g., [42,43]).

The version in Cartesian coordinates of the problem (1)–(7) is the heat equation problem with an interface [42].

The direct (forward) problem is to find the solution

u_{1} (x, t)

u_{2} (x, t)

of the IBVP (1)–(7) with given coefficients, right-hand sides, initial, boundary, and interface conditions.

2.2. Inverse Problems

We consider inverse problems where the functions

g_{0} (t)

and

g_{L} (t)

are unknown and must be identified using some overspecified data. Depending on the observation data, the inverse problems (IP) for determining the boundary functions in (4) and (5) can be formulated as follows:

IP1: Given observation

u_{1} (r_{1}^{*}, t) = Ψ_{1} (t), r_{1}^{*} \in Ω_{1}

. Find

(u_{1}, u_{2}, g_{0} (t))

;

IP2: Given observation

u_{2} (r_{2}^{*}, t) = Ψ_{2} (t), r_{2}^{*} \in Ω_{2}

. Find

(u_{1}, u_{2}, g_{L} (t))

;

IP3: Given observation

u_{1} (r_{1}^{*}, t) = Ψ_{1} (t), r_{1}^{*} \in Ω_{1}

. Find

(u_{1}, u_{2}, g_{L} (t))

;

IP4: Given observation

u_{2} (r_{2}^{*}, t) = Ψ_{2} (t), r_{2}^{*} \in Ω_{2}

. Find

(u_{1}, u_{2}, g_{0} (t))

;

IP5: Given observation

u_{i} (r_{i}^{*}, t) = Ψ_{i} (t), r_{i}^{*} \in Ω_{j}, i = 1, 2

. Find

(u_{1}, u_{2}, g_{0} (t), g_{L} (t))

3. Well-Posedness of the Direct Problem

In this section, an analytical study of the direct problem is carried out to provide a rigorous proof of the well-posedness of the direct problem (1)–(7). The main difficulties in solving the direct problem come from the singularity at

r = 0

and the interface.

Theorem 1.

Let

lim_{r \to 0} u_{1} (r, t)

be bounded,

g_{i} \in C (0, T)

f_{i} \in Q_{i T}

u_{i 0} (x) \in C (Ω_{i})

i = 1, 2

. Then, there exists a unique solution

u = {u_{1}, u_{2}}

u_{1} \in C^{2} (Ω_{1}) \cap C^{1} ({r | r = l_{1}})

u_{2} \in C^{2} (Ω_{2}) \cap C^{1} ({r | r = l_{1}}) \cap C^{1} ({r | r = R_{2}})

of the problem (1)–(7) that continuously depends on the input data.

Proof.

We consider the particular case of homogeneous Dirichlet BCs (4) and (5)

u_{1} (0, t) = g_{0} (t) = 0, u_{2} (L, t) = g_{L} (t) = 0 .

First, we study the uniqueness of the solution, assuming existence. We multiply Equation (1) by

r u_{1} (r, t)

and Equation (2) by

r u_{2} (r, t)

. Next, we integrate the new equations over the intervals

Ω_{1}

and

Ω_{2}

, respectively, and add up the results to obtain:

\frac{\partial}{\partial t} U (t) + d_{1} \int_{0}^{l_{1}} r {(\frac{\partial u_{1}}{\partial r})}^{2} d r + d_{2} \int_{l_{1}}^{l_{2}} r {(\frac{\partial u_{2}}{\partial r})}^{2} d r = \int_{0}^{l_{1}} r u_{1} f_{1} d r + \int_{l_{1}}^{l_{2}} r u_{2} f_{2} d r,

(8)

where

U (t) = \frac{1}{2} \int_{0}^{l_{1}} r u_{1}^{2} (r, t) d r + \frac{1}{2} \int_{l_{1}}^{l_{2}} r u_{2}^{2} (r, t) d r .

(9)

Let us note that for the function

u_{2} (r, t)

, it is easy to derive a Friedrich-type inequality (see e.g., [22,44]),

\int_{l_{1}}^{l_{2}} r u_{2}^{2} d r \leq C \int_{l_{1}}^{l_{2}} r {(\frac{\partial u_{2}}{\partial r})}^{2} d r, C = \frac{l_{1}^{2}}{4} (\frac{l_{2}^{2}}{l_{1}^{2}} - 2 ln \frac{l_{2}}{l_{1}} - 1) > 0 .

However, because of the singularity around

r = 0

in the Equation (1), it is not easy to obtain such an estimate for

u_{1} (r, t)

. Fortunately, in [45], the following inequality was derived:

\int_{0}^{l_{1}} r u_{1}^{2} d r + \int_{l_{1}}^{l_{2}} r u_{2}^{2} d r \leq C (d_{1} \int_{0}^{l_{1}} r {(\frac{\partial u_{1}}{\partial r})}^{2} d r + d_{2} \int_{l_{1}}^{l_{2}} r {(\frac{\partial u_{2}}{\partial r})}^{2} d r),

(10)

where

C > 0

is a constant that is not dependent on the functions

u_{1}, u_{2}

Now, using (10), we deduce that the left-hand side of the equality (8) is greater than or equal to

t_{0}

\frac{\partial U}{\partial t} \leq \frac{1}{C} (\int_{0}^{l_{1}} r u_{1}^{2} d r + \int_{l_{1}}^{l_{2}} r u_{2}^{2} d r) .

For the right-hand side of (8), using (10) again, we obtain:

\begin{matrix} \int_{0}^{l_{1}} r u_{1} f_{1} d r + \int_{l_{1}}^{l_{2}} r u_{2} f_{2} d r \leq {(\int_{0}^{l_{1}} r f_{1}^{2} d r)}^{1 / 2} {(\int_{0}^{l_{1}} r u_{1}^{2} d r)}^{1 / 2} \\ + {(\int_{l_{1}}^{l_{2}} r f_{2}^{2} d r)}^{1 / 2} {(\int_{l_{1}}^{l_{2}} r u_{2}^{2} d r)}^{1 / 2} . \\ \max (l_{1} . {\hat{f}}_{1}, \sqrt{l_{2}^{2} - l_{1}^{2}} {\hat{f}}_{2}) [{(\frac{1}{2} \int_{0}^{l_{1}} r u_{1}^{2} d r)}^{1 / 2} + \frac{1}{2} {(\int_{l_{1}}^{l_{2}} r u_{2}^{2} d r)}^{1 / 2}) \\ \leq \max (l_{1} {\hat{f}}_{1}, \sqrt{l_{2}^{2} - l_{1}^{2}} {\hat{f}}_{2}) (2 \sqrt{2} \sqrt{U}) = F (t) \sqrt{U}, \end{matrix}

(11)

where

{\hat{f}}_{i} (t) = \sqrt{U max_{Ω_{i}} | f_{1} (x, t) |}, i = 1, 2 .

Therefore, from (8) and (11), we obtain the inequalities

\frac{\partial}{\partial t} U + c U \leq F (t) U . c = \frac{1}{C} min (d_{1}, d_{2}),

which implies

U (t) \leq {(\sqrt{U_{0}} + \frac{1}{2} \int_{0}^{t} F (τ) e^{c τ} d τ)}^{2} e^{- 2 c t},

(12)

where

U_{0} = U (0) = \frac{1}{2} \int_{0}^{l_{1}} r u_{10}^{2} (r) d r + \frac{1}{2} \int_{l_{1}}^{l_{2}} r u_{20}^{2} (r) d r .

Therefore, from (9), we have:

\int_{0}^{l_{1}} r u_{1}^{2} (r, t) d r \leq 2 U (t), \int_{l_{1}}^{l_{2}} r u_{2}^{2} (r, t) d r \leq 2 U (t),

where

U (t)

satisfies (12).

Next, we will prove the boundedness of the integrals

\int_{0}^{l_{1}} r {(\frac{\partial u_{1}}{\partial r})}^{2} d r, C \int_{l_{1}}^{l_{2}} r {(\frac{\partial u_{2}}{\partial r})}^{2} d r .

For this aim, we take the square of Equations (1) and (2), and then integrate the results:

\begin{matrix} \int_{0}^{t} \int_{0}^{l_{1}} r {(\frac{\partial u_{1}}{\partial t})}^{2} d r d s - 2 d_{1} \int_{0}^{t} \int_{0}^{l_{1}} \frac{\partial u_{1}}{\partial t} \frac{\partial}{\partial r} (r \frac{\partial u_{1}}{\partial r}) d r d s \\ + \int_{0}^{t} \int_{0}^{l_{1}} d_{1}^{2} . \frac{1}{r} {(\frac{\partial}{\partial r} (r \frac{\partial u_{1}}{\partial r}))}^{2} d r d s = \int_{0}^{t} \int_{0}^{l_{1}} r . f_{1}^{2} d r d s . \\ \int_{0}^{t} \int_{l_{1}}^{l_{2}} r {(\frac{\partial u_{2}}{\partial t})}^{2} d r d s - 2 d_{2} \int_{0}^{t} \int_{l_{1}}^{l_{2}} \frac{\partial u_{2}}{\partial t} \frac{\partial}{\partial r} (r \frac{\partial u_{2}}{\partial r}) d r d s \\ + \int_{0}^{t} \int_{l_{1}}^{l_{2}} d_{2}^{2} \frac{1}{r} {(\frac{\partial}{\partial r} (r \frac{\partial u_{2}}{\partial r}))}^{2} d r d s . \end{matrix}

Summing these equalities, and using integration by parts along with the initial, boundary, and interface conditions, we obtain

\begin{matrix} \int_{0}^{t} \int_{l_{1}}^{l_{2}} r {(\frac{\partial u_{1}}{\partial t})}^{2} d r d s + \int_{0}^{t} r {(\frac{\partial u_{2}}{\partial t})}^{2} d r d s + \int_{0}^{l_{1}} r {(\frac{\partial u_{1}}{\partial r})}^{2} d r + \int_{l_{1}}^{l_{2}} r {(\frac{\partial u_{2}}{\partial r})}^{2} d r \\ + \int_{0}^{t} d_{1}^{2} \frac{1}{r} {(\frac{\partial}{\partial r} (r \frac{\partial u_{1}}{\partial r}))}^{2} d r d s + \int_{0}^{t} d_{2}^{2} \frac{1}{r} {(\frac{\partial}{\partial r} (r \frac{\partial u_{2}}{\partial r}))}^{2} \\ = \int_{0}^{t} \int_{0}^{l_{1}} r f_{1}^{2} d r d s + \int_{0}^{t} \int_{l_{1}}^{l_{2}} r f_{2}^{2} d r d s . \end{matrix}

Therefore,

\int_{Ω_{i}} r {(\frac{\partial u_{1}}{\partial r})}^{2} d r \leq \frac{1}{d_{1}} \int_{0}^{t} \int_{l_{1}}^{l_{2}} r f_{1}^{2} d r d s + \frac{1}{d_{2}} \int_{0}^{t} \int_{l_{1}}^{l_{2}} r f_{2}^{2} d r d s .

Hence, the uniqueness and continuous dependence of the solution of the problem (1)–(7) on the input data follows.

The local existence of the solution of the problem (1)–(7) follows from the general theory of parabolic equations (see e.g., [22,44]). The global existence is a direct consequence of the a priory estimates above. □

4. Numerical Solution of the Direct Problem

First, we solve the direct problem (1)–(7). We rewrite the Equations (1) and (2) in divergent form

\begin{matrix} \frac{\partial u_{1}}{\partial t} = \frac{d_{1}}{r} \frac{\partial}{\partial r} (r \frac{\partial u_{1}}{\partial r}) + f_{1} (r, t), (r, t) \in Q_{1 T} = Ω_{1} \times (0, T), \end{matrix}

(13)

\begin{matrix} \frac{\partial u_{2}}{\partial t} = \frac{d_{2}}{r} \frac{\partial}{\partial r} (r \frac{\partial u_{2}}{\partial r}) + f_{2} (r, t), (x, t) \in Q_{2 T} = Ω_{2} \times (0, T), \end{matrix}

(14)

r = 0

, we impose the condition [42]

| u_{1} (0) | < \infty,

(15)

or equivalently,

lim_{r \to 0} r \frac{\partial u_{1}}{\partial r} = 0 .

(16)

This condition can be replaced by

\frac{\partial u_{1}}{\partial r} (0, t) = 0 .

(17)

Therefore, at

r = 0

, in view of (15), we may consider either the Dirichlet boundary condition (4) or the natural boundary condition (17), based on (16).

Let us define the uniform spatial and temporal meshes

{\bar{ω}}_{h} = {r_{i} = i h, j = 0, 1, \dots, J, J h = L}, {\bar{ω}}_{τ} = {t_{n} = n τ, n = 0, 1, \dots, N, N τ = T} .

Denote by

u_{i, j}^{n}

the solution

u_{i}

at the grid node

(r_{j}, t_{n}) \in {\bar{ω}}_{h} \times {\bar{ω}}_{τ}

and by

u_{i, j + 1 / 2}^{n}

, the solution

u_{i}

at the grid node

(r_{j + 1 / 2}, t_{n})

, where

r_{j + 1 / 2} = r_{j} + h / 2

is a node of the dual mesh

0 = r_{- 1 / 2} < r_{1 / 2} < r_{3 / 2} \dots < r_{J - 1 / 2} < r_{J + 1 / 2} = L

In order to derive a conservative scheme and overcome the degeneracy at

r = 0

, we follow [42] and apply the finite volume method to the problem (13), (14), (3), (17), and (5). Let

l_{1}

be a grid node, where

l_{1} = r_{j^{*}} = j^{*} h

. At the inner spatial nodes, multiplying (13) and (14) by r, integrating over the volumes

(r_{j - 1 / 2}, r_{j + 1 / 2})

j = 1, 2, \dots, J - 1

, and dividing by

h r_{i}

, we obtain

\begin{matrix} \frac{1}{r_{i} h} \int_{r_{j - 1 / 2}}^{r_{j + 1 / 2}} r \frac{\partial u_{1}}{\partial t} d r - \frac{d_{1}}{r_{j} h} (r_{j + 1 / 2} \frac{\partial u_{1}}{\partial r} |_{j + 1 / 2} - r_{j - 1 / 2} \frac{\partial u_{1}}{\partial r} |_{j - 1 / 2}) \\ = \frac{1}{r_{j} h} \int_{r_{j - 1 / 2}}^{r_{j + 1 / 2}} f_{1} (r, t) r d r, j = 1, 2, \dots, j^{*} - 1, \end{matrix}

(18)

\begin{matrix} \frac{1}{r_{j} h} \int_{r_{j - 1 / 2}}^{r_{j}} r \frac{\partial u_{1}}{\partial t} d r + \int_{r_{j}}^{r_{j + 1 / 2}} \frac{\partial u_{2}}{\partial t} d r - \frac{1}{r_{j} h} (d_{2} r_{j + 1 / 2} \frac{\partial u_{2}}{\partial r} |_{j + 1 / 2} - d_{1} r_{j - 1 / 2} \frac{\partial u_{1}}{\partial r} |_{j - 1 / 2}) \\ = \frac{1}{r_{j} h} \int_{r_{j - 1 / 2}}^{r_{j}} f_{1} (r, t) r d r + \frac{1}{r_{j} h} \int_{r_{j}}^{r_{j + 1 / 2}} f_{2} (r, t) r d r, j = j^{*}, \end{matrix}

(19)

\begin{matrix} \frac{1}{r_{j} h} \int_{r_{j - 1 / 2}}^{r_{j + 1 / 2}} r \frac{\partial u_{2}}{\partial t} d r - \frac{d_{2}}{r_{j} h} (r_{j + 1 / 2} \frac{\partial u_{2}}{\partial r} |_{j + 1 / 2} - r_{j - 1 / 2} \frac{\partial u_{2}}{\partial r} |_{j - 1 / 2}) \\ = \frac{1}{r_{j} h} \int_{r_{j - 1 / 2}}^{r_{j + 1 / 2}} f_{2} (r, t) r d r, j = j^{*} + 1, \dots, J - 1 . \end{matrix}

(20)

Approximating the integrals in (18)–(20) by the midpoint quadrature, the fluxes in (18)–(20) by finite differences are

r_{i - 1 / 2} \frac{\partial u_{i}}{\partial r} = r_{i - 1 / 2} \frac{u_{i, j} - u_{i, j - 1}}{h},

and applying implicit time-stepping, we derive

\begin{matrix} \frac{u_{1, j}^{n + 1} - u_{1, j}^{n}}{τ} - \frac{d_{1}}{r_{j} h} (r_{j + 1 / 2} \frac{u_{1, j + 1}^{n + 1} - u_{1, j}^{n + 1}}{h} - r_{j - 1 / 2} \frac{u_{1, j}^{n + 1} - u_{1, j - 1}^{n + 1}}{h}) \\ = f_{1, j}^{n + 1}, j = 1, 2, \dots, j^{*} - 1 \end{matrix}

(21)

\begin{matrix} \frac{u_{1, j}^{n + 1} - u_{1, j}^{n}}{2 τ} + \frac{u_{2, j}^{n + 1} - u_{2, j}^{n}}{2 τ} - \frac{1}{r_{j} h} (d_{2} r_{j + 1 / 2} \frac{u_{2, j + 1}^{n + 1} - u_{2, j}^{n + 1}}{h} - d_{1} r_{j - 1 / 2} \frac{u_{2, j}^{n + 1} - u_{2, j - 1}^{n + 1}}{h}) \\ = \frac{1}{2} (f_{1, j}^{n + 1} + f_{2, j}^{n + 1}), j = j^{*}, \end{matrix}

(22)

\begin{matrix} \frac{u_{2, j}^{n + 1} - u_{2, j}^{n}}{2 τ} - \frac{d_{2}}{r_{j} h} (r_{j + 1 / 2} \frac{u_{2, j + 1}^{n + 1} - u_{2, j}^{n + 1}}{h} - r_{j - 1 / 2} \frac{u_{2, j}^{n + 1} - u_{2, j - 1}^{n + 1}}{h}) \\ = f_{2, j}^{n + 1}, j = j^{*} + 1, \dots, J - 1 . \end{matrix}

(23)

To obtain the discretization at

r = 0

, in the case of natural BCs (16) and (17), we consider the interval

I_{0} = (ε, r_{1 / 2} = h / 2)

and multiply (13) by r. Then, we integrate over the volume

I_{0}

to obtain

\int_{ε}^{r_{1 / 2}} r \frac{\partial u_{1}}{\partial t} d r - (r_{1 / 2} \frac{\partial u_{1}}{\partial r} |_{1 / 2} - ε \frac{\partial u_{1}}{\partial r} |_{ε}) = \int_{ε}^{r_{1 / 2}} f_{1} (r, t) r d r .

Using the product approximation formula, we derive

\frac{h}{4} \int_{ε}^{r_{1 / 2}} \frac{\partial u_{1}}{\partial t} d r - (r_{1 / 2} \frac{\partial u_{1}}{\partial r} |_{1 / 2} - ε \frac{\partial u_{1}}{\partial r} |_{ε}) = \frac{h}{4} \int_{ε}^{r_{1 / 2}} f_{1} (r, t) d r,

Now, we apply the rectangle quadrature formula and let

ε \to 0

. In view of (16) and (17), for the left external boundary condition, we obtain

\frac{h}{4} \frac{h}{2} \frac{\partial u_{1}}{\partial t} (0, t) - \frac{h}{2} \frac{\partial u_{1}}{\partial r} |_{1 / 2} = \frac{h}{4} \frac{h}{2} f_{1} (0, t) .

Approximating the derivatives by finite differences and using an implicit time-stepping method, we obtain

\frac{u_{1, 1}^{n + 1} - u_{1, 0}^{n + 1}}{h} = \frac{h}{4} (\frac{u_{1, 0}^{n + 1} - u_{1, 0}^{n}}{τ} - f_{1, 0}^{n + 1}) .

(24)

In the case of Dirichlet BCs (4), the approximation is trivial

u_{1, 0}^{n + 1} = g_{0} (t_{n + 1}) .

(25)

So, we use either (24) or (25) for the computations.

The numerical scheme is completed by the left external boundary condition and initial condition, corresponding to (5) and (3), respectively

\begin{matrix} u_{2, J}^{n + 1} = u_{2} (L, t_{n + 1}), \end{matrix}

(26)

\begin{matrix} u_{i, j}^{0} = u_{i 0} (r_{j}) = \{\begin{matrix} j = 0, 1, \dots, j^{*}, & i = 1, \\ j = j^{*}, j^{*} + 1, \dots, J, & i = 2, \end{matrix} u_{10} (l_{1}) = u_{20} (l_{1}) . \end{matrix}

(27)

Denote

\begin{matrix} u = [u_{1, 0}, u_{1, 1}, \dots, u_{1, j^{*}} = u_{2, j^{*}}, u_{2, j^{*} + 1}, \dots, u_{2, J}], \\ u^{0} = [u_{10} (r_{0}), u_{10} (r_{0}), \dots, u_{10} (r_{j^{*}}) = u_{20} (r_{j^{*}}), u_{10} (r_{j^{*} + 1}), \dots, u_{10} (r_{J})] \end{matrix}

and rewrite the numerical scheme (21)–(23), (24) or (25), (26), (27) in more convenient form. For

n = 0, 1, \dots, N

, we have

\begin{matrix} \begin{array}{l} P_{0} u_{0}^{n + 1} - D_{0} u_{1}^{n + 1} = F_{0}^{n}, \\ - P_{j} u_{j - 1}^{n + 1} + Q_{j} u_{j}^{n + 1} - D_{i} u_{j + 1}^{n + 1} = F_{j}^{n}, j = 1, 2, \dots, J - 1, \\ u_{J}^{n + 1} = u_{2} (L, t_{n + 1}) . \end{array} \end{matrix}

(28)

where

Q_{j} = \frac{1}{τ} + P_{j} + D_{j},

\begin{matrix} P_{j} = \frac{r_{j - 1 / 2}}{r_{j} h} \{\begin{matrix} d_{1}, & j = 1, 2, \dots, j^{*}, \\ d_{2}, & j = j^{*} + 1, \dots, J - 1, \end{matrix}, D_{j} = \frac{r_{j + 1 / 2}}{r_{j} h} \{\begin{matrix} d_{1}, & j = 1, 2, \dots, j^{*} - 1, \\ d_{2}, & j = j^{*}, \dots, J - 1, \end{matrix} \\ F_{j}^{n} = \frac{u_{j}^{n}}{τ} + \{\begin{matrix} f_{1, j}^{n}, & j = 1, 2, \dots, j^{*} - 1, \\ (f_{1, j}^{n} + f_{2, j}^{n}) / 2, & j = j^{*}, \\ f_{2, j}^{n}, & j = j^{*} + 1, \dots, J - 1, \end{matrix} \\ P_{0} = \{\begin{matrix} \frac{1}{τ} + \frac{4}{h^{2}}, & BC (24), \\ 1, & BC (25), \end{matrix} D_{0} = \{\begin{matrix} \frac{4}{h^{2}}, & BC (24), \\ 0, & BC (25), \end{matrix} \\ F_{0}^{n} = \{\begin{matrix} \frac{u_{j}^{n}}{τ} + \frac{h}{4} f_{1, 0}^{n + 1}, & BC (24), \\ g_{0} (t_{n + 1}), & BC (25) . \end{matrix} \end{matrix}

Following the arguments in [42], we may deduce that the discretization (28) has an accuracy order

O (τ + h^{2})

5. Numerical Solution of the Inverse Problems

In this section, we propose a numerical method for solving inverse problems. To obtain numerical discretization for identifying external boundary conditions, we use a decomposition technique.

We begin with solving IP1-IP5. Note that in the case of natural BC (16) and (17), the problems IP1 and IP4 are omitted, as there is nothing to determine. Similarly, IP5 reduces to IP2 and IP3.

Consider IP2 and IP3. For the numerical method we construct, it does not matter whether the measurement is in the adjacent region of

g_{L} (t)

, or not. Therefore, we consider IP2 and IP3 together.

We seek the solution in the form

\begin{matrix} u_{j}^{n + 1} = y_{j}^{n + 1} + g_{L}^{n + 1} v_{j} \end{matrix}

(29)

and substitute it in (28). As a result, we obtain two problems for the unknown functions

y_{j}^{n + 1}

and

v_{j}^{n + 1}

\begin{array}{l} P_{0} y_{0}^{n + 1} - D_{0} y_{1}^{n + 1} = F_{0}^{n}, \\ - P_{j} y_{j - 1}^{n + 1} + Q_{j} y_{j}^{n + 1} - D_{i} y_{j + 1}^{n + 1} = F_{j}^{n}, j = 1, 2, \dots, J - 1, \\ y_{J}^{n + 1} = 0, \\ u_{j}^{0} = u^{0} (j), \end{array}

(30)

\begin{array}{l} P_{0} v_{0} - D_{0} v_{1} = 0, \\ - P_{j} v_{j - 1} + Q_{j} v_{j} - D_{j} v_{j + 1} = 0, j = 1, 2, \dots, J - 1, \\ v_{J} = 1 . \end{array}

(31)

After solving (30) and (31), we find y and v in the whole computational domain. Then, from (29) and the measurement

u (r^{*}) = Ψ (t)

r^{*} \in Ω_{1} \cup Ω_{2}

, we find

g_{L}^{n + 1} = \frac{Ψ^{n + 1} - y_{r^{*}}^{n + 1}}{v_{r^{*}}},

(32)

where

y_{r^{*}}^{n + 1}

and

v_{r^{*}}

are numerical solutions

y^{n + 1}

and v at

r^{*}

. If

r^{*}

is not a grid node, then interpolation is applied.

In order to establish the correctness of the numerical algorithm (29)–(32), we have to ensure that the recovered function

g_{L}^{n + 1}

and the solution

u_{j}^{n + 1}

n = 1, 2, \dots, N

j = 1, 2, \dots, J - 1

, exist and are unique. It is enough to prove that the systems (30) and (31) have a unique solution and

v_{r^{*}} \neq 0

Theorem 2.

The systems (30) and (31) have a unique solution,

0 < v_{j} \leq 1

j = 0, 1, \dots, J

, and, therefore, the algorithm (29)–(32) is correct.

Proof.

For the proof, we follow the same line of considerations as in [27].

Since the systems (30) and (31) are non-homogeneous with a coefficient matrix that is strictly diagonally dominant and has positive elements on the main diagonal (

Q_{j} > 0

P_{j} > 0

j = 0, 1, \dots, J

) and non-positive off diagonal elements, i.e., it is an M-matrix, we conclude that the solution is unique. Moreover, its inverse consists only of non-negative entries [46]. Taking into account that the right-hand side of (31) is a vector of non-negative elements, we deduce that

v_{j} \geq 0

j = 0, 1, \dots, J

Further, we consider the

J - 1

-th equation in (31)

Q_{J - 1} v_{J - 1} = P_{J - 1} v_{J - 2} + D_{J - 1} \geq D_{J - 1} \Rightarrow v_{J - 1} > \frac{D_{J - 1}}{Q_{J - 1}} > 0 .

Then, from the

J - 2

-th equation in (31), we obtain

v_{J - 2} \geq \frac{D_{J - 2} D_{J - 1}}{Q_{J - 2} Q_{J - 1}} > 0 .

We proceed in the same manner up to the first equation and derive

v_{j} > 0

j = 0, 1, \dots, J

Finally, according to the discrete maximum principle [42], the following estimate holds

max_{0 \leq j \leq J} | v_{j} | \leq 1 .

□

In the same manner, we deal with IP1 and IP4, if Dirichlet BC (4) at

r = 0

is imposed. We consider the decomposition

\begin{matrix} u_{j}^{n + 1} = y_{j}^{n + 1} + g_{0}^{n + 1} v_{j} \end{matrix}

(33)

and derive two problems: (30) and

\begin{array}{l} P_{0} v_{0} - D_{0} v_{1} = 1, \\ - P_{j} v_{j - 1} + Q_{j} v_{j} - D_{j} v_{j + 1} = 0, j = 1, 2, \dots, J - 1, \\ v_{J} = 0 . \end{array}

(34)

Computing (30) and (31) to determine y and v throughout the entire computational domain, and using (29) along with the measurement

u (r^{*}) = Ψ (t)

r^{*} \in Ω_{1} \cup Ω_{2}

, we recover the left boundary condition

g_{0}^{n + 1} = \frac{Ψ^{n + 1} - y_{r^{*}}^{n + 1}}{v_{r^{*}}^{n + 1}},

(35)

The correctness of (33)–(35) is proved as in Theorem 2.

Further, we proceed with numerically solving IP5 in the case of Dirichlet BC (4) at

r = 0

for a more general case of measurements

u (r_{i}^{*}) = Ψ_{i} (t), r_{i}^{*} \in Ω_{1} \cup Ω_{2}, i = 1, 2 .

(36)

Introducing the decomposition

\begin{matrix} u_{j}^{n + 1} = y_{j}^{n + 1} + g_{L}^{n + 1} v_{j} + g_{0}^{n + 1} w_{j} \end{matrix}

(37)

and substituting it into (28) yields three problems for the unknown functions

y_{j}^{n + 1}

v_{j}^{n + 1}

, and

w_{j}^{n + 1}

, respectively

\begin{array}{l} y_{0}^{n + 1} = 0, \\ - P_{j} y_{j - 1}^{n + 1} + Q_{j} y_{j}^{n + 1} - D_{j} y_{j + 1}^{n + 1} = F_{j}^{n}, j = 1, 2, \dots, J - 1, \\ y_{J}^{n + 1} = 0, \\ u_{j}^{0} = u^{0} (j), \end{array}

(38)

\begin{array}{l} v_{0} = 0, \\ - P_{j} v_{j - 1} + Q_{j} v_{j} - D_{j} v_{j + 1} = 0, j = 1, 2, \dots, J - 1, \\ v_{J} = 1 . \end{array}

(39)

\begin{array}{l} w_{0} = 1, \\ - P_{j} w_{j - 1} + Q_{j} w_{j} - D_{j} w_{j + 1} = 0, j = 1, 2, \dots, J - 1, \\ w_{J} = 0 . \end{array}

(40)

Solving these problems, using (39) and taking into account (36), we obtain

\begin{matrix} Ψ_{1}^{n + 1} = y_{r_{1}^{*}}^{n + 1} + g_{L}^{n + 1} v_{r_{1}^{*}} + g_{0}^{n + 1} w_{r_{1}^{*}}, \\ Ψ_{2}^{n + 1} = y_{r_{2}^{*}}^{n + 1} + g_{L}^{n + 1} v_{r_{2}^{*}} + g_{0}^{n + 1} w_{r_{2}^{*}} . \end{matrix}

(41)

From (41), we obtain

\begin{matrix} g_{0}^{n + 1} = \frac{v_{r_{1}^{*}} (Ψ_{2}^{n + 1} - y_{r_{2}^{*}}^{n + 1}) - v_{r_{2}^{*}} (Ψ_{1}^{n + 1} - y_{r_{1}^{*}}^{n + 1})}{v_{r_{1}^{*}} w_{r_{2}^{*}} - v_{r_{2}^{*}} w_{r_{1}^{*}}}, \\ g_{L}^{n + 1} = \frac{(Ψ_{1}^{n + 1} - y_{r_{1}^{*}}^{n + 1}) w_{r_{2}^{*}} - (Ψ_{2}^{n + 1} - y_{r_{2}^{*}}^{n + 1}) w_{r_{1}^{*}}}{v_{r_{1}^{*}} w_{r_{2}^{*}} - v_{r_{2}^{*}} w_{r_{1}^{*}}} . \end{matrix}

(42)

As before, the correctness (existence of a unique solution) of the inverse problem (37)–(42) can be established if the denominator in (42) is not zero and the systems (38)–(40) have unique solutions.

Theorem 3.

The solution of the systems (38)–(40) is unique,

v_{r_{1}^{*}} w_{r_{2}^{*}} - v_{r_{2}^{*}} w_{r_{1}^{*}} \neq 0

, and, therefore, the algorithm (37)–(42) is correct.

Proof.

Since the coefficient matrix of all systems above is one and the same, as in Theorem 2, we deduce that the solution of the systems (38)–(40) is unique and

v_{j} > 0

w_{j} > 0

j = 0, 1, \dots, J

Further, we multiply each equation in the system (39) by

w_{r_{2}^{*}}

and each equation in the system (40) by

v_{r_{2}^{*}}

and subtract the resulting systems to obtain

\begin{array}{l} w_{0} v_{r_{2}^{*}} - v_{0} w_{r_{2}^{*}} = v_{r_{2}^{*}}, \\ - P_{j} (v_{j - 1} w_{r_{2}^{*}} - w_{j - 1} v_{r_{2}^{*}}) + Q_{j} (v_{j} w_{r_{2}^{*}} - w_{j} v_{r_{2}^{*}}) \\ - D_{j} (v_{j + 1} w_{r_{2}^{*}} - w_{j + 1} v_{r_{2}^{*}}) = 0, j = 1, 2, \dots, J - 1, \\ v_{J} w_{r_{2}^{*}} - w_{J} v_{r_{2}^{*}} = w_{r_{2}^{*}} . \end{array}

(43)

Note that the coefficient matrix of (43) is the same M-matrix as for the systems (39) and (40), and the right-hand side is non-negative. Therefore,

w_{0} v_{r_{2}^{*}} - v_{0} w_{r_{2}^{*}} \geq 0

v_{j} w_{r_{2}^{*}} - w_{j} v_{r_{2}^{*}} \geq 0

j = 1, 2, \dots, J

As in Theorem 2, we consequently estimate the diagonal elements of (43), starting from the last one and moving to the first equation, and conclude that

w_{0} v_{r_{2}^{*}} - v_{0} w_{r_{2}^{*}} > 0

and

v_{j} w_{r_{2}^{*}} - w_{j} v_{r_{2}^{*}} > 0

j = 1, 2, \dots, J

Therefore, since

r_{1}^{*} \in {1, 2, \dots, J - 1}

, we complete the proof. □

6. Computational Examples

In this section, we verify the accuracy of the proposed numerical method. As a test example, we consider the problem (1)–(7) with parameters and functions

d_{1} = 2, d_{2} = 10, L = 3, l_{1} = 1, T = 1, g_{0} (t) = \frac{\sqrt{2}}{2} e^{t / 2}, g_{L} (t) = e^{t / 2} \frac{68 + \sqrt{2} (1 - π)}{2} .

We take

f_{1} (r, t)

and

f_{2} (r, t)

such that the exact solution of the problem (1)–(7) is

u_{1} (r, t) = e^{t / 2} cos (π x / 4)

u_{2} (r, t) = e^{t / 2} (a x^{4} - x^{3} + c)

, where

a = (120 - \sqrt{2} π) / 160

c = (40 + \sqrt{2} (80 + π)) / 160

Example 1

(Convergence: direct problem). We test the order of convergence of the direct problem. The numerical solution is computed using (28) for a fixed ratio

τ = h^{2}

. Errors and the convergence rate of the numerical solution

u_{i}^{n}

i = 1, 2

at the final time are estimated in the maximum and

L_{2}

norms. Let

E_{i, j} = u_{i} (r_{j}, T) - u_{i, j}^{M}

and

\begin{matrix} E_{1, \infty} = E_{1, \infty} (J) = max_{0 \leq j \leq j^{*}} | E_{1, j} |, E_{2, \infty} = E_{2, \infty} (J) = max_{j^{*} \leq J} | E_{2, j} |, \\ E_{1, L_{2}} = E_{1, L_{2}} (J) = {[\sum_{j = 0}^{j^{*}} h {(E_{1, j})}^{2}]}^{1 / 2}, E_{2, L_{2}} = E_{2, L_{2}} (J) = {[\sum_{j = j^{*}}^{J} h {(E_{2, j})}^{2}]}^{1 / 2}, \\ C R_{i, \infty} = {log}_{2} \frac{E_{i, \infty} (J)}{E_{i, \infty} (2 J)}, C R_{i, L_{2}} = {log}_{2} \frac{E_{i, L_{2}} (J)}{E_{i, L_{2}} (2 J)}, i = 1, 2, \end{matrix}

The results are listed in Table 1 and Table 2. In Table 1, we provide errors and convergence rates obtained from computations imposing Dirichlet BC at

r = 0

(4), while in Table 2, we present results in the case of natural BC at

r = 0

(16) and (17). We verify that, in both cases, the accuracy is

O (τ + h^{2})

Example 2

(Inverse problem: exact measurements). We numerically recover the functions

g_{0}

and

g_{L}

by using measurements equal to the numerical solution of the direct problem at

r_{1}^{*}

and

r_{2}^{*}

for exact external BCs, i.e.,

Ψ_{1}^{n} = u^{n} (r_{1}^{*})

and

Ψ_{2}^{n} = u^{n} (r_{2}^{*})

. The errors between the exact external BCs and the numerically restored ones (

g_{0}^{n}

g_{L}^{n}

) are denoted by

E_{g_{0}} = g_{0} (t_{n}) - g_{0}^{n}

E_{g_{L}} = g_{L} (t_{n}) - g_{L}^{n}

. We will also provide the error

E_{u}

between the recovered solution

u^{N} = (u_{1}^{N}, u_{2}^{N})

, obtained by numerically solving the inverse problems, and the numerical solution, obtained by the direct problem.

First, we examine the numerical method (29)–(32) in the case of natural BC at

r = 0

for recovering the function

g_{L}^{n}

and the solution

(u_{1}, u_{2})

upon measurement at the point

r_{1}^{*} = 0.5 \in Ω_{1}

. The computations are performed for

τ = h

and

J = 120

. In Figure 1 and Figure 2, we present the exact and recovered solution u and the function

g_{L}^{n}

, respectively, along with the corresponding errors.

Next, we present computational results obtained by solving (37)–(42),

τ = h

J = 120

r_{1}^{*} = 0.5

r_{2}^{*} = 2

for recovering

g_{0}

g_{L}

, and the solution

(u_{1}, u_{2})

. In Figure 3, Figure 4 and Figure 5, we present the exact and recovered solution u, the functions

g_{0}^{n}

g_{L}^{n}

, and the corresponding errors.

The numerical results show that the solution

(u_{1}, u_{2})

and the functions

(g_{0}, g_{L})

, obtained by numerically solving the inverse problems, exactly recover the solution of the direct problem and exact external boundary conditions. Therefore, the order of convergence is the same as that for the forward problem.

Example 3

(Inverse problems: noisy data). We perform the same experiments as in Example 2, but now we add perturbations to the measurements

Ψ_{i}^{n} = u^{n} (r_{i}^{*}) + 2 ρ_{i} (σ_{i}^{n} - 0.5) u^{n} (r_{i}^{*}), i = 1, 2,

where

u^{n} (r_{i}^{*})

are numerical solutions of the direct problem at time layer n and at spatial nodes

r_{i}^{*}

ρ_{i}

are levels of noise, and

σ_{i}^{n}

are random values, uniformly distributed on the interval

[0, 1]

. The observations are at

r_{1}^{*} = 0.5

and

r_{2}^{*} = 2

. We use polynomial curve fitting of degree 5 to smooth the data.

We estimate errors of the numerical solution

(u_{1}^{n}, u_{2}^{n})

, computed by inverse problem, in the maximum and

L_{2}

norms, using the same formulas as in Example 1, but

E_{i, j}

is replaced by

E_{u_{i}}

, as in Example 2. For the errors of the recovered functions

g_{0}

g_{L}

, we use formulas

ϵ_{s, \infty} = max_{0 \leq n \leq N} | E_{g_{s}} |, ϵ_{s, L_{2}} = {[\sum_{n = 0}^{N} τ E_{g_{s}}^{2}]}^{1 / 2}, s = {0, L} .

The computations are performed for

τ = h

and

J = 120

We consider the numerical method (29)–(32) in the case of natural BC at

r = 0

, for determining the function

g_{L}^{n}

and the solution

(u_{1}, u_{2})

under measurements at the point

r_{1}^{*}

. The computational results for different levels of noise

ρ = ρ_{1}

are presented in Table 3.

In Figure 6, we depict the exact and recovered function

g_{L}

and the corresponding error for

ρ = 0.1

. The largest error appears at the initial time. Note that for

ρ = 0.1

, the relative error in the maximum norm of the restored function

g_{L}

is 7.3208 × 10⁻², while in the

L_{2}

norm, it is 1.8183 × 10⁻², which is satisfactory for this level of noise.

Next, we perform experiments with the numerical method (37)–(42) for recovering the solution

(u_{1}, u_{2})

and the functions

g_{0}

and

g_{L}

with different levels of noise

ρ_{1}

and

ρ_{2}

. The computational results are listed in Table 4.

In Figure 7 and Figure 8, we plot the exact and numerically recovered functions

g_{0}

and

g_{L}

, respectively, and the corresponding errors for

ρ_{1} = 0.01

ρ_{3} = 0.003

We observe that for the noisy measurements, the numerical approach recovers the numerical solution and external boundary conditions with satisfactory accuracy.

The computational simulations presented in this section to test the developed methods are implemented using MATLAB^® R2022a.

7. Conclusions

In the present work, we develop an efficient numerical method for solving five inverse problems for recovering the solution and the right and/or left boundary conditions in a two-layer interface parabolic equation with cylindrical symmetry. We prove the well-posedness of the direct differential problem and construct a second-order in space and first-order in time numerical approach for its solution. Due to the degeneracy of the problem at the left boundary

r = 0

, we consider both the Dirichlet and natural boundary conditions. We construct a numerical method based on a decomposition technique for numerically solving the inverse problems. The correctness of the developed approaches is investigated.

Numerical experiments show the stability and the order of convergence—second in space and first in time—of the numerical method for solving the direct problem. Since for exact observations, numerical methods for solving inverse problems recover the numerical solution of the direct problem and external boundary conditions exactly, we conclude that the numerical algorithms for solving inverse problems are stable and achieve the same order of accuracy as the numerical method for the direct problem. Simulations with noisy observations illustrate that the recovered solution and external boundary conditions are of good accuracy.

Although in this study we consider continuous functions

g_{0} (t)

and

g_{L} (t)

, the constructed numerical methods can be easily extended to a wide class of functions, for example,

L_{2}

. In future work, we plan to extend our investigation to inverse spherical geometry problems.

Author Contributions

Conceptualization, 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, M.N.K. and L.G.V.; validation, M.N.K. All authors have read and agreed to the published version of the manuscript.

Funding

This study is financed by the European Union-NextGenerationEU, through the National Recovery and Resilience Plan of the Republic of Bulgaria, project № BG-RRP-2.013-0001-C01.

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 the paper.

Conflicts of Interest

The authors declare no conflicts of interest.

References

Abdollahi, A.N.; Rostamy, D. Identifying an unknown time-dependent boundary source in time-fractional diffusion equation with a non-local boundary condition. J. Comput. Appl. Math. 2019, 355, 3–50. [Google Scholar]
Alifanov, O.M. Inverse Heat Transfer Problems, 1st ed.; Springer: Berlin/Heidelberg, Germany, 1994; 348p. [Google Scholar]
Ammari, H.; Garnier, J.; Kang, H.; Nguen, L.; Seppecher, L. Multi-Wave Medical Imaging: Mathematical Modelling and Imaging Reconstruction; World Scientific: London, UK, 2017. [Google Scholar]
Badia, A.; Ha-Duong, T. On an inverse source problem for the heat equation. Application to a pollution detection problem. J. Inverse-Ill-Posed Probl. 2002, 10, 585–599. [Google Scholar]
Bogachev, I.V.; Nedin, R.D.; Vatulin, A.O.; Vavrugan, O.V. Identification of inhomogeneous elastic properties of isotropic cylinder. ZAMM Z. Angew. Math. Mech. 2017, 97, 358–364. [Google Scholar]
Cannon, J.R. Determination of an a unknown heat source from overspecified boundary data. SIAM J. Numer. Anal. 1968, 5, 275–286. [Google Scholar]
Cheng, W.; Liu, Y.-L.; Yang, F. A modified regularization method for a spherically symmetric inverse heat conduction problem. Symmetry 2022, 14, 2102. [Google Scholar] [CrossRef]
Demir, A.; Ozbilge, E. Identification of the unknown boundary condition in a linear parabolic equation. J. Inequal. Appl. 2013, 96, 1–7. [Google Scholar]
Duckstein, S.L.; Yakowitz, S. Instability in aquifer identification: Theory and case studies. Water Resour. Res. 1980, 16, 1045–1064. [Google Scholar]
Dzhonova-Atanasova, D.; Semkov, K.; Petrova, T.; Darakchiev, S.; Stefanova, K.; Nakov, S.; Popov, R. Liquid distribution in a semi-industrial packed column-experimental theory. Food Sci. Aplied Biotechnol. 2018, 1, 19–25. [Google Scholar]
Liu, K.; Zhao, Y. Inverse moving source problems for parabolic equations. Appl. Math. Lett. 2024, 155, 109114. [Google Scholar]
Marchuk, G.I. Mathematical modelling in environmental problems. Stud. Math. Its Appl. 1986, 16, 1986. [Google Scholar]
Okubo, A. Diffusion and Ecological Problems: Mathematical Models; Springer: Berlin/Heidelberg, Germany; New York, NY, USA, 1980. [Google Scholar]
Oulmelk, A.; Afraites, L.; Hadri, A.; Zaky, M.A.; Hendy, A.S. Alternating direction multiplier method to estimate an unknown source term in the time-fractional diffusion equation. Comput. Math. Appl. 2024, 156, 195–206. [Google Scholar]
Rundell, W.; Yin, H.-M. A parabolic inverse problem with an unknown boundary condition. J. Differ. Equ. 1990, 86, 234–242. [Google Scholar]
Hasanov, A.H.; Romanov, V.G. Introduction to Inverse Problems for Differential Equations, 1st ed.; Springer: Cham, Swizerland, 2017; 261p. [Google Scholar]
Isakov, V. Inverse Problems for Partial Differential Equations, 3rd ed.; Springer: Cham, Swizerland, 2017; p. 406. [Google Scholar]
Kabanikhin, S.I. Inverse and Ill-Posed Problems; DeGruyer: Berlin, Germany, 2011. [Google Scholar]
Lesnic, D. Inverse Problems with Applications in Science and Engineering; CRC Press: Abingdon, UK, 2021; p. 349. [Google Scholar]
Samarskii, A.A.; Vabishchevich, P.N. Numerical Methods for Solving Inverse Problems in Mathematical Physics; de Gruyter: Berlin, Germany, 2007; 438p. [Google Scholar]
Razzaghi, H.; Kowsary, F.; Layeghi, M. Inverse boundary problem in estimating heat transfer coefficient of a round pulsating bubbly jet: Design of experiment. Appl. Math. Sci. Eng. 2022, 30, 210–234. [Google Scholar]
Reddy, K.D.; Reddy, B.K. Estimation of inlet conditions of fluid flow in a thick pipe using inverse technique. Fluid Mech. Fluid Power 2024, 6, 11–19. [Google Scholar]
Wei, T.; Li, Y.S. An inverse boundary problem for one-dimensional heat equation with a multilayer domain. Eng. Anal. Bound. Elem. 2009, 33, 225–232. [Google Scholar]
Rundell, W.; Xub, X.; Zuo, L. The determination of an unknown boundary condition in a fractional diffusion equation. Appl. Anal. 2013, 92, 1511–1526. [Google Scholar]
Ozbilge, E. Determination of the unknown boundary condition of the inverse parabolic problems via semigroup method. Bound Value Probl. 2013, 2013, 2. [Google Scholar]
Vasil’ev, V.I.; Su, L. Numerical method for solving boundary inverse problem for one-dimensional parabolic equation. Math. Model. 2017, 24, 108–117. [Google Scholar]
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. [Google Scholar] [CrossRef]
Koleva, M.N.; Vulkov, L.G. Numerical reconstruction of time-dependent boundary conditions to 2D heat equation on disjoint rectangles at integral observations. Mathematics 2024, 12, 1499. [Google Scholar] [CrossRef]
Koleva, M.N.; Vulkov, L.G. Reconstruction of boundary conditions of a parabolic-hyperbolic transmission problem, In New Trends in the Applications of Differential Equations in Sciences; Slavova, A., Ed.; Springer Proceedings in Mathematics & Statistics; Springer: Cham, Swizerland, 2024; Volume 449, pp. 433–443. [Google Scholar]
Abidin, M.Z.; Marwan, M.; Ullah, N.; Mohamed Zidan, A. Well-posedness in variable-exponent function spaces for the three-dimensional micropolar fluid equations. J. Math. 2023, 1, 4083997. [Google Scholar]
Wu, J.; He, X.; Li, X. Finite-time stabilization of time-varying nonlinear systems based on a novel differential inequality approach. Appl. Math. Comput. 2022, 240, 126895. [Google Scholar]
Wang, P.; Wang, Y.; Jiang, C.; Li, T. Convergence of solutions for functional integro-differential equations with nonlinear boundary conditions. Adv. Differ. Equ. 2019, 521, 521. [Google Scholar]
Johnston, S.T.; Simpson, M.J. Exact solutions for diffusive transport on heterogeneous growing domains. arXiv 2023, arXiv:2304.09451. [Google Scholar]
Chernogorova, T.; Ewing, R.E.; Iliev, O.; Lazarov, R. On the discretization of interface problems with perfect and imperfect contact. In Conference Numerical Treatment of Multiphase Flows in Porous Media, Proceedings of the International Workshop, Beijing, China, 21–25 August 2000; Springer: Berlin/Heidelberg, Germany, 2000; Volume 2–6, pp. 93–103. [Google Scholar]
Ewing, R.; Iliev, O.; Lazarov, R.; Rybak, I.; Willems, J. A simplified method for upscaling composite materials with high contrast of the conductivity. SIAM J. Sci. Comput. 2009, 31, 2568–2586. [Google Scholar]
Lazarov, R.D.; Makarov, V.L. Difference scheme of second-order of accuracy for the axisymmetric Poisson equation in generalized solutions. J. USSR Comput. Math. Math. Phys. 1981, 21, 95–107. [Google Scholar]
Zimmerman, R.A.; Jankowski, T.A.; Tartakovsky, D.M. Analytical models of axisymmetric reaction–diffusion phenomena in composite media. Int. J. Heat Mass Transf. 2016, 99, 425–431. [Google Scholar]
Zhang, Z.; Zhang, Y.-X. A fast iterative method for identifying the radiogenic source for the helium production-diffusion equation. Appl. Math. Sci. Eng. 2022, 30, 521–540. [Google Scholar]
El-Gebeily, M.; Knuliev, Y.A. Identification of well-thinning and cracks in pipes utilizing vibration models and wavelets. Inverse Probl. Eng. Appl. Math. Model. 2016, 40, 5335–5348. [Google Scholar]
Yang, L.; Yu, J.-N.; Luo, G.-W.; Deng, Z.-C. Numerical identification of source terms for a two dimensional heat conduction problem in polar coordinate system. Appl. Math. Model. 2013, 37, 939–957. [Google Scholar]
Chegis, R.Y. A scheme of increased order of accuracy in the case of cylindrical and spherical symmetry. Comput. Math. Math. Phys. 1994, 34, 323–331. (In Russian) [Google Scholar]
Samarskii, A. The Theory of Difference Schems; Marcel Dekker: New York, NY, USA, 2001. [Google Scholar]
Samarskii, A.A. Tikhovov. A.N. Equations of Mathematical Physics; Courier Corporation: Chelmsford, MA, USA, 2013; 800p. [Google Scholar]
Ladyženskaja, O.A.; Solonnikov, V.A.; Ural’cev, N.N. Linear and Quasilinear Equations of Parabolic Type. Transl. Math. Monogr. 1968, 23, 648. [Google Scholar]
Andreev, V.K. On inequalities of the Friedrichs type for combined domains. J. Sib. Fed. Univ. Math. Phys. 2009, 2, 146–157. [Google Scholar]
Varga, R.S. Matrix Iterative Analysis, 2nd ed.; Springer: Berlin/Heidelberg, Germany, 2000; p. 358. [Google Scholar]

Figure 1. (Left): Numerical solution u of the direct problem (solid red line) and recovered solution (line with blue circles); (right): error between the numerically recovered solution and numerical solution of the direct problem,

τ = h

, IP solved by (29)–(32), Example 2.

τ = h

, IP solved by (29)–(32), Example 2.

Figure 2. (Left): Exact function

g_{L} (t)

(solid red line) and recovered function

g_{L}^{n}

(line with blue circles); (right): the corresponding error,

τ = h

, IP solved by (29)–(32), Example 2.

Figure 2. (Left): Exact function

g_{L} (t)

(solid red line) and recovered function

g_{L}^{n}

(line with blue circles); (right): the corresponding error,

τ = h

, IP solved by (29)–(32), Example 2.

Figure 3. (Left): Numerical solution u of the direct problem (solid red line) and recovered solution (line with blue circles); (right): error between the numerically recovered solution and numerical solution of the direct problem,

τ = h

, IP solved by (37)–(42), Example 2.

τ = h

, IP solved by (37)–(42), Example 2.

Figure 4. (Left): Exact function

g_{0} (t)

(solid red line) and recovered function

g_{0}^{n}

(line with blue circles); (right): the corresponding error,

τ = h

, IP solved by (37)–(42), Example 2.

Figure 4. (Left): Exact function

g_{0} (t)

(solid red line) and recovered function

g_{0}^{n}

(line with blue circles); (right): the corresponding error,

τ = h

, IP solved by (37)–(42), Example 2.

Figure 5. (Left): Exact function

g_{L} (t)

(solid red line) and recovered function

g_{L}^{n}

(line with blue circles); (right): the corresponding error,

τ = h

, IP solved by (37)–(42), Example 2.

Figure 5. (Left): Exact function

g_{L} (t)

(solid red line) and recovered function

g_{L}^{n}

(line with blue circles); (right): the corresponding error,

τ = h

, IP solved by (37)–(42), Example 2.

Figure 6. (Left): Exact function

g_{L} (t)

(solid red line) and recovered function

g_{L}^{n}

(line with blue circles); (right): the corresponding error,

ρ = 0.1

, IP solved by (37)–(42), Example 3.

Figure 6. (Left): Exact function

g_{L} (t)

(solid red line) and recovered function

g_{L}^{n}

(line with blue circles); (right): the corresponding error,

ρ = 0.1

, IP solved by (37)–(42), Example 3.

Figure 7. (Left): Exact function

g_{0} (t)

(solid red line) and recovered function

g_{0}^{n}

(line with blue circles); (right): the corresponding error,

ρ_{1} = 0.01

ρ_{3} = 0.003

, IP solved by (37)–(42), Example 3.

Figure 7. (Left): Exact function

g_{0} (t)

(solid red line) and recovered function

g_{0}^{n}

(line with blue circles); (right): the corresponding error,

ρ_{1} = 0.01

ρ_{3} = 0.003

, IP solved by (37)–(42), Example 3.

Figure 8. (Left): Exact function

g_{L} (t)

(solid red line) and recovered function

g_{L}^{n}

(line with blue circles); (right): the corresponding error,

ρ_{1} = 0.01

ρ_{3} = 0.003

, IP solved by (37)–(42), Example 3.

Figure 8. (Left): Exact function

g_{L} (t)

(solid red line) and recovered function

g_{L}^{n}

(line with blue circles); (right): the corresponding error,

ρ_{1} = 0.01

ρ_{3} = 0.003

, IP solved by (37)–(42), Example 3.

Table 1. Errors and convergence rates of the recovered solution

u_{i}^{n}

i = 1, 2

, Dirichlet BC, Example 1.

Table 1. Errors and convergence rates of the recovered solution

u_{i}^{n}

i = 1, 2

, Dirichlet BC, Example 1.

J	$E_{1, \infty}$	${CR}_{1, \infty}$	$E_{2, \infty}$	${CR}_{2, \infty}$	$E_{1, L_{2}}$	${CR}_{1, L_{2}}$	$E_{2, L_{2}}$	${CR}_{2, L_{2}}$
30	5.1356 × 10⁻²		5.1356 × 10⁻²		4.3835 × 10⁻²		4.8458 × 10⁻²
60	1.2944 × 10⁻²	1.9882	1.2944 × 10⁻²	1.9882	1.1176 × 10⁻²	1.9717	1.2016 × 10⁻²	2.0118
120	3.2557 × 10⁻³	1.9913	3.2556 × 10⁻³	1.9913	2.8523 × 10⁻³	1.9702	2.9957 × 10⁻³	2.0040
240	8.1772 × 10⁻⁴	1.9933	8.1772 × 10⁻⁴	1.9933	7.2696 × 10⁻⁴	1.9722	7.4882 × 10⁻⁴	2.0002
480	2.0520 × 10⁻⁴	1.9946	2.0520 × 10⁻⁴	1.9946	1.8488 × 10⁻⁴	1.9753	1.8739 × 10⁻⁴	1.9986
960	5.1458 × 10⁻⁵	1.9956	5.1457 × 10⁻⁵	1.9956	4.6914 × 10⁻⁵	1.9785	4.6913 × 10⁻⁵	1.9980

Table 2. Errors and convergence rates of the numerical solution

u_{i}^{n}

i = 1, 2

, natural BC, Example 1.

Table 2. Errors and convergence rates of the numerical solution

u_{i}^{n}

i = 1, 2

, natural BC, Example 1.

J	$E_{1, \infty}$	${CR}_{1, \infty}$	$E_{2, \infty}$	${CR}_{2, \infty}$	$E_{1, L_{2}}$	${CR}_{1, L_{2}}$	$E_{2, L_{2}}$	${CR}_{2, L_{2}}$
30	6.2045 × 10⁻²		5.4563 × 10⁻²		6.0694 × 10⁻²		5.0790 × 10⁻²
60	1.5682 × 10⁻²	1.9842	1.3650 × 10⁻²	1.9990	1.4837 × 10⁻²	2.0324	1.2518 × 10⁻²	2.0205
120	3.9612 × 10⁻³	1.9851	3.4131 × 10⁻³	1.9998	3.6648 × 10⁻³	2.0174	3.1063 × 10⁻³	2.0107
240	1.0004 × 10⁻³	1.9853	8.5331 × 10⁻⁴	1.9999	9.1049 × 10⁻⁴	2.0090	7.7366 × 10⁻⁴	2.0054
480	2.5264 × 10⁻⁴	1.9855	2.1333 × 10⁻⁴	2.0000	2.2689 × 10⁻⁴	2.0046	1.9305 × 10⁻⁴	2.0027
960	6.3794 × 10⁻⁵	1.9856	5.3333 × 10⁻⁵	2.0000	5.6630 × 10⁻⁵	2.0024	4.8216 × 10⁻⁵	2.0014

Table 3. Errors of the recovered solution

(u_{1}^{N}, u_{2}^{N})

and function

g_{L}

for different levels of noise

ρ

, IP solved by (29)–(32), natural BC at

r = 0

, Example 3.

Table 3. Errors of the recovered solution

(u_{1}^{N}, u_{2}^{N})

and function

g_{L}

for different levels of noise

ρ

, IP solved by (29)–(32), natural BC at

r = 0

, Example 3.

$ρ$	$E_{1, \infty}$	$E_{2, \infty}$	$E_{1, L_{2}}$	$E_{2, L_{2}}$	$ϵ_{L, \infty}$	$ϵ_{L, 2}$
0.01	2.2141 × 10⁻³	1.0083 × 10⁻²	1.1401 × 10⁻³	7.7277 × 10⁻³	1.8319 × 10⁻¹	4.3869 × 10⁻²
0.03	6.1225 × 10⁻³	2.6653 × 10⁻²	3.2143 × 10⁻³	2.0666 × 10⁻²	6.6437 × 10⁻¹	1.6369 × 10⁻¹
0.05	1.0031 × 10⁻²	4.3224 × 10⁻²	5.2888 × 10⁻³	3.3605 × 10⁻²	1.1456	2.8357 × 10⁻¹
0.1	1.9802 × 10⁻²	8.4650 × 10⁻²	1.0475 × 10⁻²	6.5954 × 10⁻²	2.3486	5.8330 × 10⁻¹

Table 4. Errors of the recovered solution

(u_{1}^{N}, u_{2}^{N})

and functions

g_{0}

and

g_{L}

for different levels of noise

ρ_{1}

ρ_{2}

, IP solved by (37)–(42), Example 3.

Table 4. Errors of the recovered solution

(u_{1}^{N}, u_{2}^{N})

and functions

g_{0}

and

g_{L}

for different levels of noise

ρ_{1}

ρ_{2}

, IP solved by (37)–(42), Example 3.

$ρ_{1}$	$ρ_{2}$	$E_{1, \infty}$	$E_{2, \infty}$	$E_{1, L_{2}}$	$E_{2, L_{2}}$	$ϵ_{0, \infty}$	$ϵ_{0, 2}$	$ϵ_{L, \infty}$	$ϵ_{L, 2}$
0.03	0.001	1.051 × 10⁻¹	4.042 × 10⁻³	2.287 × 10⁻²	2.770 × 10⁻³	1.768 × 10⁻¹	3.826 × 10⁻²	7.408 × 10⁻³	1.645 × 10⁻³
0.01	0.003	2.763 × 10⁻²	1.239 × 10⁻²	6.186 × 10⁻³	8.439 × 10⁻³	6.863 × 10⁻²	1.404 × 10⁻²	2.107 × 10⁻²	4.680 × 10⁻³
0.05	0.005	1.630 × 10⁻¹	2.081 × 10⁻¹	3.568 × 10⁻²	1.417 × 10⁻²	2.980 × 10⁻¹	6.130 × 10⁻²	3.490 × 10⁻²	7.833 × 10⁻³

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.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Koleva, M.N.; Vulkov, L.G. Inverse Boundary Conditions Interface Problems for the Heat Equation with Cylindrical Symmetry. Symmetry 2024, 16, 1065. https://doi.org/10.3390/sym16081065

AMA Style

Koleva MN, Vulkov LG. Inverse Boundary Conditions Interface Problems for the Heat Equation with Cylindrical Symmetry. Symmetry. 2024; 16(8):1065. https://doi.org/10.3390/sym16081065

Chicago/Turabian Style

Koleva, Miglena N., and Lubin G. Vulkov. 2024. "Inverse Boundary Conditions Interface Problems for the Heat Equation with Cylindrical Symmetry" Symmetry 16, no. 8: 1065. https://doi.org/10.3390/sym16081065

APA Style

Koleva, M. N., & Vulkov, L. G. (2024). Inverse Boundary Conditions Interface Problems for the Heat Equation with Cylindrical Symmetry. Symmetry, 16(8), 1065. https://doi.org/10.3390/sym16081065

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Inverse Boundary Conditions Interface Problems for the Heat Equation with Cylindrical Symmetry

Abstract

1. Introduction

2. Direct and Inverse Problems

2.1. Direct (Forward) Problem

2.2. Inverse Problems

3. Well-Posedness of the Direct Problem

4. Numerical Solution of the Direct Problem

5. Numerical Solution of the Inverse Problems

6. Computational Examples

7. Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI