Identification of Boundary Conditions in a Spherical Heat Conduction Transmission Problem

Abstract

Although numerous analytical and numerical methods have been developed for inverse heat conduction problems in single-layer materials, few methods address such problems in composite materials. The following paper studies inverse interface problems with unknown boundary conditions by using interior point observations for heat equations with spherical symmetry. The zero degeneracy at the left interval $0<r<R_1$ leads to solution difficulties in the one-dimensional interface problem. So, we first investigate the well-posedness of the direct (forward) problem in special weighted Sobolev spaces. Then, we formulate three groups of unknown boundary conditions and inverse problems upon internal point measurements for the heat equation with spherical symmetry. Second-order finite difference scheme approaches for direct and inverse problems are developed. Computational test examples illustrate the theoretical statements proposed.

1. Introduction

The inverse unknown boundary condition (BC) problems for parabolic equations have many applications, such as heat–mass transfer, migration of groundwater, biology, identification and control of pollution, and environmental protection, see, e.g., [1,2,3,4]. As inverse ill-posed problems, they are of great interest in both mathematical and numerical studies [5,6,7,8].

The considered abstract problem can be physically interpreted as heat conduction within a spherical body composed of two layers, each having a different thermal conductivity. In this scenario, the time-dependent heat sources at both external boundaries are unknown, while temperature measurements are taken at internal points within the body. Convective heat transfer is neglected, as the primary factor influencing the process is its spherical geometry. The importance of this model has been explored in numerous books and research papers; see, for example, [5,9,10,11,12,13,14,15,16,17,18,19,20].

Problems involving partial differential equations (PDEs) in spherical coordinates have been extensively studied in the literature. A variety of analytical, semi-analytical, and numerical methods have been developed to solve such equations effectively.

Analytical solutions similar to the 1D heat equation in spherical coordinates are constructed in [20]. This technique is applied for solving finite line source and Stefan problems. In [13], the authors consider heat conduction equations in cylindrical and spherical geometry. They develop new analytical solutions, derived through a self-similar Ansatz. Subsequently, these solutions are accurately replicated using explicit and finite difference methods, which are unconditionally stable. In [20], the authors construct a general similarity-type solution for the 1D heat equation in spherical symmetry, and the solution is represented using Kummer functions. Semi-analytical solutions for axisymmetric reaction–diffusion equations involving first-order reaction kinetics and time-dependent boundary conditions are constructed in [21].

The authors of [9] developed a high-order numerical scheme for solving diffusion equations in cylindrical and spherical symmetry. They studied the stability, taking into account the singularity at the left boundary. Direct solvers, based on the fast Fourier transform for the Poisson equation on two-dimensional polar and spherical coordinates, were developed in [15]. A novel mesh-free numerical method for solving Stefan-type problems in spherical coordinates was developed in [22]. A finite difference scheme for 3D incompressible flows in spherical geometry is proposed in [23].

For PDEs, problems with layers or interface problems are considered in many papers; see, e.g., [24,25,26,27,28,29]. In [24,25], the authors construct semi-analytical solutions for multilayer advection–dispersion and dispersion problems. Exact solutions for diffusive transfer models on domains that exhibit spatially varied growth are proposed in [28]. The finite volume method for solving time-dependent multilayer diffusion problems is developed in [29]. The method can be applied to a variety of boundary/interface conditions. In [16,26], authors construct and analyze finite difference schemes for solving interface elliptic systems and the axisymmetric Poisson equation. The high contrast conductivity problem is theoretically and numerically investigated in [27].

In recent years, numerous studies have been devoted to developing stable and efficient methods for solving inverse problems, reflecting the increasing focus on addressing the challenges they present. Effective solving of inverse boundary identification problems for PDEs is important, as these conditions significantly influence the behavior of the solution.

The uniqueness of solutions to inverse problems for recovering boundary conditions in parabolic and time-fractional diffusion equations based on temperature observations at a fixed boundary point was explored in [30,31]. The efficient meshless method was developed in [32] to address inverse problems for recovering boundary sources in time-fractional diffusion equations. The authors of [33,34] utilized a semigroup method to investigate inverse problems related to determining the correct boundary condition in linear and quasi-linear parabolic equations, given the overspecified Dirichlet and Neumann data on the left boundary. Decomposition of the solution was applied in [35] to solve the numerical Dirichlet boundary condition inverse problem in the heat flow equation. An inverse problem for determining Dirichlet BC in the heat equation under interior measurements was developed in [36]. A deep neural network was applied in [37] to recover unknown thermo-physical parameters, such as the material’s thermal diffusivity and boundary conditions (heat flux) in a heat transfer problem.

Research on inverse problems for PDEs with interfaces or in spherical geometries is less prevalent compared to other domains. These types of problems present additional complexities due to the discontinuities at interfaces or the unique challenges posed by spherical coordinates. As a result, fewer studies focus on developing robust methods for these cases.

Numerical methods for solving boundary identification inverse problems for parabolic, parabolic-hyperbolic, and time-fractional parabolic problems on disjoint domains with interface conditions are constructed and investigated in [38,39,40]. The authors of [41] consider a one-dimensional heat equation with a multilayer domain. They propose a regularization method for solving inverse problems for recovering a moving boundary from Cauchy data.

The papers [10,11,12,14,19,42] focus on solving the inverse problem with spherical symmetry.

In [10,11], an inverse heat conduction problem with spherical geometry is studied, aimed at restoring the temperature distribution on the inner surface of a hollow sphere based on observations at a fixed position within the sphere. The authors develop a modified quasi-boundary value method to solve this problem. The regularization method of fundamental solutions is utilized in [14] to determine the data on the inner boundary from measurements on the outer boundary in radially symmetric heat problems. An iterative regularization approach is proposed in [12] to recover the source in the heat time-dependent diffusion equation on a domain with spherical symmetry. In [42], an inverse problem is investigated to determine a spatially varying source in the heat equation within a polar coordinate system using finite temperature data. Based on an eigenfunction expansion, an analytical solution for the inverse heat conduction source problem in a spherical diffusion domain is derived in [19]. Numerical regularization methods are also developed and analyzed.

In this work, we develop a numerical approach for solving boundary condition inverse problems for recovering external Dirichlet BCs on the basis of interior point observations for interface heat equations with spherical symmetry. This is a continuation of research the authors conducted in [43], where inverse boundary conditions parabolic problems with cylindrical symmetry are numerically solved.

The novelty of the present paper is proof of the well-posedness of the direct problem and the construction of the finite volume difference scheme for solving the spherical interface problem, taking into account degeneration. Also, we propose a numerical decomposition method for solving external boundary condition inverse interface problems in spherical symmetry.

The outline of this paper is as follows. In the next section, we formulate direct and inverse problems. In Section 3, we study the existence of the solution to the direct problem. A numerical approach for solving direct problems is developed in Section 4. In Section 5, we construct a computational approach for solving boundary identification inverse problems. Section 6 is devoted to the numerical experiments and validation of the proposed methods. The paper is finalized with concluding remarks.

2. Statement of the Direct and Inverse Problems

In this section, we discuss direct and inverse model problems.

2.1. Formulation of the Direct (Forward) Problems

We consider the interface initial-boundary value problem as follows: the parabolic equations

$$\begin{array}{ccc}& & {\displaystyle \frac{\partial u_1}{\partial t}}-d_1\left({\displaystyle \frac{{\partial }^2{u}_1}{\partial r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\partial u_1}{\partial r}}\right)+f_1(r,t),r\in {\mathsf{\Omega }}_1=(0,R_1),0<t<T,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & {\displaystyle \frac{{\partial }^2{u}_2}{\partial t^2}}-d_2\left({\displaystyle \frac{{\partial }^2{u}_1}{\partial r^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\partial u_2}{\partial r}}\right)+f_2(r,t),r\in {\mathsf{\Omega }}_2=(R_1,R_2),0<t<T,\hfill \end{array}$$

(2)

are subjected to initial conditions

$$\begin{array}{ccc}& & u_1(r,0)=u_{10}(r),0<r<R_1,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & u_2(r,0)=u_{20}(r),R_1<r<R_2,\hfill \end{array}$$

(4)

external boundary conditions

$$u_1(0,t)=g_0(t),u_2(R_2,t)=g_R(t),0<t<T,$$

(5)

and interface (internal boundary) conditions

$$\begin{array}{ccc}& & u_1(R_1,t)=u_2(R_1,t),0<t<T,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}& & D_1{\displaystyle \frac{\partial u_1}{\partial r}}(R_1,t)=D_2{\displaystyle \frac{\partial u_2}{\partial r}}(R_1,t).\hfill \end{array}$$

(7)

In Equations (1), (2), and (7), $d_1,d_2$ and $D_1,D_2$ are positive constants. In the case of $D_i=d_i$, $i=1,2$, we have continuity of the flux.

In (1)–(7), if all constants are given and functions are known, the direct (forward) problem is to find the function $u=(u_1(r,t),u_2(r,t))$.

2.2. Formulation of Inverse Problems

We study inverse problems for identifying unknown functions $g_0(t)$ and $g_R(t)$ on the basis of internal point observations. Depending on the measured data, the inverse problems (IP) can be formulated as follows:

