Determination of Time-Dependent Coef ﬁ cients for a Weakly Degenerate Heat Equation

: In this paper, we consider solving numerically for the ﬁ rst time inverse problems of determining the time-dependent thermal diffusivity coef ﬁ cient for a weakly degenerate heat equation, which vanishes at the initial moment of time, and/or the convection coef ﬁ cient along with the temperature for a one-dimensional parabolic equation, from some additional information about the process (the so-called over-determination conditions). Although uniquely solvable these inverse problems are still ill-posed since small changes in the input data can result in enormous changes in the output solution. The ﬁ nite difference method with the Crank-Nicolson scheme combined with the nonlinear Tikhonov regularization are employed. The resulting minimization problem is computationally solved using the MATLAB toolbox routine lsqnonlin . For both exact and noisy input data, accurate and stable numerical results are obtained.

Nomenclature a(t)¼thermal diffusivity coefficient ða j Þ j¼1;N ¼discretisation of the function a(t) b(t)¼convection coefficient ðb j Þ j¼1;N ¼discretisation of the function b(t) f¼heat source l¼length of slab p¼percentage of noise t¼time coordinate u¼temperature x¼space coordinate t j , x i ¼finite-difference grid nodes C 0 , C 1 , C 2 ¼positive constants G i,j , G i,0 ¼finite-difference quantities defined in (11) and (12), respectively F 1 , F 2 , F 3 ¼objective functionals defined in (23)-(25), respectively G 1 ¼Green's function defined in (15) H¼Heaviside function M, N¼number of grids in the x−and t−directions, respectively Q T ¼(0.l) × (0,T) T¼final time of interest A a ¼admissible set for a defined in (2) E¼C 2;1 ðQ T Þ \ C 1;0 ðQ T Þ α¼power of degeneracy β¼regularization parameter e1¼ðe1 j Þ j¼1;N ; e2¼ðe2 j Þ j¼1;N ¼random variables defined in (26) and (28), respectively φ¼initial temperature μ 1 , μ 2 ¼boundary temperatures μ 3 ¼heat flux at x=0 μ 4 ¼mass/energy of the thermal system σ 1 , σ 2 ¼standard deviations defined in (27) Δ x, Δt¼finite-difference steps Ψ¼function defined in (13) 1 Introduction Inverse problems concerning simultaneously determining several time-dependent coefficients for non-degenerate partial differential equations with fixed or even moving boundaries have been investigated in several works, see ; Huntul, Lesnic and Hussein (2017) ;Hussein, Lesnic and Ivanchov (2014) ;Hussein, Lesnic, Ivanchov et al. (2016); Ivanchov andPabyrivska (2001, 2002)], to mention only a few. For instance, [Huntul, Lesnic and Hussein (2017)], have investigated the inverse problems of simultaneous numerical reconstruction of time-dependent thermal conductivity/convection/absorption coefficient from heat moments. Determination of time-dependent coefficients and multiple free boundaries has been investigated ]. However, only a few papers are concerned with weakly or strongly degenerate parabolic equations [Hryntsiv (2009[Hryntsiv ( , 2011; Huzyk (2014Huzyk ( , 2015Huzyk ( , 2016; Ivanchov and Saldina (2005); Saldina (2005); Vlasov (2014)]. These studies are theoretical and they are important because they establish sufficient conditions for the unique solvability of the time-dependent coefficient identification problems. However, no numerical reconstruction has been attempted and it is the purpose of the present study to numerically recover the unknown coefficients in a stable and accurate manner. Therefore, in this paper inverse problems concerned with determining the time-dependent thermal diffusivity and convection coefficients for a weakly degenerate parabolic heat equation, together with the temperature from over-determination data is, for the first time, numerically solved. It is supposed that the thermal diffusivity coefficient vanishes at the initial moment of time. Here we investigate the case of weakly degeneration where the degeneracy is given by a time-dependent power law t α with α 2 (0,1). The case of strong power law t α degeneration with α≥1 can also be investigated [Ivanchov and Saldina (2006)].
The structure of the paper is as follows. The mathematical formulations of the inverse problems are described in Section 2. In Section 3, the numerical solution of the direct problem is based on finite difference method with the Crank-Nicolson scheme. The treatment for solving the degenerate parabolic equation is discussed. The numerical solutions of the inverse problems are obtained using the Tikhonov regularization method, as described in Section 4. The numerical results for a few weak degenerate inverse problems are presented and discussed in Section 5. Finally, the conclusions are presented in Section 6.

