Simultaneous Numerical Determination of Two Time-dependent Coefficients in Second Order Parabolic Equation With Nonlocal Initial and Boundary Conditions

TL;DR

This paper presents a method to accurately determine two time-dependent coefficients in a mathematical model using a stable numerical approach and regularization.

Contribution

A novel framework for simultaneously identifying two time-dependent coefficients in a parabolic equation with nonlocal conditions using regularization and optimization.

Findings

The Crank-Nicolson FDM ensures stability and accuracy in solving the direct problem.

Tikhonov regularization effectively reduces errors and stabilizes the inverse problem solution.

Numerical experiments confirm the method's reliability under exact and noisy data.

Abstract

This study establishes a mathematically consistent and computational framework for the simultaneous identification of two time-dependent coefficients in a one-dimensional second-order parabolic partial differential equation. The considered problem is governed by nonlocal initial, boundary, and integral overdetermination conditions.

The direct problem is solved using the Crank-Nicolson finite difference method (FDM), which ensures unconditional stability and second-order accuracy in both spatial and temporal discretizations. The corresponding inverse problem is reformulated as a nonlinear regularized least-squares optimization problem and efficiently solved used the MATLAB subroutine
lsqnonlin from the optimization Toolbox. Due to the intrinsic, ill-posedness of the inverse formulation, small input data errors lead to big output errors. Then, Tikhonov regularization, is employed to enhance numerical stability and robustness.

Extensive numerical experiments are carried out under exact and noisy data to evaluate the numerical accuracy and convergence behavior of the method. The results confirm that the regularization technique effectively damps numerical oscillations, minimizes reconstruction error, and ensures reliable recovery of the unknown coefficients. Sensitivity analysis further reveals the essential role of the regularization parameter in controlling the trade-off between stability and accuracy.

The proposed approach provides an accurate and computationally efficient tool for IP in heat transfer, diffusion processes, and related applied sciences.

Figures

The graphs showing analytical and computational distribution with absolute error for the direct problem (1)−(4) , for example1 (a1)−(d1), represent the size of mesh N=M={20,40,80,160} , respectively.

The graphs showing the exact and numerical values for, required output w(v1,τ) equation (5) , for example 1 (a2)−(d2), with the size of mesh N=M={20,40,80,160} , respectively.

The graphs showing the exact and numerical values for, required output w(v2,τ) equation (5) , for example 1 (a3)−(d3), with the size of mesh N=M={20,40,80,160} , respectively.

The graphs showing analytical and computational distribution with absolute error for the direct problem (1)−(4) , for example2 (a4)−(d4), represent the size of mesh N=M={20,40,80,160} , respectively.

The graphs showing the exact and numerical values for, required output w(v1,τ) equation (5) , for example 2 (a5)−(d5), with the size of mesh N=M={20,40,80,160} , respectively.

The graphs showing the exact and numerical values for, required output w(v2,τ) equation (5) , for example 2 (a6)−(d6), with the size of mesh N=M={20,40,80,160} , respectively.

The graph showing the objective function for the IP (1)−(5) , in example 3 when the size of mesh N=M=80 .

The graphs (a8) and (b8) showing the reconstructed coefficient a(τ) and b(τ) in comparison with exact value for the IP (1)−(5) , in above example when the size of mesh N=M=80 .

The graph showing the objective function when p=0.01%, no regularization for the IP (1)−(5) , in example 2, when the size of mesh N=M=80 .

The graphs (a10) and (b10) showing the exact and (unstable) numerical solution for a(τ) and b(τ) when p=0.01%, no regularization for the IP (1)−(5) , in example 3, when the size of mesh N=M=80 .

The graph showing the objective function when p=0.01% noise with regularization, β∈{10−3,…,10−7} , for the IP (1)−(5) , in example 3, when the size of mesh N=M=80 .

The graphs (a12) and (b12) showing the exact and (stable) numerical solution for a(τ) and b(τ) when p=0.01%, and regularization, β∈{10−3,…,10−7} , for the IP (1)−(5) , in example 3, when the size of mesh N=M=80 .