• IP1a: Using observed data $u_1(r_1^{*},t)={\psi }_1(t)$, $r_1^{*}\in {\mathsf{\Omega }}_1$, we determine $(u_1,u_2,g_0(t))$;
• IP1b: Using observed data $u_2(r_2^{*},t)={\psi }_2(t)$, $r_2^{*}\in {\mathsf{\Omega }}_2$, we determine $(u_1,u_2,g_0(t))$;
• IP2a: Using observed data $u_2(r_2^{*},t)={\psi }_2(t)$, $r_2^{*}\in {\mathsf{\Omega }}_2$, we determine $(u_1,u_2,g_R(t))$;
• IP2b: Using observed data $u_1(r_1^{*},t)={\psi }_1(t)$, $r_1^{*}\in {\mathsf{\Omega }}_1$, we determine $(u_1,u_2,g_R(t))$;
• IP3a: Using observed data $u_i(r_i^{*},t)={\psi }_i(t)$, $r_i^{*}\in {\mathsf{\Omega }}_i$, $i=1,2$, we determine $(u,g_0(t),g_R(t))$;
• IP3b: Using observed data $u_1(r_i^{*},t)={\psi }_i(t)$, $r_i^{*}\in {\mathsf{\Omega }}_1$, $i=1,2$, we determine $(u,g_0(t),g_R(t))$;
• IP3c: Using observed data $u_2(r_i^{*},t)={\psi }_i(t)$, $r_i^{*}\in {\mathsf{\Omega }}_2$, $i=1,2$, we determine $(u,g_0(t),g_R(t))$.

3. Analysis of the Solution of the Direct Problem

An analytical study for the direct problem is carried out in this section to give rigorous proof of the well-posedness of the direct problem (1)–(7). The main difficulties that arise in solving the direct problem are connected to the singularity at $r=0$ and the interface.

The main result in this section is the following assertion.

Theorem 1.

Assume that $\underset{r\to 0}{lim}{u}_1(r,t)$ is bounded, $g_i\in C(0,T)$, $f_i\in Q_{iT}={\mathsf{\Omega }}_i\times (0,T]$, $u_{i0}(x)\in C({\mathsf{\Omega }}_i)$, $i=1,2$. Then, there exists a unique solution $u=\{u_1,u_2\}$, $u_1\in C^2({\mathsf{\Omega }}_1)\cap C^1(\{r|r=R_1\})$, $u_2\in C^2({\mathsf{\Omega }}_2)\cap C^1(\{r|r=R_1\})\cap C^1(\{r|r=L\})$ to the problem (1)–(7) that continuously depends on the input data.

Proof. 

Without loss generality, we study the particular case of zero Dirichlet BCs (4)–(5):

$$u_1(0,t)=g_0(t)=0,u_2(L,t)=g_L(t)=0.$$

(8)

We consider the uniqueness and continuous dependence on input data of the solution, assuming existence.

We rewrite the system (1)–(2) in conservative form and multiply Equation (1) by $r^2{u}_1(r,t)$ and Equation (2) by $r^2{u}_2(r,t)$. Next, we integrate the new equations over the interval ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$, respectively, and sum up the results to obtain the following:

$${\displaystyle \frac{\partial }{\partial t}}U(t)+d_1\underset{0}{\overset{R_1}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_1}{\partial r}}\right)}^{2}dr+d_2\underset{R_1}{\overset{R_2}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_2}{\partial r}}\right)}^{2}dr=\underset{0}{\overset{R_1}{\int }}{r}^2{u}_1{f}_{1}dr+\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_2{f}_{2}dr,$$

(9)

where

$$U(t)={\displaystyle \frac{1}{2}}\underset{0}{\overset{R_1}{\int }}{r}^2{u}_1^2(r,t)dr+{\displaystyle \frac{1}{2}}\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_2^2(r,t)dr.$$

(10)

Let us note that for the function $u_2(r,t)$, it is easy to derive a Friedrich-type inequality; see, e.g., [43,44]

$$\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_2^{2}dr\le C\underset{R_1}{\overset{R_2}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_2}{\partial r}}\right)}^{2}dr,C={\displaystyle \frac{R_1^2}{4}}\left({\displaystyle \frac{R_2^2}{R_1^2}}-2ln{\displaystyle \frac{R_2}{R_1}}-1\right)>0.$$

Due to the the singularity at $r=0$ in (1), the Friedrich-type inequality has the following form; see, e.g., [43]:

$$\underset{0}{\overset{R_1}{\int }}{r}^2{u}_1^{2}dr+\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_2^{2}dr\le C\left(d_1\underset{0}{\overset{R_1}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_1}{\partial r}}\right)}^{2}dr+d_2\underset{R_1}{\overset{R_2}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_2}{\partial r}}\right)}^{2}dr\right),$$

(11)

where $C>0$ is a constant, independent of the functions $u_1,u_2$.

Now, using the last two inequalities, we deduce that the left-hand side of equality (9) is greater or equal to $t_0$:

$${\displaystyle \frac{\partial U}{\partial t}}\le {\displaystyle \frac{1}{C}}\left(\underset{0}{\overset{R_1}{\int }}{r}^2{u}_1^{2}dr+\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_2^{2}dr\right).$$

For the right-hand side of (9), again using (11), we get the following:

$$\begin{array}{ccc}& & \underset{0}{\overset{R_1}{\int }}{r}^2{u}_1{f}_{1}dr+\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_2{f}_{2}dr\le {\left(\underset{0}{\overset{R_1}{\int }}{r}^2{f}_1^{2}dr\right)}^{1/2}{\left(\underset{0}{\overset{R_1}{\int }}{r}^2{u}_1^{2}dr\right)}^{1/2}\hfill \\ & & +{\left(\underset{R_1}{\overset{R_2}{\int }}{r}^2{f}_2^{2}dr\right)}^{1/2}{\left(\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_2^{2}dr\right)}^{1/2}.\hfill \\ & & max(R_1.{\widehat{f}}_1,\sqrt{R_2^2-R_1^2}{\widehat{f}}_2)\left[{\left({\displaystyle \frac{1}{2}}\underset{0}{\overset{R_1}{\int }}{r}^2{u}_1^{2}dr\right)}^{1/2}+{\displaystyle \frac{1}{2}}{\left(\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_2^{2}dr\right)}^{1/2}\right)\hfill \\ & & \le max(R_1{\widehat{f}}_1,\sqrt{R_2^2-R_1^2}{\widehat{f}}_2)(2\sqrt{2}\sqrt{U})=F(t)\sqrt{U},\hfill \end{array}$$

(12)

where

$${\widehat{f}}_i(t)=\sqrt{U\underset{{\mathsf{\Omega }}_i}{max}\left|f_1(x,t)\right|},i=1,2.$$

Therefore, from (9) and (12), we derive the inequalities

$${\displaystyle \frac{d}{dt}}U+cU\le F(t)U.c={\displaystyle \frac{1}{C}}min(d_1,d_2),$$

which implies

$$U(t)\le {\left(\sqrt{U_0}+{\displaystyle \frac{1}{2}}\underset{0}{\overset{t}{\int }}F(\tau )e^{c\tau }d\tau \right)}^2{e}^{-2ct},$$

(13)

where

$$U_0=U(0)={\displaystyle \frac{1}{2}}\underset{0}{\overset{R_1}{\int }}{r}^2{u}_{10}^2(r)dr+{\displaystyle \frac{1}{2}}\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_{20}^2(r)dr.$$

Therefore, from (10), we have the following:

$$\underset{0}{\overset{R_1}{\int }}{r}^2{u}_1^2(r,t)dr\le 2U(t),\underset{R_1}{\overset{R_2}{\int }}{r}^2{u}_2^2(r,t)dr\le 2U(t),$$

where $U(t)$ satisfies (13).

Next, we will prove the boundedness of the integrals

$$\underset{0}{\overset{R_1}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_1}{\partial r}}\right)}^{2}dr\mathrm{and}C\underset{R_1}{\overset{R_2}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_2}{\partial r}}\right)}^{2}dr.$$

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

$$\begin{array}{ccc}& & \underset{0}{\overset{t}{\int }}\underset{0}{\overset{R_1}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_1}{\partial t}}\right)}^{2}drds-2d_1\underset{0}{\overset{t}{\int }}\underset{0}{\overset{R_1}{\int }}{\displaystyle \frac{\partial u_1}{\partial t}}{\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial u_1}{\partial r}}\right)drds\hfill \\ & & +\underset{0}{\overset{t}{\int }}\underset{0}{\overset{R_1}{\int }}{d}_1^2.{\displaystyle \frac{1}{r^2}}{\left({\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial u_1}{\partial r}}\right)\right)}^{2}drds=\underset{0}{\overset{t}{\int }}\underset{0}{\overset{R_1}{\int }}{r}^2.f_1^{2}drds.\hfill \\ & & \underset{0}{\overset{t}{\int }}\underset{R_1}{\overset{R_2}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_2}{\partial t}}\right)}^{2}drds-2d_2\underset{0}{\overset{t}{\int }}\underset{R_1}{\overset{R_2}{\int }}{\displaystyle \frac{\partial u_2}{\partial t}}{\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial u_2}{\partial r}}\right)drds\hfill \\ & & +\underset{0}{\overset{t}{\int }}\underset{R_1}{\overset{R_2}{\int }}{d}_2^2{\displaystyle \frac{1}{r^2}}{\left({\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial u_2}{\partial r}}\right)\right)}^{2}drds.\hfill \end{array}$$

Summing these equalities, using integration by parts and the initial, boundary, and interface conditions, we obtain the following equality:

$$\begin{array}{ccc}& & {\displaystyle \frac{D_1^2}{d_1^2}}\underset{0}{\overset{t}{\int }}\underset{R_1}{\overset{R_2}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_1}{\partial t}}\right)}^{2}drds+{\displaystyle \frac{D_2^2}{d_1^2}}\underset{0}{\overset{t}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_2}{\partial t}}\right)}^{2}drds+{\displaystyle \frac{D_1^2}{d_1^2}}\underset{0}{\overset{R_1}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_1}{\partial r}}\right)}^{2}dr\hfill \\ & & +{\displaystyle \frac{D^2}{d_1^2}}\underset{R_1}{\overset{R_2}{\int }}{r}^2{\left({\displaystyle \frac{\partial u_2}{\partial r}}\right)}^{2}dr+\underset{0}{\overset{t}{\int }}{D}_1^2{\displaystyle \frac{1}{r^2}}{\left({\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial u_1}{\partial r}}\right)\right)}^{2}drds\hfill \\ & & +\underset{0}{\overset{t}{\int }}{D}_2^2{\displaystyle \frac{1}{r^2}}{\left({\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial u_2}{\partial r}}\right)\right)}^2=\underset{0}{\overset{t}{\int }}\underset{0}{\overset{R_1}{\int }}{r}^2{f}_1^{2}drds+\underset{0}{\overset{t}{\int }}\underset{R_1}{\overset{R_2}{\int }}{r}^2{f}_2^{2}drds.\hfill \end{array}$$

Therefore,

$$\underset{{\mathsf{\Omega }}_i}{\int }{r}^2{\left({\displaystyle \frac{\partial u_1}{\partial r}}\right)}^{2}dr\le {\displaystyle \frac{D_1}{d_1}}\underset{0}{\overset{t}{\int }}\underset{R_1}{\overset{R_2}{\int }}{r}^2{f}_1^{2}drds+{\displaystyle \frac{D_2}{d_2}}\underset{0}{\overset{t}{\int }}\underset{R_1}{\overset{R_2}{\int }}{r}^2{f}_2^{2}drds.$$

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

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

4. Finite Difference Solution of the Direct Problem

We start with the numerical discretization of the direct problem. We represent Equations (1) and (2) in divergent form:

$$\begin{array}{ccc}& & {\displaystyle \frac{\partial u_1}{\partial t}}-{\displaystyle \frac{d_1}{r^2}}{\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial u_1}{\partial r}}\right)=f_1(r,t),r\in {\mathsf{\Omega }}_1,0<t\le T,\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & {\displaystyle \frac{\partial u_2}{\partial t}}-{\displaystyle \frac{d_2}{r^2}}{\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial u_2}{\partial r}}\right)=f_2(r,t),r\in {\mathsf{\Omega }}_2,0<t\le T.\hfill \end{array}$$

(15)

At $r=0$, the solution degenerates and according to [46],

$$\underset{r\to 0}{lim}{r}^2{\displaystyle \frac{\partial u_1}{\partial r}}=0,$$

(16)

which ensures the solution $u_1(r,t)$ is bounded to

$$\left|u_1(0)\right|<\infty .$$

(17)

Consequently, taking into account (17), at $r=0$, we may impose either the Dirichlet boundary condition (5) or the natural boundary condition (16).

Consider the spatial and temporal meshes with uniform steps:

$${\overline{\omega }}_h=\{r_j=jh,j=0,1,\dots ,J,Jh=R_2\},{\overline{\omega }}_{\tau }=\{t_n=n\tau ,n=0,1,\dots ,N,N\tau =T\}.$$

The solution $u_i$, $i=1,2$ at the mesh point $(r_j,t_n)\in {\overline{\omega }}_h\times {\overline{\omega }}_{\tau }$ is denoted by $u_{i,j}^n$, while by $u_{i,j+1/2}^n$, we denote the solution $u_i$ at the mesh point $(r_{j+1/2},t_n)$, $r_{j+1/2}=r_j+h/2$, where $r_{j+1/2}$ are nodes of the dual spatial grid $0=r_{-1/2}<r_{1/2}<r_{3/2}\cdots <r_{J-1/2}<r_{J+1/2}=L$.

To develop a conservative scheme and address degeneracy at $r=0$, we adopt the method from [46] and implement the finite volume approach for systems (3)–(5), (14) and (15), taking into account (16). Let $R_1$ be a grid node, where $R_1=r_{j^{*}}=j^{*}h$.

At the internal spatial points $j\ne j^{*}$, we multiply (14) and (15) by $r^2$, integrate over the volumes first for $(r_{j-1/2},r_{j+1/2})$, $j=1,2,\dots j^{*}-1$, then over $(r_{j-1/2},r_{j+1/2})$, $j=j^{*}+1,\dots ,J-1$, and divide the resulting expression by $h{r}_i^2$ to derive the following:

$$\begin{array}{lll}& & {\displaystyle \frac{1}{r_i^{2}h}}\underset{r_{j-1/2}}{\overset{r_{j+1/2}}{\int }}{r}^2{\displaystyle \frac{\partial u_1}{\partial t}}dr-{\displaystyle \frac{d_1}{r_j^{2}h}}\left(r_{j+1/2}^2{\displaystyle \frac{\partial u_1}{\partial r}}{|}_{r_{j+1/2}}-r_{j-1/2}^2{\displaystyle \frac{\partial u_1}{\partial r}}{|}_{r_{j-1/2}}\right)\hfill \\ & & ={\displaystyle \frac{1}{r_j^{2}h}}\underset{r_{j-1/2}}{\overset{r_{j+1/2}}{\int }}{f}_1(r,t)r^{2}dr,j=1,2,\dots ,j^{*}-1,\hfill \end{array}$$

(18)

$$\begin{array}{lll}& & {\displaystyle \frac{1}{r_j^{2}h}}\underset{r_{j-1/2}}{\overset{r_{j+1/2}}{\int }}{r}^2{\displaystyle \frac{\partial u_2}{\partial t}}dr-{\displaystyle \frac{d_2}{r_j^{2}h}}\left(r_{j+1/2}^2{\displaystyle \frac{\partial u_2}{\partial r}}{|}_{r_{j+1/2}}-r_{j-1/2}^2{\displaystyle \frac{\partial u_2}{\partial r}}{|}_{r_{j-1/2}}\right)\hfill \\ & & ={\displaystyle \frac{1}{r_j^{2}h}}\underset{r_{j-1/2}}{\overset{r_{j+1/2}}{\int }}{f}_2(r,t)r^{2}dr,j=j^{*}+1,\dots ,J-1.\hfill \end{array}$$

(19)

To obtain discretization at the interface node $j=j^{*}$, we use an approach analogous to that in [29]. Multiplying (14) and (15) by $r^2$, integrating over volumes $(r_{j-1/2},r_j)$ and $(r_j,r_{j+1/2})$, and dividing the resulting expressions by $h{r}_i^2$, we get the following:

$$\begin{array}{ccc}& & {\displaystyle \frac{1}{r_j^{2}h}}\underset{r_{j-1/2}}{\overset{r_j}{\int }}{r}^2{\displaystyle \frac{\partial u_1}{\partial t}}dr-{\displaystyle \frac{d_1}{r_j^{2}h}}\left(r_j^2{\displaystyle \frac{\partial u_1}{\partial r}}{|}_{r_j}-r_{j-1/2}^2{\displaystyle \frac{\partial u_1}{\partial r}}{|}_{r_{j-1/2}}\right)\hfill \\ & & ={\displaystyle \frac{1}{r_j^{2}h}}\underset{r_{j-1/2}}{\overset{r_j}{\int }}{f}_1(r,t)r^{2}dr,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& & {\displaystyle \frac{1}{r_j^{2}h}}{\int }_{r_j}^{r_{j+1/2}}{r}^2{\displaystyle \frac{\partial u_2}{\partial t}}dr-{\displaystyle \frac{d_2}{r_j^{2}h}}\left(r_j^2{\displaystyle \frac{\partial u_2}{\partial r}}{|}_{r_j}-r_{j-1/2}^2{\displaystyle \frac{\partial u_1}{\partial r}}{|}_{r_{j-1/2}}\right)\hfill \\ & & ={\displaystyle \frac{1}{r_j^{2}h}}{\int }_{r_j}^{r_{j+1/2}}{f}_2(r,t)r^{2}dr.\hfill \end{array}$$

(21)

Multiplying (20) by $D_1/d_1$ and (21) by $D_2/d_2$ and summing up the resulting equations and taking into account the interface conditions (6) and (7), we obtain the following:

$$\begin{array}{ccc}& & {\displaystyle \frac{D_1}{d_1{r}_j^{2}h}}\underset{r_{j-1/2}}{\overset{r_j}{\int }}{r}^2{\displaystyle \frac{\partial u_1}{\partial t}}dr+{\displaystyle \frac{D_2}{d_2{r}_j^{2}h}}{\int }_{r_j}^{r_{j+1/2}}{r}^2{\displaystyle \frac{\partial u_2}{\partial t}}dr\hfill \\ & & -{\displaystyle \frac{1}{r_j^{2}h}}\left(D_2{r}_{j+1/2}^2{\displaystyle \frac{\partial u_2}{\partial r}}{|}_{r_{j+1/2}}-D_1{r}_{j-1/2}^2{\displaystyle \frac{\partial u_1}{\partial r}}{|}_{r_{j-1/2}}\right)\hfill \\ & & ={\displaystyle \frac{D_1}{d_1{r}_j^{2}h}}\underset{r_{j-1/2}}{\overset{r_j}{\int }}{f}_1(r,t)r^{2}dr+{\displaystyle \frac{D_2}{d_2{r}_j^{2}h}}{\int }_{r_j}^{r_{j+1/2}}{f}_2(r,t)r^{2}dr.\hfill \end{array}$$

(22)

Next, we use midpoint quadrature for the approximation of the integrals in (18), (19), and (22), replace the fluxes with finite differences discretization

$$r_{j\pm 1/2}^2{\displaystyle \frac{\partial u_i}{\partial r}}{|}_{r_{j\pm 1/2}}=r_{j\pm 1/2}^2{\displaystyle \frac{\pm u_{i,j\pm 1}\mp u_{i,j}}{h}},$$

and involve implicit time-stepping to derive