Mathematical formulations of the inverse problems
We consider the convection-diffusion equation in a finite slab of length l>0 over a time duration T>0 satisfying the parabolic partial differential equation (PDE) @u @t ðx; tÞ ¼ aðtÞ @ 2 u @x 2 ðx; tÞ þ bðtÞ @u @x ðx; tÞ þ f ðx; tÞ; ðx; tÞ 2 Q T :¼ ð0; lÞ Â ð0; T Þ; where u is the unknown temperature, f is a given heat source, and a and b are the time-dependent thermal diffusivity and convection coefficients, respectively, which may be known or unknown. For simplicity, we have assumed that no reaction term c(x,t)u(x,t) is present in (1). In nondegenerate uniformly parabolic equations, the diffusivity a(t)≥a 0 >0 for all t 2 [0,T], for some positive known constant a 0 . However, in certain porous media practical applications related to clogging or sink-hole [Arbogast and Taicher (2016)], the diffusivity coefficient may vanish, making the PDE (1) degenerate and non-uniformly parabolic. In this paper, we consider such as a degeneracy and assume that the coefficient a belongs to the class A a :¼ a 2 C½0; T j aðtÞ > 0; t 2 ð0; T and there exists the finite lim t&0 aðtÞ t a > 0 ; where α 2 (0,1) is the given degree of weakly power law degeneration. The case α≥1 corresponding to strong degeneration will be investigated in a separate work. Eq. (1) is subjected to the initial condition uðx; 0Þ ¼ 'ðxÞ; x 2 ½0; l; (3) and the Dirichlet boundary conditions uð0; tÞ ¼ l 1 ðtÞ; uðl; tÞ ¼ l 2 ðtÞ; t 2 ½0; T: We next formulate three inverse problems with respect to whether the coefficients a(t) and/or b(t) are known or unknown and state sufficient conditions for the uniqueness of solution.
Condition (B) requires that the heat flux (5) does not vanish at any instant t 2 (0,T], but it behaves like the power law degeneracy t α for small t. A two-dimensional variant of the IP1 has been considered theoretically [Vlasov (2014)], but its numerical simulation in the context of our investigation is deferred to a future work.
Theorem 2.2 Assume that the following conditions are satisfied: Then, the IP2 given by Eqs.
Theorem 2.3 Let the assumptions (B) and (D) hold and assume also that the following condition is satisfied: Then, the IP3 given by Eqs.
Away from t=0 the heat Eq. (1) is non-degenerate and it can be approximated as usual using the Crank-Nicolson FDM given by (10) and (11), to read as follows: 4 Numerical solution for the inverse problems We wish to obtain stable reconstructions of the unknown coefficients a(t) and/or b(t) together with the temperature u(x,t), by minimizing the nonlinear Tikhonov regularization function or, or, respectively, where u solves (1), (3) and (4) for given a and b, and β≥0 is regularization parameter to be prescribed. The minimization of F 1 , or F 2 , or F 3 is performed using the MATLAB toolbox routine lsqnonlin, which does not require the user to supply the gradient of the objective function. This routine attempts to find the minimum of a sum of squares by starting from a given initial guess. Furthermore, within lsqnonlin, we use the Trust Region Reflective (TRR) algorithm [Coleman and Li (1996)], which is based on the interior-reflective Newton method. In the numerical computation, we take the parameters of the routine as follows: • Maximum number of iterations, (MaxIter)=400.
The IP1, IP2 and IP3 are solved subject to both exact and noisy measurements (5) and (6). The noisy data are numerically simulated as follows: where ε1 j and ε2 j are random variables generated from a Gaussian normal distribution with mean zero and standard deviations σ1 and σ2 given by where p represents the percentage of noise. We use the MATLAB function normrnd to generate the random variables e1 ¼ ðe1 j Þ j¼1;N and e2 ¼ ðe2 j Þ j¼1;N as follows: 5 Numerical results and discussion In this section, we present examples for IP1, IP2 and IP3 in order to test the accuracy and stability of the numerical methods introduced in Section 3 based on the FDM combined with the minimization of the objective function F 1 , or F 2 , or F 3 , as described in Section 4. Furthermore, we add noise to the input data (5) and (6) to simulate the real situation of noisy measurements, by using Eqs. (26)-(28). To assess the accuracy of the approximate solutions, we introduce the root mean square l 2 − errors (rmse) defined as follows: For simplicity, we take l=T=1 in all examples. We take the lower and upper simple bounds for a(t) to be 0 and 10 2 , and for b(t) to be −10 2 and 10 2 , respectively. These bounds allow a wide search range for the unknowns. In the FDM, we take M=N=40. We also take α=0.5 as a typical degree of weak power law degeneracy in (2). In what follows, the expressions (34), (1), (3)-(5) with unknown thermal diffusivity a(t), and input data bðtÞ ¼ 0; One can observe that the assumptions (A) and (B) of Theorem 2.1 are satisfied and thus the solution of IP1 is unique, if it exists. It can easily be checked by direct substitution that the analytical solution for the temperature u(x,t) is and for the thermal diffusivity a(t) is We take the initial guess for a(t) as a 0 (t)=t/6 for t 2 (0,1], knowing that since a 2 A a we must have a(0)=0.
Before we attempt any finite-difference numerics it is important to calculate the function Ψ(x) given by Eq. (13) since its value is needed in initiating the FDM time-marching procedure in Eq. (12). With the data (31) and (32), Eqs. (2), (13) and (16) It is sufficient to estimate the term corresponding to n=0 in the above, because the integration function remaining after eliminating the summand corresponding to n=0 has no singularities [Ivanchov (2003)], and thus, when multiplied with a(t), it has the zero limit as t a 0. Considering only the term corresponding to n=0 in (37) we obtain (after dropping the constant 1=2 ffiffiffi p p and using that e τ ≤ e t ≤ e T ≤ e) that for some positive constants C 1 and C 2 , where we have used the definition of hðtÞ ¼ R t 0 aðsÞds and that a 2 A 1=2 . Using the change of variable z ¼ x= in the integral on the right-hand side of (38) we obtain that So, from (36)-(38), it follows that for example 1.1, we have Ψ(x) ≡ 0 in (12), (13) and (36).
We attempt to recover the unknown thermal diffusivity a(t) and the temperature u(x,t) for exact input data, i.e., p=0 in (25), as well as for p 2 {1,2,3}% noisy data. The unregularized objective function F 1 given by (23) with β=0, as a function of the number of iterations, is plotted in Fig. 1. From this figure, it can be seen that the objective function F 1 is rapidly decreasing to a very low value of O(10 −25 ) in about 10 iterations (in less than 3 minutes CPU time). The related numerical results for the thermal diffusivity a(t) are presented in Fig. 1. From this figure it can be seen that there is good agreement between the numerical results and the analytical solution (35) for exact data, i.e., p=0, and consistent with the errors in the input data for p>0. The numerical solution for the thermal diffusivity converges to the analytical solution (35), as the percentage of noise p decreases, with rmse(a) 2 {0.0001, 0.0044, 0.0088, 0.0132} for p 2 {0,1,2,3}%, respectively.
Finally, Fig. 2 shows the absolute error between the exact solution (34) and the numerical solutions for the temperature u(x,t) for various amounts of noise p 2 {0,1,2,3}%. From this figure it can be seen that the numerical solution is stable and furthermore, its accuracy improves as the noise level p decreases.
Example 1.2 Consider now a more complicated test example for the temperature u(x,t) given by the oscillatory analytical solution uðx; tÞ ¼ In addition of being a more complicated test example than example 1.1, condition (B) in Theorem 2.1 is only partially fulfilled since the above flux μ 3 (t) vanishes at a point in the interval (0,1]. However, similarly as computed for example 1.1, we obtain that Ψ (x)≡0. Analogous to Fig. 1, Fig. 3 presents the unregularized objective function F 1 with β=0 and the numerical results for the thermal diffusivity for various percentages of noise p 2 {0,1,2,3}%. It can be seen that a rapid decrease in the objective function F 1 is realised in about 10 iterations. Moreover, on comparing Figs. 1 and 3, it can be observed that the numerical retrieved solutions for the more complicated example 1.2 is less stable than that for example 1.1. However, this instability can be mitigated by including regularization with β>0, as illustrated in the next subsections for the problems IP2 and IP3.