Tables

Table 1.. Numerical information for no regularization β=0 and noise p=0 .

M=N	80
RMSE a( τ )	0.0241
RMSE b( τ )	2.1592 ×10−4
Computational time	3594 Sec
No. of Iterations	16
Obj. Function	1.2221× 10−10

Table 2.. Numerical information for various regularization parameters and noise p=0.01%,(N=80 ).

β	10−3	10−4	10−5	10−6	10−7	10−8	10−9	10−10
RMSE a( τ )	0.3302	0.1220	0.0714	0.0617	0.1189	0.2408	0.2835	0.2865
RMSE b( τ )	0.0552	0.0225	0.0259	0.0344	0.0365	0.0373	0.0377	0.0377
Computational time	84320 Sec	9211 Sec	1453 Sec	1451 Sec	3139 Sec	6277 Sec	5814 Sec	5506 Sec
No. of Iterations	401	43	6	6	14	22	22	22

Full text

1. Introduction

Inverse problems (IPs), arise naturally in many scientific and engineering disciplines as they are typically modeled by differential equations. In the context of ordinary differential equations, the direct problem refers to obtaining solutions for a given system, whereas IP focuses on reconstructing the governing system from observed characteristics. Historically, the role of IPs was recognized in celestial mechanics. ^ 1 ^ In practical applications, IPs are commonly encountered when one seeks to determine unknown causes from their observed outcomes, in contrast to direct problems where effects are predicted from known causes. Compared with direct problems, IPs, are often more challenging to solve due to their ill-posedness, solutions may fail to exist, may not be unique, or may not change continuously with the input data. ^ 2 ^ Such problems appear in almost every scientific and technological domain, particularly in models derived from social and physical systems. Since most of these models are expressed through differential and integral equations, their analysis requires not only solving the equations but also interpreting the system’s behavior under varying conditions, provided that sufficient information is available to fully describe the model. IPs, associated with such equations arise in a wide range of scientific and engineering applications, as well as in the modeling of social processes. Many physical phenomena are described by mathematical models in the form of initial and boundary value problems for partial differential equations (PDEs). Frequently, these models involve differential and integral equations. ^ 3 ^ In boundary type IPs, the goal is to recover unknown boundary conditions, which usually leads to classical ill-posedness because the absence of continuous dependence of the solution on the input data. Numerical solution (NS) of direct mathematical physics problems is presently a well-studied matter. In solving multi-dimensional boundary value problems, difference methods and the finite element method are widely used. Partial differential equations give a basis for applied mathematical models, and their solution can be found in mathematical physics equations and some relations, including additional relations, boundary, and initial conditions, among other elements. ^ 4 ^ The following reconstruction problems serve as examples of IP applications in daily life, tomography, initial condition estimation of transient problems in conductive heat transfer, detection of non-metallic materials beneath the surface using reflected radiation methods, intensity, position estimation of illuminated radiation from a biological source using experimental radiation measures, and thermal source intensity estimation with functional dependence in space. ^ 5– 9 ^ On the other side, there are three different kinds of differential equations, Second-order equations are the most crucial for applications, Elliptic, Parabolic, and hyperbolic equations are examples of these equations. Yuldashev and other authors studied elliptic type integrate differential equations in ^ 10, 11 ^ while hyperbolic and parabolic were studied in ^ 12– 14 ^ and ^ 15– 17 ^ respectively. The alternating direction explicit method used for reconstructing solutions is also efficient and unconditionally stable, in ^ 18 ^ is restructured into a nonlinear regularized least-square optimization problem, and is effectively resolved using the MATLAB subroutine lsqnonlin from the optimization toolbox. The Crank-Nicolson (C-N) FDM together with the TR was used in ^ 19– 21 ^ the effectiveness of the computational method was shown together with the proof of the solution’s existence and uniqueness. To find a stable and accurate approximate solution of finite differences. In ^ 22– 26 ^ the study shows that for the IP with the time-dependent source coefficient in heat equations, there’s a smooth solution pair as extended as a single data point at an observation point and by providing an explicit formula for the time-dependent source coefficient, additional information at the observation point improves our comprehension of the behavior of the problem. This paper’s IP classical solution for a second order parabolic equation has been shown to be exist and unique by Elvin Azizbayov. ^ 27 ^ As a result, the main objective of the current work is a numerical realization of such a problem. In this paper, the article is divided into five sections, Section 2, includes a mathematical formulation of the IP under study. Section 3, a numerical method for solving the forward problem using C-N FDM. In Section 4, the IPs mathematical solution is shown with an initial guess, while, in Section 5, discuss the numerical problem and explain the obtained results.