$$\begin{array}{ccc}& & {\displaystyle \frac{u_{1,j}^{n+1}-u_{1,j}^n}{\tau }}-{\displaystyle \frac{d_1}{r_j^{2}h}}\left(r_{j+1/2}^2{\displaystyle \frac{u_{1,j+1}^{n+1}-u_{1,j}^{n+1}}{h}}-r_{j-1/2}^2{\displaystyle \frac{u_{1,j}^{n+1}-u_{1,j-1}^{n+1}}{h}}\right)\hfill \\ & & =f_{1,j}^{n+1},j=1,2,\dots ,j^{*}-1,\hfill \end{array}$$

(23)

$$\begin{array}{ccc}& & {\displaystyle \frac{D_1}{d_1}}{\displaystyle \frac{u_{1,j}^{n+1}-u_{1,j}^n}{2\tau }}+{\displaystyle \frac{D_2}{d_2}}{\displaystyle \frac{u_{2,j}^{n+1}-u_{2,j}^n}{2\tau }}-{\displaystyle \frac{1}{r_j^{2}h}}\left(D_2{r}_{j+1/2}^2{\displaystyle \frac{u_{2,j+1}^{n+1}-u_{2,j}^{n+1}}{h}}\right.\hfill \\ & & \left.-D_1{r}_{j-1/2}^2{\displaystyle \frac{u_{1,j}^{n+1}-u_{1,j-1}^{n+1}}{h}}\right)={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{D_1}{d_1}}{f}_{1,j}^{n+1}+{\displaystyle \frac{D_2}{d_2}}{f}_{2,j}^{n+1}\right),j=j^{*},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}& & {\displaystyle \frac{u_{2,j}^{n+1}-u_{2,j}^n}{2\tau }}-{\displaystyle \frac{d_2}{r_j^{2}h}}\left(r_{j+1/2}^2{\displaystyle \frac{u_{2,j+1}^{n+1}-u_{2,j}^{n+1}}{h}}-r_{j-1/2}^2{\displaystyle \frac{u_{2,j}^{n+1}-u_{2,j-1}^{n+1}}{h}}\right)\hfill \\ & & =f_{2,j}^{n+1},j=j^{*}+1,\dots ,J-1.\hfill \end{array}$$

(25)

Let $j=0$ ($r=0$). We obtain the approximation for natural BCs (16). Then, we introduce the interval $I_0=(\epsilon ,r_{1/2}=h/2)$ and multiply (14) by $r^2$. Integrating over the interval $I_0$, we get the following:

$$\underset{\epsilon }{\overset{r_{1/2}}{\int }}{r}^2{\displaystyle \frac{\partial u_1}{\partial t}}dr-\left(r_{1/2}^2{d}_1{\displaystyle \frac{\partial u_1}{\partial r}}{|}_{r_{1/2}}-{\epsilon }^2{d}_1{\displaystyle \frac{\partial u_1}{\partial r}}{|}_{\epsilon }\right)=\underset{\epsilon }{\overset{r_{1/2}}{\int }}{r}^2{f}_1(r,t)dr.$$

(26)

Applying the product approximation formula for the integrals in (26), solving the integrals with integrand function $r^2$ exactly, and letting $\epsilon \to 0$ and in view of (16), we obtain the following:

$${\displaystyle \frac{h^3}{24}}{\displaystyle \frac{\partial u_1}{\partial t}}(0,t)-{\displaystyle \frac{d_1{h}^2}{4}}{\displaystyle \frac{\partial u_1}{\partial r}}{|}_{r_{1/2}}={\displaystyle \frac{h^3}{24}}{f}_1(0,t).$$

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

$${\displaystyle \frac{u_{1,0}^{n+1}-u_{1,0}^n}{\tau }}-{\displaystyle \frac{6d_1}{h}}{\displaystyle \frac{u_{1,1}^{n+1}-u_{1,0}^{n+1}}{h}}=f_{1,0}^{n+1}.$$

(27)

The discretization of Dirichlet BCs (4) (if this is the case) is standard

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

(28)

Therefore, we may impose either (27) or (28) for the numerical simulations.

The finite difference scheme is finalized using the discretizations of the right external BC (5) and the initial condition (3), (4)

$$u_{2,J}^{n+1}=u_2(R_2,t_{n+1}),$$

(29)

