The inverse problem of reconstructing the time-dependent thermal conductivity and free boundary coefficients along with the temperature in a two-dimensional parabolic equation with initial and boundary conditions and additional measurements is, for the first time, numerically investigated. This inverse problem appears extensively in the modelling of various phenomena in engineering and physics. For instance, steel annealing, vacuum-arc welding, fusion welding, continuous casting, metallurgy, aircraft, oil and gas production during drilling and operation of wells. From literature we already know that this inverse problem has a unique solution. However, the problem is still ill-posed by being unstable to noise in the input data. For the numerical realization, we apply the alternating direction explicit method along with the Tikhonov regularization to find a stable and accurate numerical solution of finite differences. The root mean square error (rmse) values for various noise levels p for both smooth and non-smooth continuous time-dependent coefficients Examples are compared. The resulting nonlinear minimization problem is solved numerically using the MATLAB subroutine lsqnonlin. Both exact and numerically simulated noisy input data are inverted. Numerical results presented for two examples show the efficiency of the computational method and the accuracy and stability of the numerical solution even in the presence of noise in the input data.
Free boundary problems for the parabolic partial differential equations have significant applications in various fields of engineering, physics, chemistry, see [1–5] to mention only a few. A Stefan problem is a moving boundary value problem that concerns the distribution of heat in a period of transforming medium. The authors of [6] considered two fairly different free boundary problems of nonlinear diffusion. Chen et al. [7] recast the formulation of various shock diffraction/reflection problems as a free boundary problem. Shidfar et al. [8] studied the inverse moving boundary problem using the least-squares method. In [9], the authors investigated the numerical solution of inverse Stefan problems using the method of fundamental solutions. Recently, the authors of [10] numerically solved the inverse heat problems of determining the time-dependent coefficients with free boundaries.
There is also an analysis of a one-dimensional inverse problem of reconstructing the timewise heat source with a moving boundary [11]. The authors in [12] estimated free boundary coming from two new scenarios, aggregation processes and nonlocal diffusion. Snitko [13,14], theoretically, and Huntul [15], numerically, investigated the inverse problem of determining the time-dependent reaction coefficient in a two-dimensional parabolic problem with a free boundary.
The challenge associated with free boundary problems arises from the finding that the solution domain is unknown. Only a few studies focus on the time-dependent free boundary in higher dimensions [16–19]. These papers are theoretical and very important because they present conditions for the unique solvability of the unknown coefficients.
This work examines the inverse problem to recover the time-dependent thermal conductivity coefficient and free boundaries from the heat flux and the nonlocal integral observations as over-specification conditions. The inverse problem presented in this paper has already been showed to be locally uniquely solvable by Barans’ka et al. [20], but no numerical determination has been realised so far, therefore, the main goal of this work is to attempt numerical solution of this problem.
The arrangement of this paper is systematized as follows. Section 2 describes the formulations of the inverse problem. The solution of the direct problem using the alternating direction explicit (ADE) method is presented in Section 3. The ADE direct solver is coupled with the Tikhonov regularization method in Section 4. In Section 5, computational results and discussions are presented. Finally, Section 6 highlights the conclusions.
Formulation of the Inverse Problem
Consider the inverse problem of reconstructing time-dependent thermal conductivity coefficient a(t) > 0, and the free boundaries l(t) > 0 and h(t) > 0 in the two-dimensional parabolic equation
where f(x1,x2,t) is a known heat source, u=u(x1,x2,t) is the unknown temperature in the moving domain ΩT:={(x1,x2,t)∣0<x1<l(t),0<x2<h(t),0<t<T<∞}, with the initial condition
u(x1,x2,0)=φ(x1,x2),(x1,x2)∈[0,l0]×[0,h0],
where l0 = l(0) and h0 = h(0) are given positive numbers, the Dirichlet boundary conditions
where φ and μi for i=1,7¯ are given functions. The function μ5(t) in Eq. (5) represents heat flux boundary condition. The function μ6(t) in Eq. (6) corresponds to the specification of mass/energy, [21,22], whilst the function μ7(t) in Eq. (7) represents the first-order moment specification, [23].
Introducing the new variables y1 = x1/l(t) and y2 = x2/h(t), see [20], we recast the problem given by Eqs. (1) to (7) into the problem given below for the unknowns l(t), h(t), a(t) and v(y1,y2,t):=u(y1l(t),y2h(t),t), namely,
Then, it is possible to indicate a time T0∈(0,T], determined by the input data, such that there exists a solution (l(t),h(t),a(t),v(y1,y2,t))∈C1[0,T0]×C1[0,T0]×C[0,T0]×C2,1(Ω¯T0) with l(t) > 0, h(t) > 0, a(t) > 0 for t∈[0,T0], to problem Eqs. (8) to (14).
Theorem 2Suppose that the following conditions are fulfilled:
Then, it is possible to indicate a time T1∈(0,T], determined by the input data, such that problem Eqs. (8) to (14) has at most one solution (l(t),h(t),a(t),v(y1,y2,t))∈C1[0,T1]×C1[0,T1]×C[0,T1]×C2,1(Ω¯T1) with l(t) > 0, h(t) > 0, a(t) > 0 for t∈[0,T1].
Numerical Solution of the Forward Problem
Consider the direct (forward) problem Eqs. (8) to (11). When l(t), h(t), a(t), f(x1,x2,t), μi(t), for i=1,4¯ and φ(y1l0,y2h0) are given in the direct problem v(y1,y2,t) is to be found along with the quantities of interest μi(t) for i=5,7¯. Denote v(y1,y2,tn)=vi,jn, l(tn) = ln, h(tn) = hn, a(tn) = an and f(y1l(tn),y2h(tn),tn)=fi,jn, where y1i=iΔy1, y2j=jΔy2, tn=nΔt, Δy1=1/M1, Δy2=1/M2, Δt=T/N, i=0,1,…,M1, j=0,1,…,M2, n=0,1,…,N.
The ADE method, [3,24], which is unconditionally stable, is described for solving numerically the direct problem Eqs. (8) to (11). Let ṽi,jn and ũi,jn satisfy
In this section, our aim is to obtain stable and accurate reconstructions of the free boundary l(t), h(t), the thermal conductivity a(t) together with the temperature v(y1,y2,t), satisfying Eqs. (8) to (14). The inverse problem can be formulated as a nonlinear least-squares minimization of the least-squares objective function given by
where v(y1,y2,t) solves Eqs. (8) to (11) for given (l(t),h(t),a(t)), respectively. The minimization of the function Eq. (27) is performed using the MATLAB subroutine lsqnonlin [25].
The inverse problem given by Eqs. (8) to (14) is solved subject to both analytical and noisy (perturbed) measurements Eqs. (12) to (14). The perturbed data are numerically formulated as follows
μiϵ(tn)=μi(tn)+ϵn,n=1,N¯,i=5,7¯,
where ϵn are random variables with mean zero and standard deviations σi, given by
σi=p× maxt∈[0,T]|μi(t)|,i=5,7¯,
where p represents the percentage of noise. We use the MATLAB function normrnd to generate the random variables ϵ=(ϵn)n=1,N¯ as follows:
ϵ=normrnd(0,σi,N),i=1,2,3.
In the case of perturbed data Eq. (28), we replace in Eq. (27)μi(tn) by μiϵ(tn) for i=5,7¯.
Numerical Results and Discussion
The accuracy is measured by rmse:
rmse(l)=[TN∑n=1N(lnumerical(tn)-lexact(tn))2]1/2,
rmse(h)=[TN∑n=1N(hnumerical(tn)-hexact(tn))2]1/2,
rmse(a)=[TN∑n=1N(anumerical(tn)-aexact(tn))2]1/2.
Let us take T = 1, for simplicity. The lower bounds and upper bounds for l(t) > 0, h(t) > 0 and a(t) > 0 are 10−9 for (LB) and 102 for (UB), respectively.
Example 1
Consider the inverse problem given by Eqs. (1) to (7), with the smooth unknown timewise terms l(t), h(t) and a(t), and we solve it with:
We observe that the conditions (A1)–(A6) of Theorems 1 and 2 are fulfilled and thus, the uniqueness condition of the solution is guaranteed. It can be easily verified that the analytical solution of Eqs. (1) to (4) is
First, let us solve the forward problem Eqs. (1) to (4) with the input data (34) when l(t), h(t), a(t) are known and given by Eq. (37). Tab. 1 reveals that the exact and approximate solutions for Eqs. (12) to (14), which exactly is given by Eq. (35), obtained with number of grids M1 = M2 = 10 and N∈{20,40,80}, are in good agreement. The analytical Eq. (38) and approximate solutions for v(y1,y2,t) is plotted in Fig. 1. It is clear from Tab. 1 and Fig. 1 that the accuracy of the approximate solution increases, as the mesh sizes decreases.
The numerical and analytical Eq. (35) solutions for μi(t), i=5,7¯, with M1 = M2 = 10 and various N∈{20,40,80}, for direct problem
t
0.1
0.2
0.3
…
0.8
0.9
N
rmse (μi)
μ5(t)
1.9011
1.8021
1.7032
…
1.2043
1.1059
20
1.3E −3
1.9001
1.8002
1.7003
…
1.2002
1.1001
40
9.3E −4
1.9000
1.8000
1.7000
…
1.2000
1.1000
80
3.6E −5
1.9000
1.8000
1.7000
…
1.2000
1.1000
exact
0
μ6(t)
2.2571
2.7787
3.3689
…
7.5195
8.6382
20
5.8E −3
2.2589
2.7800
3.3783
…
7.5147
8.6318
40
2.3E −3
2.2596
2.7810
3.3781
…
7.5128
8.6294
80
8.1E −4
2.2611
2.7817
3.3700
…
7.5125
8.6280
exact
0
μ7(t)
1.0712
1.4233
1.8530
…
5.6705
6.9071
20
4.4E −3
1.0710
1.4231
1.8525
…
5.6651
6.9001
40
7.1E −3
1.0711
1.4232
1.8527
…
5.6728
6.9023
80
9.3E −4
1.0746
1.4278
1.8583
…
5.6738
6.9079
Exact
0
The solutions v(y1,y2,1) and absolute errors with M1 = M2 = 10 for: (a) N = 20, (b) N = 40 and (c) N = 80, for direct problem
In the problem Eqs. (8) to (14), we take the initial guesses for l,h and a, as follows:
l0(tn)=l0=1,h0(tn)=h0=1,a0(tn)=a(0)=2,n=1,N¯.
We take a mesh size with M1 = M2 = 10, N = 40 and we start the examination for recovering the timewise terms l(t), h(t) and a(t), when p = 0 in the measurements Eqs. (12) to (14), as in Eq. (31). The function Eq. (27) is illustrated in Fig. 2a, where a convergence, which is monotonically decreasing, is obtained in 10 iterations to achieve a low specified tolerance of O(10−28). Figs. 2b–2d illustrates the corresponding approximate results for the functions l(t), h(t) and a(t). A good agreement between the analytical Eq. (37) and approximate solutions can be observed from these figures with rmse(l) = 4.4E −3, rmse(h) = 3.8E −3 and rmse(a) = 6.4E −3.
(a) The objective function Eq. (27)vs. the number of iterations, and the exact Eq. (37) and numerical solutions for: (b) l(t), (c) h(t) and (d) a(t), with no noise, for Example 1
Next, the stability of the numerical solution is investigated with respect to the noisy data Eq. (28). We include p∈{0.1%,0.5%,1%,2%} noise to the over-determination conditions μ5, μ6 and μ7, as in Eq. (28) to test the stability. The function Eq. (27) is plotted in Fig. 3a and convergence is again rapidly achieved. The approximate results for l(t), h(t) and a(t) are depicted in Fig. 3. Accurate and stable results are obtained for l(t), h(t) and a(t) by observing Figs. 3b–3d and the rmse values, the minimum value of Eq. (27) at final iteration and the number of iterations are given in Tab. 2. As noise p is increased the approximate results for l(t), h(t) and a(t) start to build up oscillations. In order to retrieve stability, we penalise the function Eq. (27) by adding λ(∥l′(t)∥2+∥h′(t)∥2+∥a(t)∥2) to it since the theory provide l∈C1[0,T], h∈C1[0,T] and a∈C[0,T], where 0$]]>λ>0 is the Tikhonov’s regularization parameter to be chosen. Then, in discretised form of Tikhonov functional recasts as
(a) The objective function Eq. (27)vs. the number of iterations, and the exact Eq. (37) and numerical solutions for: (b) l(t), (c) h(t) and (d) a(t), with p∈{0.1%,0.5%,1%,2%} noise, for Example 1
The values of rmseEqs. (31) to (33), the function Eq. (27) at final iteration and the number of iterations, for p∈{0,0.1%,0.5%,1%,2%} noise, for Example 1
p (%)
rmse(l)
rmse(h)
rmse(a)
Minimum value of (27)
Iter
0
0.0044
0.0038
0.0064
2.5E −28
10
0.1
0.0093
0.0079
0.0122
3.2E −28
10
0.5
0.0427
0.0358
0.0571
3.8E −28
10
1
0.0847
0.0717
0.1141
3.0E −28
21
2
0.1698
0.1448
0.1141
2.3E −28
21
For p=5% noise, Fig. 4 illustrates the analytical solution Eq. (37) and the approximate solutions obtained by minimizing the functional Eq. (40) for various regularization parameters. The rmse(l,h,a) values are {0.4822,0.4313,0.3049,0.1574,0.0566,0.1040}, {0.4388,0.3720,0.2857,0.1575,0.0646,0.1039} and {0.6357,0.5994,0.4222,0.2375,0.1503,0.2363} for λ∈{0,10-5,10-4,10-3,10-2,10-1}, see Tab. 3 for more information. It can be observed that the approximate unregularized solution obtained with λ=0 demonstrates instability, however, on including regularization with λ=10-3 to λ=10-2, we get a stable solution which is consistent in accuracy with the p=5% noise desecrating the input data Eq. (28).
The numerical and analytical Eq. (37) solutions for: (a) l(t), (b) h(t) and (c) a(t), for p=5% noise, with various λ∈{0,10-2,10-1}, for Example 1
The values of rmseEqs. (31) to (33) and the number of iterations, for p=5% noise, without and with regularization, for Example 1
p (%)
λ
rmse(l)
rmse(h)
rmse(a)
Iter
5
0
0.4822
0.4388
0.6357
50
5
10−5
0.4313
0.3720
0.5994
10
5
10−4
0.3049
0.2857
0.4222
10
5
10−3
0.1574
0.1575
0.2375
10
5
10−2
0.0566
0.0646
0.1503
10
5
10−1
0.1040
0.1039
0.2363
10
Example 2
The smooth timewise coefficients l(t), h(t), a(t) given by Eq. (37) has been recovered in previous example. In this example, let us examine the numerical scheme for recovering a non-smooth test:
l(t)=11+t,h(t)=11+t,a(t)=1+ cos2(2πt),t∈[0,1]
and the analytical solution for the temperature u(x1,x2,t) given by Eq. (36). The rest of the input data are given as
With this data, the conditions (A1)–(A6) of Theorems 1 and 2 are holds, thus the uniqueness of the solution is also ensured. The initial guesses for l,h and a, for this example has been taken as in Eq. (39).
As in Example 1, we take Δy1=Δy2=0.1 and Δt=0.025 and we first choose the case when p = 0 in the data μ5, μ6 and μ7, in Eq. (28). The analytical Eq. (41) and numerical solutions for l(t), h(t) and a(t) is illustrated in Fig. 5, where the reconstructed timewise terms are in excellent agreement with their corresponding analytical solutions, obtaining with rmse(l) = 6.3E −3, rmse(h) = 1.7E −3 and rmse(a) = 9.8E −3, see Tab. 4.
The analytical Eq. (41) and numerical solutions for: (a) l(t), (b) h(t) and (c) a(t), with no noise and no regularization, for Example 2
The values of rmseEqs. (31) to (33), the function Eq. (27) at final iteration and the number of iterations, for p∈{0,0.5%,1%,2%} noise, for Example 2
p (%)
rmse(l)
rmse(h)
rmse(a)
Minimum value of (27)
Iter
0
0.0063
0.0017
0.0098
3.5E −29
15
0.5
0.0146
0.0083
0.0276
1.8E −28
8
1
0.0385
0.0206
0.0704
9.4E −29
10
2
0.0625
0.0423
0.1180
5.8E −29
10
Next, we add p∈{0.5%,1%,2%} noise into the measured data Eqs. (12) to (14), in order to investigate the stability of the numerical results. The approximate results for the free boundaries l(t), h(t) and the thermal conductivity a(t) are illustrated in Fig. 6. From this figure it can be observed that reasonable and stable numerical results are obtained. The numerical solutions for l(t), h(t) and a(t) obtained with the p=5% noise have been found highly oscillatory and unstable when no regularization was employed, i.e., λ=0, obtaining rmse(l) = 0.1317, rmse(h) = 0.2439 and rmse(a) = 0.2779, see Fig. 7. Therefore, regularization is required in order to restore the stability of the solution in the components l(t), h(t) and a(t). Including regularization, i.e., λ∈{10-3,10-2}, in Eq. (40) alleviates this instability, as shown in the regularized approximate results in Fig. 7 and Tab. 5. The absolute errors between the analytical Eq. (43) and approximate v(y1,y2,t) with the p=5% noise, and without and with regularization parameters, is illustrated in Fig. 8. From this figure and Tab. 5 it can be noticed that the temperature v(y1,y2,t) component is stable and accurate. Overall, the same conclusions can be carried out about the stable determination for the unknown time-dependent coefficients.
The analytical Eq. (41) and numerical solutions for: (a) l(t), (b) h(t) and (c) a(t), with p∈{0.5%,1%,2%} noise, for Example 2
The analytical Eq. (41) and numerical solutions for: (a) l(t), (b) h(t) and (c) a(t), for p=5% noise and without and with regularization, for Example 2
The values of rmseEqs. (31) to (33) and the number of iterations, for p=5% noise, without and with regularization, for Example 2
p (%)
λ
rmse(l)
rmse(h)
rmse(a)
Iter
5
0
0.1317
0.2439
0.2779
18
5
10−5
0.0979
0.0730
0.2217
10
5
10−4
0.0628
0.0460
0.1815
10
5
10−3
0.0293
0.0208
0.1514
10
5
10−2
0.0712
0.0623
0.1776
10
5
10−1
0.2728
0.2964
0.5437
10
The absolute errors between the analytical Eq. (43) and numerical v(y1,y2,t) with (a) λ=0, (b) λ=10-3 and (c) λ=10-2, for p=5%, for Example 2
Conclusions
The inverse problem concerning the reconstruction of the time-dependent thermal conductivity a(t), and free boundary coefficients l(t) and h(t) along with the temperature u(x1,x2,t) in a two-dimensional parabolic equation from the over-specification conditions has been solved for the first time numerically. The forward solver based on the ADE was employed. The inverse problem approach based on a nonlinear least-squares minimization problem using the MATLAB optimization subroutine was developed. The Tikhonov regularization has been employed in order to obtain stable and accurate solutions since the inverse problem is ill-posed. The numerical results for the inverse problem shows that stable and accurate approximate results have been obtained. Extension to three-dimensions is in principle straightforward.
The comments and suggestions made by the referees are gratefully acknowledged.
Funding Statement: The author received no specific funding for this study.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
ReferencesM. J.Huntul and D.Lesnic, “Determination of time-dependent coefficients and multiple free boundaries,” , vol. 5, pp. 15–43, 2017.M. J.Huntul and D.Lesnic, “Time-dependent reaction coefficient identification problems with a free boundary,” , vol. 20, pp. 99–114, 2019.M. J.Huntul and D.Lesnic, “Determination of a time-dependent free boundary in a two-dimensional parabolic problem,” , vol. 5, no. 4, pp. 118, 2019.M. S.Hussein, D.Lesnic, M. I.Ivanchov and H. A.Snitko, “Multiple time-dependent coefficient identification thermal problems with a free boundary,” , vol. 99, pp. 24–50, 2016.B. T.Johansson, D.Lesnic and T.Reeve, “A method of fundamental solutions for the one-dimensional inverse stefan problem,” , vol. 35, pp. 4367–4378, 2011.P.Broadbridge, P.Tritscher and A.Avagliano, “Free boundary problems with nonlinear diffusion,” , vol. 18, pp. 15–34, 1993.G. Q.Chen and M.Feldman, “Free boundary problems in shock reflection/diffraction and related transonic flow problems,” , vol. 373, pp. 20140276, 2015.A.Shidfar and G. R.Karamali, “Numerical solution of inverse heat conduction problem with nonstationary measurements,” , vol. 168, no. 1, pp. 540–548, 2005.Y. C.Hon and M.Li, “A computational method for inverse free boundary determination problem,” , vol. 73, no. 9, pp. 1291–1309, 2008.M.Huntul and M.Tamisr, “Simultaneous identification of timewise terms and free boundaries for the heat equation,” , vol. 38, no. 1, pp. 442–462, 2020.I. G.Malyshev, “Inverse problems for the heat-conduction equation in a domain with a moving boundary,” , vol. 27, pp. 568–572, 1976.J. A.Carrillo and J. L.Vázquez, “Some free boundary problems involving non-local diffusion and aggregation,” , vol. 373, pp. 20140275, 2015.H. A.Snitko, “Inverse problem for a parabolic equation with unknown minor coefficient in a free boundary domain,” , vol. 77, pp. 218–230, 2012.H. A.Snitko, “Determination of the minor coefficients in a parabolic equation in a free boundary domain,” , vol. 81, pp. 142–158, 2016.M. J.Huntul, “Recovering the timewise reaction coefficient for a two-dimensional free boundary problem,” , vol. 7, pp. 66–85, 2019.M. I.Ivanchov, “A problem with free boundary for a two-dimensional parabolic equation,” , vol. 183, pp. 17–28, 2012.H. A.Snitko, “Coefficient inverse problem for a parabolic equation in a domain with free boundary,” , vol. 167, pp. 30–46, 2010.H. A.Snitko, “Inverse problem of finding time-dependent functions in the minor coefficient of a parabolic equation in the domain with free boundary,” , vol. 203, pp. 40–54, 2014.H. A.Snitko, “Inverse coefficient problem for a two-dimensional parabolic equation in a domain with free boundary,” , vol. 68, pp. 1108–1120, 2016.I. E.Barans’ka and I. M.Ivanchov, “Inverse problem for a two-dimensional heat-conduction equation in a domain with free boundary,” , vol. 4, no. 4, pp. 457–484, 2007.J. R.Cannon and J.van der Hoek, “The one phase Stefan problem subject to the specification of energy,” , vol. 86, no. 1, pp. 281–291, 1982.J. R.Cannon and J.van der Hoek, “Diffusion subject to the specification of mass,” , vol. 115, no. 2, pp. 517–529, 1986.D.Mugnolo and S.Nicaise, “The heat equation under conditions on the moments in higher dimensions,” , vol. 288, no. 2–3, pp. 295–308, 2015.Z.Barakat, M.Ehrhardt and M.Gunther, “Alternating direction explicit methods for convection diffusion equations,” , vol. 84, pp. 309–325, 2015.Mathworks, “Documentation Optimization Toolbox-Least Squares (Model Fitting) Algorithms, 2016. [Online]. Available: www.mathworks.com.