2. Mathematical formulation

As a mathematical model, we consider $[eqn]$ be a rectangular region that is defined by $[eqn]$ and $[eqn]$ be a fixed number, the one-dimensional IP of determining of unknown functions $[eqn]$ for the upcoming parabolic equation,

$$[eqn]$$

with, the nonlocal initial conditions (ICs)

$$[eqn]$$

with, periodic boundary condition,

$$[eqn]$$

and, nonlocal integral condition

$$[eqn]$$

and the overdetermination conditions

$$[eqn]$$

where, $[eqn]$ , are fixed numbers, $[eqn]$ are given functions, $[eqn]$ are the sought functions ( $[eqn]$ is space component and $[eqn]$ is time component).

2.1 Definition 1

Consider the triplet $[eqn]$ represents a classical solution, to the IP $[eqn]$ if the functions $[eqn]$ and $[eqn]$ satisfies the following conditions:

• 1)The function $[eqn]$ and its derivatives $[eqn]$ are continuous in the domain $[eqn]$
• 2)The functions $[eqn]$ and $[eqn]$ are continuous on $[eqn]$

Equation (1) and conditions $[eqn]$ are satisfied in the classical (usual) sense.

2.2 Existence of the inverse problem classical solution

2.2.1 Theorem 1[28]: Let the assumptions $[eqn]$ and the condition

$$[eqn]$$

where,

$$[eqn]$$

where, $[eqn]$ and $[eqn]$ ; $[eqn]$

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where,

$$[eqn]$$

holds, then problem $[eqn]$ , has a unique solution in the ball $[eqn]$ of the space $[eqn]$ , where, $[eqn]$

3. FDM scheme for direct (Forward) problem

The direct problem will be resolved in this section, i.e., when unknown, $[eqn]$ , $[eqn]$ and $[eqn]$ are assumed to be given, in order, to solve this problem, Crank-Nicolson (C-N) FDM scheme had been used for finding the solutions of the nonlocal problem given by Equation (1)–(4) $[eqn]$ Where the domain $[eqn]$ had been divided into $[eqn]$ mesh with spatial step size $[eqn]$ and the time step size $[eqn]$ where $[eqn]$ and $[eqn]$ are given positive integers. The grid points have been given by

$$[eqn]$$ $$[eqn]$$

These quantities’ discretized form is given as follows,

$$[eqn]$$

Applying the C-N FDM, for the discretizes Equation (1) $[eqn]$ which approximated as,

$$[eqn]$$

where, $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and, $[eqn]$ simplifying the above Equation (28), get the following difference equation,

$$[eqn]$$

where, $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and, $[eqn]$ .

The C-N FDM discretizes Equations (2)–(4) which approximated as,

$$[eqn]$$

where $[eqn]$ is integral in Equation (2), and we have,

$$[eqn]$$

Finally, the trapezoidal rule and C-N discretizes integral condition ( $[eqn]$ ) as;

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

From Equation (29) $[eqn]$ with the nonlocal initial conditions $[eqn]$ (assuming $[eqn]$ ), and using Equations (30)–(33), with each time step $[eqn]$ , we can be described in matrix form as

$$[eqn]$$ $$[eqn]$$

with,

$$[eqn]$$

and, $[eqn]$ , where,

$$[eqn]$$ $$[eqn]$$

the Equation (4) desired output data will be approximated using trapezoidal rule as