$$u_{i,j}^0=\left\{\begin{array}{cc}{u}_{10}(r_j),\hfill & j=0,1,\dots ,j^{*},\hfill \\ u_{20}(r_j),\hfill & j=j^{*},j^{*}+1,\dots ,J,\hfill \end{array}\right.u_{10}(R_1)=u_{20}(R_1).$$

(30)

We introduce the vectors

$$\begin{array}{ccc}& & u=[u_1,u_2,\dots ,u_J]=[u_{1,0},u_{1,1},\dots ,u_{1,j^{*}}=u_{2,j^{*}},u_{2,j^{*}+1},\dots ,u_{2,J}],\hfill \\ & & u^0=[u_1^0,u_2^0,\dots ,u_J^0]==[u_{10}(r_0),\dots ,u_{10}(r_{j^{*}})=u_{20}(r_{j^{*}}),u_{20}(r_{j^{*}+1}),\dots ,u_{20}(r_J)]\hfill \end{array}$$

and represent the difference scheme (23)–(25), (27) or (28), (29), and (30) in the equivalent form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & M_0{u}_0^{n+1}-S_0{u}_1^{n+1}=F_0^n,\hfill \\ \hfill -& Q_j{u}_{j-1}^{n+1}+M_j{u}_j^{n+1}-S_j{u}_{j+1}^{n+1}=F_j^n,j=1,2,\dots ,J-1,\hfill \\ & u_J^{n+1}=g_R(t_{n+1}),\hfill \end{array}\end{array}$$

(31)

where $n=0,1,\dots ,N$ and

$$\begin{array}{ccc}& & Q_j={\displaystyle \frac{r_{j-1/2}^2}{r_j^2{h}^2}}\left\{\begin{array}{cc}{d}_1,\hfill & j=1,2,\dots ,j^{*}-1,\hfill \\ D_1,\hfill & j=j^{*},\hfill \\ d_2,\hfill & j=j^{*}+1,\dots ,J-1,\hfill \end{array}\right.,S_j={\displaystyle \frac{r_{j+1/2}^2}{r_j^2{h}^2}}\left\{\begin{array}{cc}{d}_1,\hfill & j=1,2,\dots ,j^{*}-1,\hfill \\ D_2,\hfill & j=j^{*},\hfill \\ d_2,\hfill & j=j^{*}+1,\dots ,J-1,\hfill \end{array}\right.\hfill \\ & & M_j=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\tau }+Q_j+S_j,}\hfill & j=1,2,\dots ,j^{*}-1,\hfill \\ {\displaystyle \frac{1}{2\tau }\left(\frac{D_1}{d_1}+\frac{D_2}{d_2}\right)+Q_j+S_j,}\hfill & j=j^{*},\hfill \\ {\displaystyle \frac{1}{\tau }+Q_j+S_j,}\hfill & j=j^{*}+1,\dots ,J-1,\hfill \end{array}\right.\hfill \\ & & F_j^n=\left\{\begin{array}{cc}{\displaystyle \frac{u_j^n}{\tau }+f_{1,j}^n,}\hfill & j=1,2,\dots ,j^{*}-1,\hfill \\ {\displaystyle \left(\frac{D_1}{d_1}+\frac{D_2}{d_2}\right)\frac{u_j^n}{2\tau }+\frac{1}{2}\left(\frac{D_1}{d_1}{f}_{1,j}^n+\frac{D_2}{d_2}{f}_{2,j}^n\right),}\hfill & j=j^{*},\hfill \\ {\displaystyle \frac{u_j^n}{\tau }+f_{2,j}^n,}\hfill & j=j^{*}+1,\dots ,J-1,\hfill \end{array}\right.\hfill \\ & & M_0=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\tau }+\frac{6d_1}{h^2},}\hfill & \mathrm{BC}(27),\hfill \\ 1,\hfill & \mathrm{BC}(28),\hfill \end{array}\right.S_0=\left\{\begin{array}{cc}{\displaystyle \frac{6d_1}{h^2},}\hfill & \mathrm{BC}(27),\hfill \\ 0,\hfill & \mathrm{BC}(28),\hfill \end{array}\right.\hfill \\ & & F_0^n=\left\{\begin{array}{cc}{\displaystyle \frac{u_j^n}{\tau }+f_{1,0}^{n+1},}\hfill & \mathrm{BC}(27),\hfill \\ g_0(t_{n+1}),\hfill & \mathrm{BC}(28).\hfill \end{array}\right.\hfill \end{array}$$

The coefficient matrix of the system (31) is a tri-diagonal M-matrix. Namely it is strictly diagonally dominant with positive main diagonal entries and non-positive off-diagonal elements. With the same line of considerations as in [29,46], we may conclude that the order of convergence of the discrete scheme (31) is $O(\tau +h^2)$.

5. Numerical Solution of the Inverse Problems

In this section, we develop an approach for numerically solving inverse problems IP1a–IP3c. On a discrete level, we apply the decomposition of the solution with respect to the unknown boundary function(s) to construct discretization for recovering external boundary conditions.

First, we consider IP1a and IP1b. In this case, we impose left BCs (5), consider the following representation of the solution $u_j^{n+1}$

$$\begin{array}{ccc}\hfill & & u_j^{n+1}=y_j^{n+1}+g_0^{n+1}{v}_j,\hfill \end{array}$$

(32)

and substitute $u_j^{n+1}$ into (31). The resulting problems for the unknown functions $y_j^{n+1}$ and $v_j^{n+1}$ are

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & y_0^{n+1}=0,\hfill \\ \hfill -& Q_j{y}_{j-1}^{n+1}+M_j{y}_j^{n+1}-S_i{y}_{j+1}^{n+1}=F_j^n,j=1,2,\dots ,J-1,\hfill \\ & y_J^{n+1}=g_R(t_{n+1}),\hfill \\ & u_j^0=u^0(j),\hfill \end{array}\end{array}$$

(33)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & v_0^{n+1}=1,\hfill \\ \hfill -& Q_j{v}_{j-1}+M_j{v}_j-S_j{v}_{j+1}=0,j=1,2,\dots ,J-1,\hfill \\ & v_J=0.\hfill \end{array}\end{array}$$

(34)

We solve (33) and (34) in order to obtain y and v in the discrete domain ${\overline{\omega }}_h\times {\overline{\omega }}_{\tau }$. Further, using (32) and the measurement $u(r^{*},t)=\psi (t)$, $r^{*}\in {\mathsf{\Omega }}_1\cup {\mathsf{\Omega }}_2$, we determine

$$g_0^{n+1}={\displaystyle \frac{{\psi }^{n+1}-y_{r^{*}}^{n+1}}{v_{r^{*}}}}.$$

(35)

Here, $y_{r^{*}}^{n+1}=y^{n+1}(r^{*})$ and $v_{r^{*}}=v^{n+1}(r^{*})$ if $r^{*}$ is a mesh point; otherwise, we apply interpolation.

Similarly, we develop the numerical algorithm for solving IP2a and IP2b. In this case, both natural (16) and Dirichlet (5) BCs can be considered at $r=0$. We then introduce the decomposition

$$\begin{array}{ccc}\hfill & & u_j^{n+1}=y_j^{n+1}+g_R^{n+1}{v}_j,\hfill \end{array}$$

(36)

and substitute it into (31) to get the following two problems:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & M_0{y}_0-S_0{y}_1=F_1^n,\hfill \\ \hfill -& Q_j{y}_{j-1}^{n+1}+M_j{y}_j^{n+1}-S_i{y}_{j+1}^{n+1}=F_j^n,j=1,2,\dots ,J-1,\hfill \\ & y_J^{n+1}=0,\hfill \\ & u_j^0=u^0(j),\hfill \end{array}\end{array}$$

(37)

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & M_0{v}_0-S_0{v}_1=0,\hfill \\ \hfill -& Q_j{v}_{j-1}+M_j{v}_j-S_j{v}_{j+1}=0,j=1,2,\dots ,J-1,\hfill \\ & v_J=1.\hfill \end{array}\end{array}$$

(38)

Solving (37) and (38) to find y and v in the entire computational region, and from (36) along with the given observations $u(r^{*},t)=\psi (t)$, $r^{*}\in {\mathsf{\Omega }}_1\cup {\mathsf{\Omega }}_2$, we find the right BC

$$g_R^{n+1}={\displaystyle \frac{{\psi }^{n+1}-y_{r^{*}}^{n+1}}{v_{r^{*}}^{n+1}}}.$$

(39)

Next, we solve the inverse problems IP3a–IP3c numerically. To recover $g_0$, we need to impose the Dirichlet BC (5) at $r=0$. Alternatively, we may consider the natural BC (16) and use the decomposition (36) to restore $g_R$, applying (37)–(39). In this approach, $g_0$ is determined as a part of the recovered solution u, namely $g_0^n=u_0^n$. We consider the case of Dirichlet BC (5) at $r=0$.

We define data points that encompass the related measurements for IP3a–IP3c

$$u(r_i^{*})={\psi }_i(t),r_i^{*}\in {\mathsf{\Omega }}_1\cup {\mathsf{\Omega }}_2,i=1,2.$$

(40)

Further, we introduce the decomposition

$$\begin{array}{ccc}\hfill & & u_j^{n+1}=y_j^{n+1}+g_R^{n+1}{v}_j+g_0^{n+1}{w}_j\hfill \end{array}$$

(41)

and replace it with (31) to obtain the following three problems, corresponding to the unknown functions $y_j^{n+1}$, $v_j^{n+1}$, and $w_j^{n+1}$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & y_0^{n+1}=0,\hfill \\ \hfill -& Q_j{y}_{j-1}^{n+1}+M_j{y}_j^{n+1}-S_j{y}_{j+1}^{n+1}=F_j^n,j=1,2,\dots ,J-1,\hfill \\ & y_J^{n+1}=0,\hfill \\ & u_j^0=u^0(j),\hfill \end{array}\end{array}$$

(42)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & v_0=0,\hfill \\ \hfill -& Q_j{v}_{j-1}+M_j{v}_j-S_j{v}_{j+1}=0,j=1,2,\dots ,J-1,\hfill \\ & v_J=1.\hfill \end{array}\end{array}$$

(43)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & w_0=1,\hfill \\ \hfill -& Q_j{w}_{j-1}+M_j{w}_j-S_j{w}_{j+1}=0,j=1,2,\dots ,J-1,\hfill \\ & w_J=0.\hfill \end{array}\end{array}$$

(44)

Solving these problems, from (41) and in view of (40), we get the following:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\psi }_1^{n+1}=y_{r_1^{*}}^{n+1}+g_R^{n+1}{v}_{r_1^{*}}+g_0^{n+1}{w}_{r_1^{*}},\hfill \\ & {\psi }_2^{n+1}=y_{r_2^{*}}^{n+1}+g_R^{n+1}{v}_{r_2^{*}}+g_0^{n+1}{w}_{r_2^{*}}.\hfill \end{array}\end{array}$$

(45)

Consequently, from (45), we derive

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & g_0^{n+1}={\displaystyle \frac{v_{r_1^{*}}({\psi }_2^{n+1}-y_{r_2^{*}}^{n+1})-v_{r_2^{*}}({\psi }_1^{n+1}-y_{r_1^{*}}^{n+1})}{v_{r_1^{*}}{w}_{r_2^{*}}-v_{r_2^{*}}{w}_{r_1^{*}}}},\hfill \\ & g_R^{n+1}={\displaystyle \frac{({\psi }_1^{n+1}-y_{r_1^{*}}^{n+1})w_{r_2^{*}}-({\psi }_2^{n+1}-y_{r_2^{*}}^{n+1})w_{r_1^{*}}}{v_{r_1^{*}}{w}_{r_2^{*}}-v_{r_2^{*}}{w}_{r_1^{*}}}}.\hfill \end{array}\end{array}$$

(46)

6. Computational Test Examples

In this section, we illustrate the efficiency of the developed methods for solving both direct and inverse problems. The test example is (1)–(7), where

$$\begin{array}{cc}& {\displaystyle R_2=3,R_1=1,T=1,d_1=1,d_2=10,D_1=0.05,D_2=3,}\\ & {\displaystyle g_0(t)=e^{-t/2},g_R(t)=\sqrt{2}\frac{40-9\pi }{80}{e}^{-t/2},}\\ & {\displaystyle u_{10}(r)=cos\frac{\pi x}{4},u_{20}(r)=\frac{\sqrt{2}}{2}\left(\frac{\pi (x^3-x^4)}{240}+1\right).}\end{array}$$

The functions $f_1(r,t)$ and $f_2(r,t)$ are chosen such that $u=(u_1,u_2)$, where ${\displaystyle u_1(r,t)=e^{-t/2}}$ $cos{\displaystyle \frac{\pi x}{4}}$, ${\displaystyle u_2(r,t)=\frac{\sqrt{2}}{2}\left(\frac{\pi (x^3-x^4)}{240}+1\right)e^{-t/2}}$ is the exact solution of the problem (1)–(7).

We involve the following notations: $u(r,t)=\left(u_1(r,t),u_2(r,t)\right)$ is the exact solution of the problem (1)–(7), $g_0(t_n)$ and $g_R(t)$ are exact boundary functions, $u_j^n=(u_{1,j}^n,u_{1,j}^n)$ is the numerical solution of the direct problem at point $(x_j,t_n)$ computed by (31), $g_0^n$ and $g_R^n$ are numerically recovered boundary functions at time $t_n$, and $U_{i,j}^n=(U_{1,j}^n,U_{2,j}^n)$ is the numerical solution at grid node $(r_j,t_n)$ obtained by solving inverse problems IP1a–IP3c.

The error of the numerical solution of the direct problem at final lime is defined by

$${\mathcal{E}}_{i,j}=u_{i,j}^N-u_i(r_j,T),i=1,2,j=0,1,\dots ,J,n=0,1,\dots ,N,$$

while for the inverse problems, the errors of the recovered solution and boundary functions are computed as follows:

$${\mathcal{E}}_{i,j}=U_{i,j}^N-u_{i,j}^N,i=1,2,j=0,1,\dots ,J,{\epsilon }_{s,n}=g_s^n-g_s(t_n),n=0,1,\dots ,N,s=\{0,R\}.$$

Consequently, the errors in maximal and $L_2$ norms and the corresponding orders of convergence of the numerical solution of the direct and inverse problems IP1a–IP3c are given by

$$\begin{array}{cc}& {\displaystyle E_{1,\infty }=E_{1,\infty }(J)=\underset{0\le j\le j^{*}}{max}|{\mathcal{E}}_{1,j}|,E_{2,\infty }=E_{2,\infty }(J)=\underset{j^{*}\le jJ}{max}\left|{\mathcal{E}}_{2,j}\right|,}\\ & {\displaystyle E_{1,L_2}=E_{1,L_2}(J)={\left[\sum _{j=0}^{j^{*}}h{({\mathcal{E}}_{1,j})}^2\right]}^{1/2},E_{2,L_2}=E_{2,L_2}(J)={\left[\sum _{j=j^{*}}^{J}h{({\mathcal{E}}_{2,j})}^2\right]}^{1/2},}\\ & {\displaystyle C{R}_{i,\infty }={log}_2\frac{E_{i,\infty }(J)}{E_{i,\infty }(2J)},C{R}_{i,L_2}={log}_2\frac{E_{i,L_2}(J)}{E_{i,L_2}(2J)},i=1,2,}\\ & {\displaystyle e_{s,\infty }=e_{s,\infty }(N)=\underset{0\le n\le N}{max}\left|{\epsilon }_{s,n}\right|,e_{s,L_2}=e_{s,L_2}(N)={\left[\sum _{n=0}^{N}h{\epsilon }_{s,n}^2\right]}^{1/2},}\\ & {\displaystyle c{r}_{s,\infty }={log}_2\frac{e_{s,\infty }(N)}{e_{s,\infty }(2N)},c{r}_{s,L_2}={log}_2\frac{e_{s,L_2}(N)}{e_{s,L_2}(2N)}.}\end{array}$$

Example 1

(Direct problem: convergence). First, we verify the convergence rate of the numerical solution for the direct problem computed using (31). Computations are carried out with a fixed ratio $\tau =h^2$. We test both the Dirichlet BC (5) and the natural BC (16) at $r=0$. The results are presented in

Table 1. Errors and convergence rates, direct problem, and Dirichlet BC: Example 1.

J	${\mathit{E}}_{1,\mathbf{\infty }}$	${\mathit{CR}}_{1,\mathbf{\infty }}$	${\mathit{E}}_{2,\mathbf{\infty }}$	${\mathit{CR}}_{2,\mathbf{\infty }}$	${\mathit{E}}_{1,{\mathit{L}}_2}$	${\mathit{CR}}_{1,{\mathit{L}}_2}$	${\mathit{E}}_{2,{\mathit{L}}_2}$	${\mathit{CR}}_{2,{\mathit{L}}_2}$
30	1.9891 $\times {10}^{-3}$		8.7711 $\times {10}^{-4}$		1.3559 $\times {10}^{-3}$		8.0890 $\times {10}^{-4}$
60	6.2112 $\times {10}^{-4}$	1.6792	2.1954 $\times {10}^{-4}$	1.9983	3.7743 $\times {10}^{-4}$	1.8449	1.9933 $\times {10}^{-4}$	2.0208
120	1.8709 $\times {10}^{-4}$	1.7311	5.4910 $\times {10}^{-5}$	1.9993	1.0188 $\times {10}^{-4}$	1.8894	4.9466 $\times {10}^{-5}$	2.0107
240	5.4842 $\times {10}^{-5}$	1.7704	1.3731 $\times {10}^{-5}$	1.9997	2.6861 $\times {10}^{-5}$	1.9232	1.2320 $\times {10}^{-5}$	2.0054
480	1.5869 $\times {10}^{-5}$	1.7891	3.4331 $\times {10}^{-6}$	1.9998	6.9625 $\times {10}^{-6}$	1.9479	3.0743 $\times {10}^{-6}$	2.0027
960	4.5344 $\times {10}^{-6}$	1.8072	8.5833 $\times {10}^{-7}$	1.9999	1.7831 $\times {10}^{-6}$	1.9652	7.6787 $\times {10}^{-7}$	2.0014
1920	1.2758 $\times {10}^{-6}$	1.8296	2.1459 $\times {10}^{-7}$	1.9999	4.5287 $\times {10}^{-7}$	1.9772	1.9188 $\times {10}^{-7}$	2.0007

and

Table 2. Errors and convergence rates, direct problem, and natural BC: Example 1.

J	${\mathit{E}}_{1,\mathbf{\infty }}$	${\mathit{CR}}_{1,\mathbf{\infty }}$	${\mathit{E}}_{2,\mathbf{\infty }}$	${\mathit{CR}}_{2,\mathbf{\infty }}$	${\mathit{E}}_{1,{\mathit{L}}_2}$	${\mathit{CR}}_{1,{\mathit{L}}_2}$	${\mathit{E}}_{2,{\mathit{L}}_2}$	${\mathit{CR}}_{2,{\mathit{L}}_2}$
30	2.4116 $\times {10}^{-3}$		8.7772 $\times {10}^{-4}$		1.6661 $\times {10}^{-3}$		8.0927 $\times {10}^{-4}$
60	7.5982 $\times {10}^{-4}$	1.6663	2.1963 $\times {10}^{-4}$	1.9987	4.4026 $\times {10}^{-4}$	1.9200	1.9939 $\times {10}^{-4}$	2.0210
120	2.2986 $\times {10}^{-4}$	1.7249	5.4925 $\times {10}^{-5}$	1.9995	1.1379 $\times {10}^{-4}$	1.9520	4.9475 $\times {10}^{-5}$	2.0108
240	6.7518 $\times {10}^{-5}$	1.7674	1.3733 $\times {10}^{-5}$	1.9998	2.9009 $\times {10}^{-5}$	1.9717	1.2322 $\times {10}^{-5}$	2.0055
480	1.9403 $\times {10}^{-5}$	1.7990	3.4334 $\times {10}^{-6}$	1.9999	7.3351 $\times {10}^{-6}$	1.9836	3.0745 $\times {10}^{-6}$	2.0028
960	5.4826 $\times {10}^{-6}$	1.8233	8.5837 $\times {10}^{-7}$	2.0000	1.8457 $\times {10}^{-6}$	1.9906	7.6789 $\times {10}^{-7}$	2.0014
1920	1.5288 $\times {10}^{-6}$	1.8425	2.1460 $\times {10}^{-7}$	2.0000	4.6313 $\times {10}^{-7}$	1.9947	1.9188 $\times {10}^{-7}$	2.0007

. Clearly, in both cases (Dirichlet and natural BCs at $r=0$), the accuracy is $O(\tau +h^2)$, and the precision is almost the same.

Example 2

(Inverse problems IP1: exact measurements). We test the efficiency of the numerical approach (32)–(35) for solving IP1a and IP1b. The measurements at points $r_1^{*}$ and $r_2^{*}$ are obtained from the numerical solution of the direct problem, namely ${\psi }_i^n=u^n(r_i^{*})$, where $r_i^{*}$ and $i=1,2$ are grid nodes. All runs are conducted for $J=240$ and $\tau =h$. In this case, the condition number of the coefficient matrix in (31) is $8.5\times {10}^6$. Note that the coefficient matrix of the discrete systems for all methods for recovering $g_R$ or/and $g_0$ is almost the same as the one in (31). The only difference is in the first and last row.

In Figure 1 and Figure 2, we plot the exact $g_0(t_n)$ and recovered boundary functions $g_0^n$, $n=0,$$1,\dots ,N$, the numerical solution $u_{i,j}^N$ and the recovered solution $U_{i,j}^N$, and the corresponding errors ${\epsilon }_{0,n}$, $n=0,1,\dots ,N$, $\mathcal{E}=({\mathcal{E}}_{1,j_1},{\mathcal{E}}_{2,j_2})$, $j_1=0,1,\dots ,j^{*}$, $j_2=j^{*},j^{*}+1,\dots ,J$ in the case $r_1^{*}=0.5$ (i.e., IP1a). Similar types of graphics, but for $r_2^{*}=2.5$ (i.e., IP1b), are shown in Figure 3 and Figure 4.

We observe that if the measurement is close to $r=0$, the solution and boundary functions are recovered almost exactly; therefore, the numerical algorithm (32)–(35) for solving IP1a has the same order of convergence as the discrete scheme (31) for solving the direct problem. In this case, when the measurements are away from the boundary $r=0$ (i.e., IP1b), the boundary function $g_0$ and the numerical solution U of the inverse problem are restored with lower (but sufficiently good) precision.

Example 3

(Inverse problems IP2: exact measurements). We verify the performance of the numerical method (36)–(39), using the natural BC (16) for solving IP2a and IP2b. We study the influence of the location of the measurement on the accuracy of the recovered boundary function $g_R$ and the solution U. For the simulations, we generate the observations, as shown in Example 2. In

Table 3. Errors of the recovered solution U and function $g_L$, IP2a and IP2b. Example 3.

IP	${\mathit{r}}_{\mathit{i}}^{*}$	${\mathit{E}}_{1,\mathbf{\infty }}$	${\mathit{E}}_{2,\mathbf{\infty }}$	${\mathit{E}}_{1,{\mathit{L}}_2}$	${\mathit{E}}_{2,{\mathit{L}}_2}$	${\mathit{e}}_{\mathit{R},\mathbf{\infty }}$	${\mathit{e}}_{\mathit{R},2}$
IP2a	0.25	1.8374 $\times {10}^{-13}$	8.0447 $\times {10}^{-11}$	4.0118 $\times {10}^{-14}$	3.2759 $\times {10}^{-11}$	2.3996 $\times {10}^{-10}$	9.3420 $\times {10}^{-11}$
IP2a	0.5	3.4361 $\times {10}^{-14}$	7.0244 $\times {10}^{-12}$	7.6707 $\times {10}^{-15}$	3.1521 $\times {10}^{-12}$	1.8443 $\times {10}^{-11}$	9.7559 $\times {10}^{-12}$
IP2a	0.95	3.3307 $\times {10}^{-16}$	7.9675 $\times {10}^{-13}$	1.3441 $\times {10}^{-16}$	2.5261 $\times {10}^{-13}$	7.4146 $\times {10}^{-12}$	3.0400 $\times {10}^{-12}$
IP2b	1.5	3.5527 $\times {10}^{-15}$	5.3402 $\times {10}^{-14}$	2.7998 $\times {10}^{-15}$	1.9208 $\times {10}^{-14}$	1.3461 $\times {10}^{-13}$	5.2055 $\times {10}^{-14}$
IP2b	2	6.1062 $\times {10}^{-15}$	2.1649 $\times {10}^{-15}$	5.4508 $\times {10}^{-15}$	1.5495 $\times {10}^{-15}$	9.0206 $\times {10}^{-15}$	4.1202 $\times {10}^{-15}$
IP2b	2.5	2.3315 $\times {10}^{-15}$	2.1649 $\times {10}^{-15}$	1.6722 $\times {10}^{-15}$	1.1177 $\times {10}^{-15}$	2.7200 $\times {10}^{-15}$	7.0581 $\times {10}^{-16}$

, we present errors of the recovered solution U and function $g_R$, $J=240$, $\tau =h$.

Our results indicate that the boundary function is restored almost exactly; similarly, the solution $U_j^n$, determined by IP2, recovers almost exactly the numerical solution of the direct problem $u_j^n$. The position of the measurements slightly influences the accuracy of $g_R$ and U. Specifically, when the measurements are near the right boundary, the precision of the recovered data is higher compared to when the measurements are distant from $r=R_2$. In other words, the closer the measurements are to the boundary, the greater the precision of the corresponding recovered boundary data. The reason is that the solution v to the problem (38) is strictly positive $0<v\le 1$ (see [35,43]), increasing function. Consequently, far from the right boundary, the value of $v_{r^{*}}$ can become very small, potentially introducing additional error into formula (39) due to division.

Note that when using numerical methods to solve IP2a and IP2b with natural BC, we achieve a higher precision in also restoring $g_0^n=U_0^n$ and $n=0,1,\dots ,N$ compared to the precision obtained with IP1a and IP1b. Therefore, employing natural BC for solving inverse problems results in more successful recovery than algorithms that utilize Dirichlet boundary conditions at $r=0$.

Example 4

(Inverse problems IP3: exact measurements). In this example, we illustrate the efficiency of the numerical approach (41)–(46) for solving IP3a–IP3c. The mesh and the measurements are generated, as in Example 2. The computational results are given in

Table 4. Errors of the recovered solution U and functions $g_0$, $g_R$, IP3a–IP3c. Example 4.

IP	${\mathit{r}}_1^{*}$	${\mathit{r}}_2^{*}$	${\mathit{E}}_{1,\mathbf{\infty }}$	${\mathit{E}}_{2,\mathbf{\infty }}$	${\mathit{E}}_{1,{\mathit{L}}_2}$	${\mathit{E}}_{2,{\mathit{L}}_2}$	${\mathit{e}}_{0,\mathbf{\infty }}$	${\mathit{e}}_{0,2}$	${\mathit{e}}_{\mathit{R},\mathbf{\infty }}$	${\mathit{e}}_{\mathit{R},2}$
IP3a	0.25	1.5	1.7 $\times {10}^{-12}$	3.6 $\times {10}^{-14}$	1.9 $\times {10}^{-13}$	1.4 $\times {10}^{-14}$	3.5 $\times {10}^{-12}$	1.1 $\times {10}^{-12}$	1.5 $\times {10}^{-13}$	5.6 $\times {10}^{-14}$
IP3a	0.5	2	9.3 $\times {10}^{-12}$	2.8 $\times {10}^{-15}$	1.1 $\times {10}^{-12}$	1.3 $\times {10}^{-15}$	2.4 $\times {10}^{-11}$	1.1 $\times {10}^{-11}$	1.3 $\times {10}^{-14}$	4.8 $\times {10}^{-15}$
IP3a	0.95	1.05	4.9 $\times {10}^{-8}$	4.9 $\times {10}^{-13}$	5.7 $\times {10}^{-9}$	1.6 $\times {10}^{-13}$	2.1 $\times {10}^{-7}$	8.4 $\times {10}^{-8}$	8.4 $\times {10}^{-12}$	3.1 $\times {10}^{-12}$
IP3b	0.5	0.75	3.1 $\times {10}^{-12}$	4.6 $\times {10}^{-12}$	3.5 $\times {10}^{-13}$	1.7 $\times {10}^{-12}$	1.7 $\times {10}^{-11}$	7.5 $\times {10}^{-12}$	1.1 $\times {10}^{-11}$	5.1 $\times {10}^{-12}$
IP3b	0.25	0.95	3.2 $\times {10}^{-13}$	3.1 $\times {10}^{-13}$	3.6 $\times {10}^{-14}$	2.3 $\times {10}^{-13}$	2.4 $\times {10}^{-12}$	7.7 $\times {10}^{-13}$	7.6 $\times {10}^{-12}$	3.0 $\times {10}^{-12}$
IP3b	0.25	0.75	5.9 $\times {10}^{-13}$	2.4 $\times {10}^{-12}$	6.8 $\times {10}^{-14}$	8.2 $\times {10}^{-13}$	2.1 $\times {10}^{-12}$	7.9 $\times {10}^{-13}$	1.3 $\times {10}^{-11}$	4.4 $\times {10}^{-12}$
IP3c	1	2	4.4 $\times {10}^{-6}$	3.4 $\times {10}^{-15}$	5.1 $\times {10}^{-7}$	2.3 $\times {10}^{-15}$	1.9 $\times {10}^{-5}$	7.9 $\times {10}^{-6}$	1.8 $\times {10}^{-14}$	6.5 $\times {10}^{-15}$
IP3c	1.5	2.5	6.8 $\times {10}^{-7}$	8.4 $\times {10}^{-15}$	7.7 $\times {10}^{-8}$	2.1 $\times {10}^{-15}$	7.8 $\times {10}^{-6}$	2.8 $\times {10}^{-6}$	1.4 $\times {10}^{-15}$	5.7 $\times {10}^{-16}$
IP3c	1.5	2.95	4.0 $\times {10}^{-7}$	6.2 $\times {10}^{-15}$	4.6 $\times {10}^{-8}$	1.7 $\times {10}^{-15}$	8.9 $\times {10}^{-6}$	3.3 $\times {10}^{-6}$	1.1 $\times {10}^{-16}$	3.8 $\times {10}^{-17}$

.

We observe that the solution and boundary functions determined by IP3 exhibit high accuracy. The results indicate that lower (but at the same time good enough) precision is achieved in recovering $g_0$ using IP3a when the measurements are taken far from $r=0$ and IP3c, where both measurements are located within ${\mathsf{\Omega }}_2$. As noted previously, this can be explained by the values of the functions v and w, which solve problems (43)–(44), respectively.

When the measurements are far from $z=0$, the precision of recovering $g_0$ is lower compared to that of $g_R$ when the measurements are far from the right boundary. Thus, the accuracy of the restored $g_0$ is influenced by the measurement locations, whereas the precision of $g_R$ remains almost unaffected by their position. This also influences the accuracy of the recovered solution U. Specifically, the precision of the determined $g_0^n$ primarily impacts the accuracy of the solution $U_1^n$, while the precision of $g_R^n$ affects the accuracy of the solution $U_2^n$.

Example 5

(Inverse problems: noisy data). We will demonstrate the accuracy of the numerically recovering boundary function and the solution for perturbed observations. As we mentioned before, the use of natural BC at degenerate boundaries leads to better precision of the restored solution and boundary functions. Moreover, IP2a and IP2b also recover $g_0$ as a part of the solution. Therefore, we examine numerical methods for solving IP2a and IP2b with natural BC at $r=0$ and IP3a–IP3c. We repeat the simulations from Examples 3 and 4 but add noise to the observed data:

$${\psi }_{i,j}^n=\left(1+2{\rho }_i({\sigma }_i^n-0.5)\right)u_{i,j_i}^n,r_i^{*}=r_{j_i},i=1,2,j=1,2,\dots ,J-1,$$

where ${\rho }_i$ represents the noise levels, and ${\sigma }_i^n$ denotes random values that are uniformly distributed over the interval $[0,1]$. In order to smooth the measured data, we apply polynomial curve fitting of degree 5.

The simulation results for IP2 and IP3 with a variety of noise levels are shown in

Table 5. Errors of the recovered solution and boundary functions $g_0$ and $g_R$ for different levels of noise, IP3. Example 5.

${\mathit{\rho }}_1$	${\mathit{\rho }}_2$	${\mathit{E}}_{1,\mathbf{\infty }}$	${\mathit{E}}_{2,\mathbf{\infty }}$	${\mathit{E}}_{1,{\mathit{L}}_2}$	${\mathit{E}}_{2,{\mathit{L}}_2}$	${\mathit{e}}_{0,\mathbf{\infty }}$	${\mathit{e}}_{0,2}$	${\mathit{e}}_{\mathit{R},\mathbf{\infty }}$	${\mathit{e}}_{\mathit{R},2}$
		IP3a, $r_1^{*}=0.5$, $r_2^{*}=2$
0.001	0.003	2.80 $\times {10}^{-1}$	8.25 $\times {10}^{-4}$	3.21 $\times {10}^{-2}$	6.46 $\times {10}^{-4}$	1.18	1.89 $\times {10}^{-1}$	3.47 $\times {10}^{-3}$	5.27 $\times {10}^{-4}$
0.002	0.03	1.26	8.21 $\times {10}^{-3}$	1.45 $\times {10}^{-1}$	6.44 $\times {10}^{-3}$	2.62	5.19 $\times {10}^{-1}$	3.49 $\times {10}^{-2}$	5.31 $\times {10}^{-3}$
		IP3b, $r_1^{*}=0.25$, $r_2^{*}=0.75$
0.001	0.003	7.95 $\times {10}^{-2}$	3.10 $\times {10}^{-3}$	9.10 $\times {10}^{-3}$	2.49 $\times {10}^{-3}$	7.95 $\times {10}^{-2}$	2.20 $\times {10}^{-2}$	8.11 $\times {10}^{-1}$	1.41 $\times {10}^{-1}$
0.002	0.03	5.28 $\times {10}^{-1}$	2.96 $\times {10}^{-2}$	6.05 $\times {10}^{-2}$	2.46 $\times {10}^{-2}$	5.27 $\times {10}^{-1}$	1.83 $\times {10}^{-1}$	8.11	1.41
0.03	0.002	1.01	3.21 $\times {10}^{-3}$	1.16 $\times {10}^{-1}$	2.41 $\times {10}^{-3}$	1.51	3.43 $\times {10}^{-1}$	5.39 $\times {10}^{-1}$	9.38 $\times {10}^{-2}$
		IP3c, $r_1^{*}=1.5$, $r_2^{*}=2.5$
0.001	0.003	fails	5.14 $\times {10}^{-4}$	fails	3.74 $\times {10}^{-4}$	fails	fails	8.17 $\times {10}^{-4}$	1.57 $\times {10}^{-4}$
0.0005	0.1	fails	1.56 $\times {10}^{-2}$	fails	1.12 $\times {10}^{-2}$	fails	fails	2.76 $\times {10}^{-2}$	5.16 $\times {10}^{-3}$

and

Table 6. Errors of the recovered solution and boundary functions $g_0$ and $g_R$ for different levels of noise, IP2. Example 5.

${\mathit{\rho }}_{\mathit{i}}$	${\mathit{E}}_{1,\mathbf{\infty }}$	${\mathit{E}}_{2,\mathbf{\infty }}$	${\mathit{E}}_{1,{\mathit{L}}_2}$	${\mathit{E}}_{2,{\mathit{L}}_2}$	${\mathit{e}}_{0,\mathbf{\infty }}$	${\mathit{e}}_{0,2}$	${\mathit{e}}_{\mathit{R},\mathbf{\infty }}$	${\mathit{e}}_{\mathit{R},2}$
	IP2a, $r_2^{*}=1.05$
0.001	5.53 $\times {10}^{-5}$	1.10 $\times {10}^{-4}$	5.25 $\times {10}^{-5}$	6.24 $\times {10}^{-5}$	2.41 $\times {10}^{-4}$	1.64 $\times {10}^{-4}$	2.20 $\times {10}^{-2}$	3.49 $\times {10}^{-3}$
0.01	2.80 $\times {10}^{-4}$	1.27 $\times {10}^{-3}$	9.89 $\times {10}^{-5}$	9.54 $\times {10}^{-4}$	7.07 $\times {10}^{-4}$	5.17 $\times {10}^{-4}$	1.38 $\times {10}^{-1}$	2.20 $\times {10}^{-2}$
0.1	1.59 $\times {10}^{-1}$	3.39 $\times {10}^{-2}$	7.59 $\times {10}^{-3}$	3.26 $\times {10}^{-2}$	3.70 $\times {10}^{-3}$	2.03 $\times {10}^{-3}$	4.47 $\times {10}^{-1}$	7.19 $\times {10}^{-2}$
	IP2a, $r_2^{*}=1.5$
0.001	9.24 $\times {10}^{-5}$	3.64 $\times {10}^{-4}$	3.46 $\times {10}^{-5}$	2.81 $\times {10}^{-4}$	2.02 $\times {10}^{-4}$	1.62 $\times {10}^{-4}$	4.17 $\times {10}^{-3}$	6.48 $\times {10}^{-4}$
0.01	6.27 $\times {10}^{-4}$	2.19 $\times {10}^{-3}$	2.34 $\times {10}^{-4}$	1.75 $\times {10}^{-3}$	3.33 $\times {10}^{-4}$	1.78 $\times {10}^{-4}$	3.92 $\times {10}^{-2}$	5.89 $\times {10}^{-3}$
0.1	1.24 $\times {10}^{-3}$	5.27 $\times {10}^{-3}$	7.27 $\times {10}^{-4}$	2.84 $\times {10}^{-3}$	5.19 $\times {10}^{-3}$	3.16 $\times {10}^{-3}$	7.52 $\times {10}^{-1}$	1.16 $\times {10}^{-1}$
	IP2b, $r_1^{*}=0.5$
0.001	6.16 $\times {10}^{-5}$	2.63 $\times {10}^{-4}$	5.52 $\times {10}^{-5}$	1.51 $\times {10}^{-4}$	2.53 $\times {10}^{-4}$	1.82 $\times {10}^{-4}$	1.39	2.64 $\times {10}^{-1}$
0.01	4.02 $\times {10}^{-4}$	2.32 $\times {10}^{-3}$	3.45 $\times {10}^{-4}$	1.34 $\times {10}^{-3}$	6.65 $\times {10}^{-4}$	3.84 $\times {10}^{-4}$	8.98	1.73
	IP2b, $r_1^{*}=0.95$
0.001	1.94 $\times {10}^{-5}$	8.77 $\times {10}^{-5}$	1.34 $\times {10}^{-5}$	4.72 $\times {10}^{-5}$	2.32 $\times {10}^{-4}$	1.93 $\times {10}^{-4}$	2.08 $\times {10}^{-2}$	3.49 $\times {10}^{-3}$
0.01	5.73 $\times {10}^{-4}$	2.70 $\times {10}^{-3}$	5.15 $\times {10}^{-4}$	1.61 $\times {10}^{-3}$	3.85 $\times {10}^{-4}$	1.92 $\times {10}^{-4}$	2.47 $\times {10}^{-2}$	4.11 $\times {10}^{-3}$
0.1	2.75 $\times {10}^{-3}$	2.65 $\times {10}^{-3}$	2.22 $\times {10}^{-3}$	1.67 $\times {10}^{-3}$	3.19 $\times {10}^{-3}$	1.82 $\times {10}^{-3}$	2.29	3.85 $\times {10}^{-1}$

. In Figure 5, Figure 6 and Figure 7, we plot exact and recovered boundary functions, numerical solutions u, U of the direct and inverse problems, recovered by IP3b, $r_1^{*}=0.25$, $r_2^{*}=0.75$, ${\rho }_1=0.001$, ${\rho }_2=0.003$, and the corresponding errors. In Figure 8, Figure 9 and Figure 10, we illustrate the same type of graphics, but obtained by IP2a for $r_2^{*}=1.5$, ${\rho }_2=0.01$.

Regarding IP3, our numerical results show that the method is more efficient when the measurements are taken in the region ${\mathsf{\Omega }}_1$, i.e., IP3b. In this case, both functions $g_0$ and $g_R$ and the solution U are recovered with optimal accuracy. The cases of big errors in the maximum norm are due to the higher error of the corresponding recovered functions close to the initial time. Then, moving away from $t=0$, the function is recovered with satisfactory accuracy, and this affects the error in the $L_2$ norm. This is illustrated in Figure 5 and Figure 6. Consequently, the determined solution U at the final time is accurate enough, see Figure 7.

If the measurements are located in both domains ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$ (IP3a), recovery is successful for small levels of noise. Those determined by the IP3c function $g_R$ and solution $U_2$ are of high precision, but the restoration of function $g_0$ and solution $U_1$ fails even at lower noise levels.

Much more successful is the determination of the boundary functions and the solution performed with IP2 in comparison to IP3, see Table 6. A good precision of the recovered functions $g_0$, $g_R$ and solution U is obtained when the measurements are within the domain ${\mathsf{\Omega }}_2$, i.e., IP2a, including the case of high noise levels. If the measurements are located in the domain ${\mathsf{\Omega }}_1$, they must be at points nearer to the interface, such that to be closer to the domain ${\mathsf{\Omega }}_2$ and, respectively, to the right boundary, in order to attain satisfactory precision of the restored functions. Again, the largest error of the restored functions $g_R$ appears close to the initial time. See Figure 9.

Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 illustrate the impact of noisy data on solution accuracy. Figure 1, Figure 2, Figure 3 and Figure 4 show a close fit between the exact and recovered data, while Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 display the effect of increased noise levels in the measurements on both the recovered boundary function and the solution, resulting in larger errors.

7. Conclusions

In this study, we construct an effective numerical technique for addressing inverse problems aimed at retrieving the solution along with the right and/or left Dirichlet boundary conditions in an interface parabolic equation in spherical geometry. We establish the existence and the uniqueness of the solution of the direct differential problem and develop a numerical method that is second-order in space and first-order in time for solving it. Given the degeneracy at the left boundary $r=0$, we examine both Dirichlet and natural boundary conditions. We build a method rooted in a decomposition strategy to solve the inverse problems numerically.

Numerical experiments demonstrate that the order of convergence of the numerical method for solving the direct problem is second in space and first in time. Since the numerical methods for solving inverse problems recover almost exactly the numerical solution of the direct problem and the boundary functions when provided with exact observations, we conclude that the numerical algorithms for addressing inverse problems achieve the same order of accuracy as the numerical method used for the direct problem.

Utilizing natural boundary conditions for numerically solving inverse problems results in a more effective recovery compared to algorithms that employ Dirichlet boundary conditions at the degenerate boundary $r=0$. The Dirichlet boundary condition at the right boundary is recovered with greater precision than that at the left boundary. Moreover, the recovery of the Dirichlet boundary condition at the degenerate boundary is significantly influenced by the location of the measurement(s).

Computational results with noisy data demonstrate that the numerical method is able to recover Dirichlet boundary conditions and the numerical solution of the direct problem with optimal accuracy, especially when the natural boundary condition is imposed on the degenerate boundary.

The limitation of the proposed numerical approach for solving the inverse problem lies in the need to control the measurement locations by selecting appropriate positions or adjusting the time and space steps to prevent excessively small values of $v_{r^{*}}$ and $w_{r^{*}}$.

It is of interest to apply a collocation method to heat conduction problems with spherical symmetry, using a meshless localized space–time radial basis function approach, as proposed in [47], for simulating transient heat conduction in functionally graded materials with variable material properties.