Example 2 (for IP2)-Finding b(t) when a(t) is known
Consider the IP2 given by Eqs. (1), (3), (4) and (6) with unknown convection coefficient b(t), and input data (31), (35) and One can remark that the assumptions (C) and (D) of Theorem 2.2 are satisfied and thus the solution of IP2 is unique, if it exists. In fact, the analytical solution for the unknown convection coefficient b(t) is given by (32) and for the temperature u(x,t) by (34). Although the analytical solution (32) for b(t) is trivial, the numerical rmse(b) given by (30) can still be calculated and is meaningful.
We consider first the case where there is no noise (i.e., p=0) included in the input data μ 4 (t) in (39). The objective function F 2 , as a function of the number of iterations is displayed in Fig. 4. From this figure, it can be seen that the decreasing convergence of the objective function F 2 is very fast and is achieved in 5 iterations (in 2 minutes CPU time) to reach a stationary value of O(10 −29 ). The corresponding numerical results of the time-dependent convection coefficient b(t) are depicted in Fig. 4 and accurate results of O(10 −4 ) error can be observed.
Next, we investigate the stability of the IP2 with respect to noise. We include p=1% noise to the data (39) simulated numerically, via Eq. (26) for μ 4 (t). The rmse values (30) versus the number for the unknown convection coefficient b(t) are presented in Fig. 5 with and without regularization, versus the number of iterations. It can be seen that the rmse values settle rapidly to a stationary value O(10 −3 ) after 2 to 3 iterations when regularization is included, but in case of no regularization they increase with the number of iterations, as expected since the unregularized solution is unstable. In Fig. 6 and Tab. 1 we present the unknown convection coefficient b(t) and the rmes(b) given by Eq. (30), the number of iterations and computational time. It can be seen that the numerical results for the i.e., β=0, is employed, or even if β is too small such as 10 −3 . At the other extreme, if β is too large this would penalise too much the solution norm resulting in an increased residual in the functional (24) that is minimized. Clearly, one can observe the effect of the regularization parameter β>0 in decreasing the oscillatory unstable behaviour of the convection coefficient b(t). Overall, the numerical results obtained with β 2 {1,2} seem stable and accurate, see Fig. 6 and Tab. 1. A further discussion on how the regularization parameter can be chosen is explained later on in Section 6.