$$[eqn]$$

3.1 Stability and convergence analysis

The adopted C-N finite difference scheme used to solve the direct problem is unconditionally stable and second order accurate in both spatial and temporal discretization. A brief von Neumann stability analysis shows that the amplification factor satisfies $[eqn]$ , for all mesh ratios, which guarantees the unconditional stability of the scheme. The local truncation error is of order $[eqn]$ confirming that the method provides second order accuracy and reliable convergence toward the analytical solution as the mesh is refined. Hence, the C-N formulation ensures both numerical stability and accuracy for the present problem.

3.2 Examples for direct problem

3.2.1 Example 1: Assume an illustration for the direct problem given by $[eqn]$ with $[eqn]$ and the input data are,

$$[eqn]$$

and, the analytical solution is,

$$[eqn]$$

We investigate the accuracy of the solution using diverse mesh grid size, $[eqn]$ . In Figure 1 (a)–(d), an excellent agreement between the exact and numerical solutions can be clearly observed, indicating high accuracy. It can be noted that as the number of mesh points increases, the accuracy of the obtained solution improves, which reflects the convergence and stability of the proposed numerical scheme. Moreover, the absolute error graph at each mesh shows that the maximum magnitude of error didn’t exceed $[eqn]$ , as depicted in Figure 1 (a)–(d).

*The graphs showing analytical and computational distribution with absolute error for the direct problem

(1)−(4) , for example1

(a1)−(d1), represent the size of mesh

N=M={20,40,80,160} , respectively.*

The obtained results clearly indicate that the mesh size has a noticeable influence on the numerical accuracy. As the spatial and temporal mesh are refined (i.e., as $[eqn]$ and $[eqn]$ increase), the numerical solution becomes smoother and shows a closer agreement with the analytical one. This behavior confirms the expected second-order convergence of the C-N finite difference scheme in both time and space. Beyond a certain refinement level (e.g., $[eqn]$ ), further mesh subdivision produces only negligible changes, implying that the proposed schemeis numerically stable and nearly mesh-independent. Whilst, Figure 2 and Figure 3, show the numerical outcomes for the required information $[eqn]$ evaluated at various mesh sizes such that (s.t.), observed that an excellent agreement is obtained.

*The graphs showing the exact and numerical values for, required output

w(v1,τ) equation

(5) , for example 1

(a2)−(d2), with the size of mesh

N=M={20,40,80,160} , respectively.*

*The graphs showing the exact and numerical values for, required output

w(v2,τ) equation

(5) , for example 1

(a3)−(d3), with the size of mesh

N=M={20,40,80,160} , respectively.*

3.2.2. Example 2: Assume an illustration for the direct problem given by $[eqn]$ with $[eqn]$ and the input data are,

$$[eqn]$$

and, the analytical solution is,

$$[eqn]$$

The accuracy of the numerical solution is examined using several mesh grid size, $[eqn]$ . Figure 4 $[eqn]$ , clearly demonstrates an excellent agreement between the analytical and numerical solutions, confirming the reliability of the proposed approach. It can be noted that as the number mesh becomes finer, the numerical accuracy improves, indicating the expected convergence and stability of the C-N finite difference scheme. Furthermore, the absolute error plots at different mesh sizes reveal that the maximum error magnitude remains below $[eqn]$ , as illustrated in Figure 4 $[eqn]$ .

*The graphs showing analytical and computational distribution with absolute error for the direct problem

(1)−(4) , for example2

(a4)−(d4), represent the size of mesh

N=M={20,40,80,160} , respectively.*

Similarly, in Example 2, the refinement of the spatial and temporal mesh has a direct impact on the numerical accuracy. As the grid is successively refined, the computed solution aligns more closely with the analytical one, and the absolute error distribution becomes smoother across the domain. This confirms that the proposed C-N finite difference formulation maintains its second-order convergence and numerical stability. It can therefore be concluded that the results are consistent and nearly insensitive to further mesh refinement beyond the tested resolutions, indicating satisfactory mesh-independence of the method. Whilst, Figure 5 and Figure 6, show the numerical outcomes for the required information $[eqn]$ evaluated at various mesh sizes s.t. observed that an excellent agreement is obtained.

*The graphs showing the exact and numerical values for, required output

w(v1,τ) equation

(5) , for example 2

(a5)−(d5), with the size of mesh

N=M={20,40,80,160} , respectively.*

*The graphs showing the exact and numerical values for, required output

w(v2,τ) equation

(5) , for example 2

(a6)−(d6), with the size of mesh

N=M={20,40,80,160} , respectively.*

4. Inverse problem solution for

Equations (1)–(5)

For the nonlinear IP $[eqn]$ we seek precise and stable identification of $[eqn]$ and $[eqn]$ , that mean potential term $[eqn]$ , and $[eqn]$ is unknown. The one-dimensional second-order parabolic equation together with $[eqn]$ satisfies the problem given by Equation (1)–(5), the problem is reformulated as a nonlinear least-squares minimization task. The resulting minimization problem is solved numerically using the lsqnonlin routine available in MATLAB’s optimization Toolbox, which implements a trust-region reflective algorithm suitable for nonlinear least squares problems with bound constraints. Due to the ill-posedness of such problems, especially in the presence of noise data, Tikhonov Regularization (TR) is employed to stabilize the solution. The regularization parameter $[eqn]$ , is selected carefully to balance the trade-off between measured and the numerically computed solution. For this purpose, MATLAB routine lsqnonlin is employed, to minimize the Tikhonov Regularized objective functional, the procedure begins with a suitably selected initial guess. The TR functional is derived based on condition $[eqn]$ , and the functional error is incorporated as follows:

$$[eqn]$$

where $[eqn]$ is a regularization parameter. The discretization form of $[eqn]$ is

$$[eqn]$$

The unregularized case, i.e., $[eqn]$ produces the regular a nonlinear least-squares functional, which is inherently unstable when dealing with noisy data. The MATLAB routine lsqnonlins is used to minimize $[eqn]$ under some physical constraint. The following parameters had been used for the subroutine:

• •(Maxlter) maximum number of iterations $[eqn]$ .
• •Solution and Objective function tolerance $[eqn]$
• •The lower and upper bounds on the component of the vector $[eqn]$ are $[eqn]$ and $[eqn]$ , respectively.

The IP $[eqn]$ are resolved via both precise and noisy measurement $[eqn]$ . By including a random error, the noisy data is numerically simulated:

$$[eqn]$$

where $[eqn]$ , is random Gaussian normal distribution vectors with mean zero and standard deviations $[eqn]$ , given by

$$[eqn]$$

where, $[eqn]$ is the percentage of noise. We use the MATLAB function “ normrnd”, as:

$$[eqn]$$

to generate the random variables $[eqn]$ and, $[eqn]$

5. Results and discussion

To evaluate the precision and stability of the numerical methods, a set of numerical experiments is conducted. These experiments are designed to evaluate the accuracy of the computed results by simulating realistic measurement conditions, where noise is introduced into the input data. To quantitatively evaluate the discrepancy between the exact and the numerically computed solutions, the root mean squares errors (RMSE) are utilized by the following expression

$$[eqn]$$

we take, $[eqn]$ for simplicity.

5.1. Example 3: IP

Consider the IP $[eqn]$ with the following input data,

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ , and, the analytical solution for this IP is provided as follows,

$$[eqn]$$

The time-dependent coefficient $[eqn]$ , and $[eqn]$ , are then reconstructed (in case, $[eqn]$ , is taken), and consider the case of noise-free in the measurements, see Table 1, this mean $[eqn]$ in $[eqn]$ and, no regularization applied, Figure 7 shows the objective function $[eqn]$ without regularization i.e. $[eqn]$ while, in Figure 9, noise without regularization.

**Table 1.: Numerical information for no regularization

β=0 and noise

p=0 .**

*The graph showing the objective function for the IP

(1)−(5) , in example 3 when the size of mesh

N=M=80 .*

*The graphs