Example 3 (for IP3)-Finding a(t) and b(t) together
Consider the IP3 given by Eqs.
We start first with the case of exact data, i.e., p=0. Fig. 7 illustrates the exact and numerical coefficients a(t) and b(t) plotted after 6 iterations of minimization of the objective function F 3 in (25) without regularization, i.e., β=0. From this figure, it can be seen that a very good  We next include p=1% in the input data (33) and (39). From the previous IP2 analysis, we anticipate that regularization is needed in order to achieve stable and accurate results. The numerical results for the thermal diffusivity a(t) and the convection coefficient b(t) for p=1% noise are presented in Fig. 8. From this figure it can be seen that stable and reasonable accurate numerical results are obtained for β 2 {10 −2 ,10 −1 }.

Conclusions
This paper has presented the determination of time-dependent thermal diffusivity coefficient and/or the convection coefficient for a weakly degenerate heat equation from heat flux and/or mass/energy measurement/specification/over-determination. Three coefficient identification problems (termed IP1, IP2 and IP3) have been investigated. The uniqueness of solution holds under easy verifiable sufficient conditions on the input data, as proven in the previous theoretical literature [Hryntsiv (2009);Huzyk (2014); Ivanchov and Saldina (2005); Saldina (2005)], but without numerical realisation.
The resulting inverse problems have been reformulated as constrained regularized minimization problems which were solved using the MATLAB optimization toolbox routine lsqnonlin. The nonlinear Tikhonov regularization has been employed in order to obtain stable and accurate results because the inverse problems under investigation are ill-posed and sensitive to noise. The uncertainty analysis has been considered through the random generation of noise in Eq. (26) in order to test the stability of the proposed method and to stress the need of using regularization. The numerically obtained results  (35) and (32)) and numerical solutions for the thermal diffusivity a(t) and the convection coefficient b(t), with p=1% noise for example 3, with and without regularization are stable and accurate. A complete statistical analysis which may include confidence/ credibility intervals for the retrieved coefficients is beyond the scope of the present paper and is deferred to a future work.
The main difficulty in regularization when we solve the IP2 or IP3 is how to choose an appropriate regularization parameter β which compromises between accuracy and stability. However, one can use techniques, such as the L-curve method or Morozov's discrepancy principle, to find such a parameter, but in our work we have used trial and error. As mentioned in Dennis et al. [Dennis, Dulikravich and Yoshimura (2004)], the regularization parameter β is selected based on experience by first choosing a small value and gradually increasing it until any numerical oscillations in the unknown coefficients disappear. It would also be interesting the use of the time step Δt as a self-regularization parameter [Joachimiak, Joachimiak, Cialkowski et al. (2019)], but this investigation is deferred to a future work.
As this is the first investigation into numerically solving the IP1, IP2 and IP3, there are no other studies/methods to compare our results with. However, we stress that when solving parabolic PDEs which are degenerate at the initial time t=0, the numerical challenge is how to calculate the function Ψ(x) in (13) in order to initiate, via (12), the time-marching Crank-Nicolson FDM process. Another way could be to employ the backward FDM, but this would ignore the degeneracy at t=0 altogether. Probably better and deferred to a possible future work would be to employ the method of Green's functions or the boundary element method, after the transformation hðtÞ ¼ R t 0 aðsÞds [Ivanchov and Saldina (2005)].