(a8) and

(b8) showing the reconstructed coefficient

a(τ) and

b(τ) in comparison with exact value for the IP

(1)−(5) , in above example when the size of mesh

N=M=80 .*

*The graph showing the objective function when

p=0.01%, no regularization for the IP

(1)−(5) , in example 2, when the size of mesh

N=M=80 .*

We investigate the time-dependent coefficient $[eqn]$ and $[eqn]$ under both exact and perturbed measurements, as demonstrated in Figures 7 $[eqn]$ 12. The associated results are unstable, as illustrated in Figure 9- Figure 10. This reveals that the problem is not properly posed, a noticeable instability appears in the results, which is anticipated due to the ill-posed nature of the investigated IP. The unregularized solutions exhibit noticeable instability even in the noise-free case, demonstrating the sensitivity of the problem to small perturbations in the input data. The TR approach is used by incorporating the penalty term $[eqn]$ into the classical least-squares formulation, as presented in Equation (37). In Figure 11, $[eqn]$ noise strategy is employed to enhance the solution’s robustness with regularization parameters, $[eqn]$ , yielding numerically stable regularized reconstructions of $[eqn]$ and $[eqn]$ , as shown in Table 2. These results indicate that smaller values of $[eqn]$ are effective when the data is not contaminated.

*The graphs

(a10) and

(b10) showing the exact and (unstable) numerical solution for

a(τ) and

b(τ) when

p=0.01%, no regularization for the IP

(1)−(5) , in example 3, when the size of mesh

N=M=80 .*

*The graph showing the objective function when

p=0.01% noise with regularization,

β∈{10−3,…,10−7} , for the IP

(1)−(5) , in example 3, when the size of mesh

N=M=80 .*

*The graphs

(a12) and

(b12) showing the exact and (stable) numerical solution for

a(τ) and

b(τ) when

p=0.01%, and regularization,

β∈{10−3,…,10−7} , for the IP

(1)−(5) , in example 3, when the size of mesh

N=M=80 .*

**Table 2.: Numerical information for various regularization parameters and noise

p=0.01%,(N=80 ).**

Next, we will study the proportion of noise contamination with a ratio of $[eqn]$ noise, in $[eqn]$ via $[eqn]$ , the associated results are unstable, as illustrated. This reveals that the problem is not properly, posed and small error in the data input cause large errors in the output solution, thus, the problem needs to be regularized. As a consequence, for the situation of $[eqn]$ noise, the best selection for the regularization parameters is $[eqn]$ for $[eqn]$ and $[eqn]$ for $[eqn]$ as in Table 2, which results in the lowest RMSE values.

6. Conclusions

In this study, the identification of the unknown time-dependent coefficient $[eqn]$ and $[eqn]$ in a one-dimensional parabolic PDE are investigated. The problem is formulated from a numerically evaluated nonlocal integral subjected to initial and boundary conditions. The direct problem is solved using the Crank-Nicolson (C-N) scheme, while the IP is transformed into a nonlinear least squares optimization problem using the lsqnonlin routine in MATLAB. An example test cases are employed in numerical experiments to evaluate the accuracy and stability of the proposed method. According to the results, using TR significantly enhances the solution’s stability. When the root-mean square error (RMSE) numbers are their lowest, the most effective approach is identified. Furthermore, numerical results have been provided to illustrate the accuracy and stability of the numerical method. It found that applying TR method stabilizes the solution. Cases of numerical results exist with noise (with and without) regularization. In the first case, consider when $[eqn]$ noise and without regularization as in Figure 9- Figure 10 in example 3, the results are unstable and highly oscillated, on the contrary. Other case, when $[eqn]$ noise introduced with $[eqn]$ regularization, obtain are stable and accurate results for the reconstructed coefficient $[eqn]$ and $[eqn]$ , which are (stable) numerical solutions. As for the remaining cases, we observe that the results are stable and steady when regularization values are $[eqn]$ for $[eqn]$ and, $[eqn]$ for $[eqn]$ in example 3, as shown in Table 2.