=Paper=
{{Paper
|id=Vol-2642/paper2
|storemode=property
|title=On Fast Automatic Differentiation To Solving The Inverse Coefficient Problems
|pdfUrl=https://ceur-ws.org/Vol-2642/paper2.pdf
|volume=Vol-2642
|authors=Yuri G. Evtushenko,,Vladimir I. Zubov,,Alla F. Albu
}}
==On Fast Automatic Differentiation To Solving The Inverse Coefficient Problems==
On Fast Automatic Differentiation To Solving The Inverse Coefficient Problems Yuri G. Evtushenko Vladimir I. Zubov Alla F. Albu Federal Research Center Federal Research Center Federal Research Center ”Computer Science and ”Computer Science and ”Computer Science and Control” of RAS Control” of RAS Control” of RAS Moscow, Russia 119333 Moscow, Russia 119333 Moscow, Russia 119333 yuri-evtushenko@yandex.ru vladimir.zubov@mail.ru alla.albu@yandex.ru Abstract The inverse problem under consideration is to determine a temperature- dependent thermal conductivity coefficient from experimental observa- tions of the temperature field in the studied substance and (or) the heat flux on the surface of the object. The study is based on the Dirichlet boundary value problem for the nonstationary heat equation stated in the general n-dimensional formulation. For the numerical solution of the problem an algorithm based on the modern Fast Automatic Differ- entiation technique is proposed. 1 Introduction Recently, the need to create new materials has increased and this has led to a growing interest in inverse problems. This is due to the fact that some characteristics of these materials are unknown in advance, and they can be determined by solving inverse problems. In the description and mathematical modeling of many thermal processes, the classical heat equation is often used. The density of the substance, its specific thermal capacity, and the thermal conductivity coefficient appearing in this equation are assumed to be known functions of the coordinates and temperature. Additional boundary value conditions makes it possible to determine the dynamics of the temperature field in the substance under examination. However, often happens that the thermal conductivity coefficient depends only on the temperature, and this dependence is not known. In this case, the problem of determining the dependence of the thermal conductivity coefficient on the temperature based on experimental measurements of the temperature field arises. This problem also arises when a complex thermal process should be described by a simplified mathematical model. For example, in studying and modeling heat propagation in complex porous composite materials, where the radiation heat transfer plays a considerable role, both the convective and radiative heat transfer must be taken into account. The thermal conductivity coefficients in this case typically depend on the temperature. To estimate these coefficients, various models of the medium are used. As a result, one has to deal with a complex non-linear model that describes the heat propagation in the composite material (see [Dul’nev1974], [Alifanov2009], [Zverev2008]). However, another approach is possible: a simplified model is constructed in which the radiative heat transfer is Copyright ⃝ c by the paper’s authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). In: Sergei S. Goncharov, Yuri G. Evtushenko (eds.): Proceedings of the Workshop on Applied Mathematics and Fundamental Computer Science 2020 , Omsk, Russia, 30-04-2020, published at http://ceur-ws.org 1 not taken into account, but its effect is modeled by an effective thermal conductivity coefficient that is determined based on experimental data. Determining the thermal conductivity coefficient of a substance is an important problem, and it has been stud- ied for a long time. This is confirmed by a large number of publications (e.g., see [Kozdoba1982], [Alifanov1988], [Vabishchevich1990], [Samarskii1997], [Samarskii2003], [Marchuk1992]). In these works, the inverse coefficient problems are studied theoretically, and numerical methods for their solution are developed. In [Czel2012] simultaneous inverse identification of the temperature-dependent volumetric heat capacity and thermal conductivity of a solid material based on transient temperature histories was studied. There, it is pro- posed to use a genetic method for solving inverse problems, which, as we know, refers to zero-order optimization methods. In other sources, inverse problems are considered only for the case when the desired coefficient of thermal conductivity does not depend on temperature. The inverse coefficient problem considered here is nonlinear. The study is based on the Dirichlet boundary value problem for the nonstationary heat equation stated in the general n-dimensional formulation. 2 Statement Of The Problem Consider a restricted domain Q ⊂ Rn with piecewise-smooth boundary S = ∂Q. This domain is filled with the substance being investigated. The distribution of the temperature field in Q at each moment is described by the following initial boundary value (mixed) problem: ∂T (x, t) C(x) = divx (K(T (x, t))∇x T (x, t)), x ∈ Q, 0 < t ≤ Θ, (1) ∂t T (x, 0) = w0 (x), x ∈ Q, (2) T (x, t) = ws (x, t), x ∈ S, 0 ≤ t ≤ Θ. (3) Here x = (x1 , ..., xn ) are the Cartesian coordinates; t is time; T (x, t) is the temperature of the material at the point with the coordinates x at time t; C(x) is the volumetric heat capacity of the material; K(T ) is the thermal conductivity coefficient; w0 (x) is the given temperature at the initial time t = 0; ws (x, t) is the given temperature on the boundary of the object. The volumetric heat capacity of a substance C(x) is considered a known function of the coordinates. If the dependence of the thermal conductivity coefficient K(T ) on the temperature T is known, then we can solve the mixed problem (1)-(3) to find the temperature distribution T (x, t) in Q × (0, Θ]. Problem (1)-(3) below we will call primal problem. The inverse coefficient problem is reduced to the following variational problem: find the dependence K(T ) on T under which the temperature field T (x, t), obtained( by solving the mixed ) problem (1)-(3), is close to the field Y (x, t) obtained experimentally, and the heat flux −K(T (x, t)) ∂T∂n (x,t) on the boundary of the domain is close to the experimentally flux P (x, t). The quantity ∫Θ∫ 2 Φ(K(T )) = [T (x, t) − Y (x, t)] · µ(x, t)dx dt + 0 Q ∫Θ∫ [ ]2 ∫b ∂T (x, t) + β(x, t) −K(T (x, t)) − P (x, t) ds dt + ε (K ′ (T ))2 dT (4) ∂n 0 S a can be used as the measure of difference between these functions. Here, ε ≥ 0, β(x, t) ≥ 0, µ(x, t) ≥ 0 are given weight parameters; Y (x, t) is a given temperature field; P (x, t) is a given heat flux at the boundary S of the domain Q, ∂T∂n is the derivative of the temperature along the direction of the outer normal to the boundary of the domain; [a; b] is the interval on which the function K(T ) will be restored. The last term in functional (4) is used to obtain a smooth solution of variational problem with moderate perturbations of experimental data. Thus, the inverse problem is reduced to the following optimal control problem: find the optimal control K(T ) and the corresponding solution T (x, t) of the problem (1)-(3) that minimizes functional (4). Optimal control problems of this type are usually solved numerically with the help of some gradient descent, which requires knowledge of the gradient of functional (4). Moreover, it is of utmost importance to deter- mine the exact value of this gradient. In [Zubov2016] an analytical expression of the functional gradient was 2 obtained for a one-dimensional spatial problem, and in [Albu2018] an analytical gradient was obtained for a two-dimensional problem. When determining the functional gradient, the classical apparatus was used in these works ([Sirazetdinov1997], [Kraiko1979], [Vasil’ev2002]). However, it is almost impossible to use an analytical expression for the functional gradient in the numerical solution of the problem due to its complexity. Accordingly, the gradient of the functional in the proposed algorithm was computed by applying the efficient FAD-technique. The time cost of this technique was at most three times the CPU time required for computing the function. The efficiency of the FAD-technique as applied to computing the function gradient was ensured by using the solution of the adjoint problem. The construction of an approximation of the adjoint problem that is consistent with the chosen approximation of the primal problem is the basic advantage of the FAD. The FAD-technique delivers canonical formulas calculated by which the gradient of the cost function is exact for the selected approximation of the optimal control problem. 3 Algorithm For Solving A Primal Problem A major element of the algorithm proposed for solving the inverse coefficient problem is the solution of the primal problem. For the numerical solution of the primal problem, we introduced grids (generally nonuniform) in time and space. At each grid node of the computational domain Q × [0, Θ], all functions were specified by their values at this point. The heat equation was approximated by a suitable finite-difference scheme. In [Zubov2016], the problem of identifying the thermal conductivity coefficient in the one-dimensional case is investigated. There is considered a layer of material of width L is considered. The inverse problem was reduced to the following variational problem: find such dependence K(T ) on temperature at which the temperature field and the heat flux on the left boundary of the sample (left hand side), obtained by solving the direct problem, differ little from the data obtained experimentally. To approximate the direct problem in the one-dimensional case, the domain [0, L] × [0, Θ] is decomposed by the grid lines {xei }Ii=0 and {e tj }Jj=0 into rectangles. A two-layer implicit scheme was used to discretize the thermal conductivity equation: ( ) ( ) Tij = bji K(Tij ) + K(Ti+1 j j ) (Ti+1 − Tij ) − aji K(Tij ) + K(Ti−1 j ) (Tij − Ti−1 j )+ ( ) +Tij−1 + cji K(Tij−1 ) + K(Ti+1 j−1 j−1 ) (Ti+1 − Tij−1 )− ( ) −dji K(Tij−1 ) + K(Ti−1j−1 ) (Tij−1 − Ti−1 j−1 ), (5) aji = στ j /(Ci hi−1 (hi + hi−1 )), bji = στ j /(Ci hi (hi + hi−1 )), cji = (1 − σ)τ j /(Ci hi (hi + hi−1 )), dji = (1 − σ)τ j /(Ci hi−1 (hi + hi−1 )), ei , τ j = e ei+1 − x σ – weight parameter, hi = x tj − e tj−1 , i = 1, I − 1, j = 1, J. Ti0 = (w0 )i , i = 0, I, T0j = w1j , TIj = w2j , j = 0, J. The resulting system of nonlinear algebraic equations was solved iteratively using the Gaussian elimination. This approach was used to solve the mixed problem (1)-(3) and the function (more precisely, its approximation Tij ) was found. When studying the inverse problem in the two-dimensional case, we consider a rectangular plate of length L and width R. For the numerical solution of the problem, we introduced grids in time and space (generally nonuniform). Two spatial grids were also introduced: basic and auxiliary. To construct the basic spatial grid on the segment N [0, L], a system of nodal points {xn }n=0 was chosen so that x0 = 0, xN = L and xn < xn+1 for all 0 ≤ n < N . Distance between xn and xn+1 was denoted as hxn , i.e. hxn = xn+1 − xn , n = 0, N − 1. Similarly, on the segment I [0, R] a system of nodal points {yi }i=0 was chosen so that y0 = 0, yI = R and yi < yi+1 for all 0 ≤ i < I. In this y case hi = yi+1 − yi , i = 0, I − 1.The points of the basic grid are the set of points {xn , yi } where n = 0, N and i = 0, I. Auxiliary grid {e xn , yei }, n = 0, N + 1 and i = 0, I + 1 was built in a similar way: e0 = x0 , x eN +1 = xN , x en = xn−1 + hxn−1 /2, x n = 1, N , 3 ye0 = y0 , yeI+1 = yI , yei = yi−1 + hyi−1 /2, i = 1, I. The lines x = x en , y = yei divide the domain Q into computational cells. An elementary cell is assigned the xn , yei ), and its area is indices (n, i) if the cell vertex nearest to the point (x0 , y0 ) coincides with the grid point (e denoted by Sni , i.e. Sni = (e xn+1 − x yi+1 − yei ), n = 1, N − 1, i = 1, I − 1. en )(e j We also denote the average temperature in the cell with indices (n, i) through Tni (t), and through Tni - j its value at time t . The heat equation was approximated using the implicit alternating direction scheme (see [Peaceman1955]) with splitting in x and y, so two subproblems were formed: x − direction j+ 1 j j+ 1 j+ 1 j+ 1 j+ 1 qni (Tni 2 −Tni )/τ j = (K(Tn+1,i 2 ) + K(Tni 2 ))(Tn+1,i 2 −Tni 2 )(hyi + hyi−1 )/(8hxn )− j+ 1 j+ 1 j+ 1 j+ 1 −(K(Tni 2 ) + K(Tn−1,i 2 ))(Tni 2 −Tn−1,i 2 )(hyi + hyi−1 )/(8hxn−1 )+ j j j j +(K(Tn,i+1 ) + K(Tni ))(Tn,i+1 −Tni )(hxn + hxn−1 )/(8hyi )− j j j j −(K(Tni ) + K(Tn,i−1 ))(Tni −Tn,i−1 )(hxn + hxn−1 )/(8hyi−1 ), (6) y − direction j+1 j+ 21 j+1 j+1 j+1 j+1 qni (Tni −Tni )/τ j = (K(Tn,i+1 ) + K(Tni ))(Tn,i+1 −Tni )(hxn + hxn−1 )/(8hyi )− j+1 j+1 j+1 j+1 −(K(Tni ) + K(Tn,i−1 ))(Tni −Tn,i−1 )(hxn + hxn−1 )/(8hyi−1 )+ j+ 1 j+ 1 j+ 1 j+ 1 +(K(Tn+1,i 2 ) + K(Tni 2 ))(Tn+1,i 2 −Tni 2 )(hyi + hyi−1 )/(8hxn )− j+ 1 j+ 1 j+ 1 j+ 1 −(K(Tni 2 ) + K(Tn−1,i 2 ))(Tni 2 −Tn−1,i 2 )(hyi + hyi−1 )/(8hxn−1 ), (7) (n = 1, N−1, i = 0, I−1, j = 1, J), qni = Sni · ρni Cni . j j+ 1 In (6)–(7) along with the basic values Tni of the grid function, the intermediate values Tni 2 of the temperature function are introduced, which can be formally considered as the values of the temperature at the time moment tj + τ j /2 (see [Samarskii2001]). The initial and boundary conditions (2)–(3) are written in the difference form: 0 Tni = (w0 )ni , (n = 0, N , i = 0, I), j T0i = (w1 )ji , TNj i = (w2 )ji , (i = 0, I, j = 1, J), j j Tn0 = (w3 )jn , TnI = (w4 )jn , (n = 0, N , j = 1, J). j+ 1 The system of nonlinear algebraic equations (6) is solved with respect to unknown values Tni 2 . The system of equations (6) is divided into a number of subsystems. Each subsystem is characterized by a fixed index i and is solved independently of the other subsystems. These subsystems of nonlinear equations are solved by an iterative method, and at each iteration the tridiagonal Gaussian elimination is used. The system of nonlinear j+1 j+ 1 algebraic equations (7) is solved with respect to unknown values Tni . The calculated values Tni 2 are used as the initial conditions for the second subproblem. When considering the problem in three-dimensional space, the sample of the studied substance has the form of a straight parallelepiped of length L, width R and height H. The basic and auxiliary spatial grids are constructed in the same way as in the case of a two-dimensional problem. In this case, the choice of a difference scheme for solving a three-dimensional non-stationary heat equation becomes more complicated. To discretize the initial boundary value problem (1)–(3), we recommend using one of the following schemes of alternating directions: the Pisman-Reckford scheme, the Douglas-Reckford scheme, or a locally one-dimensional scheme. The Douglas-Reckford scheme ([Douglas1956]) and the locally one-dimensional scheme ([Samarskii2003]) are stable and are most often used by researchers in practice. The Pisman-Reckford scheme ([Peaceman1955], [Thibault1985]) for the three-dimensional heat equation is conditionally stable. In [Gao1996], using the von Neumann criterion, was analyzed the stability of the Pisman- Reckford difference scheme for the n-dimensional heat-conductivity equation and showed that for n ≥ 3 it is conditionally stable even if the coefficients of the equation do not depend on temperature. It requires working with a much smaller time step than the first two schemes. However, sometimes it is the Pisman-Reckford scheme 4 may be the most effective. Therefore, when solving each specific problem, it is necessary to conduct research on the choice of the difference scheme. The algorithm for solving finite-difference equations, which are obtained as a result of discretization of the three-dimensional heat equation using one of the alternating directions schemes indicated here, is similar to the algorithm described above, which is used in the two-dimensional case. 4 Numerical Solution Of The Optimal Control Problem The formulated optimal control problem was solved numerically. The temperature interval [a, b] is partitioned by the points Te0 = a, Te1 , Te2 , . . . , TeM = b into M parts (they can be equal or of different lengths). Each point Tem (m = 0, . . . , M ) is assigned a number km = K(Tem ). The function K(T ) to be found is approximated by a continuous { }M piecewise linear function with the nodes at the points (Tem , km ) so that K(T ) = km−1 + km −km−1 (T −Tem−1 ) Tem −Tem−1 m=0 for Tem−1 ≤ T ≤ Tem , (m = 1, . . . , M ). If the temperature at the point fell outside the boundaries of the segment [a, b], then the linear extrapolation was used to determine the function K(T ). The cost functional was approximated by a function F (k0 , k( 1 , . . . , kM ) of the finite) number of variables using ∂T (x,t) the rectangles method. An approximation of the expression −K(T (x, t)) ∂n by which the heat flux at the object boundary is calculated is presented in [Zubov2016] (one-dimensional case) and [Albu2018] (two- dimensional case). The optimal control problem was solved numerically using gradient methods to minimize the cost function. The effective Fast Automatic Differentiation technique was used to calculate the functional gradient. Formulas for calculating the gradient components and the discrete conjugate problem with which these com- ponents are calculated, in the case of studying a one-dimensional spatial problem, were obtained in [Zubov2016]. In [Albu2018-2], the inverse coefficient problem was considered on the basis of the first boundary-value prob- lem for the two-dimensional unsteady heat equation. The root-mean-square deviation of the temperature field calculated as a result of solving the problem from the experimental data was chosen as the cost functional there (it was assumed that in functional (4) β(x, t) ≡ 0). In [Albu2018], the root-mean-square deviation of the calcu- lated heat flux on the body surface from the experimentally obtained flux was chosen as the cost functional (it was assumed that in functional (4) µ(x, t) ≡ 0). The discrete conjugate problem and the formula for calculating the gradient of the cost function in both cases were obtained based on the approximation of the primal problem using the Peaceman–Rachford scheme. Numerous test calculations have shown the effectiveness of the proposed algorithm. The thermal conductivity coefficient of a substance is reproduced with high accuracy. The examples below illustrate this. The considered examples also revealed some features of the considered inverse problem and the proposed algorithm. 1) If the parameter β(x, t) ≡ 0 in the cost functional, then the inverse coefficient problem may have a non- unique solution. For example, in [Zubov2016]√ we consider a one-dimensional problem (0 ≤ x ≤ 1) for which Θ = 1, C(x) ≡ 1. The function Y (x, t) = 2(t + 1.5 − x) is the solution of a mixed problem (1)–(3) for w0 (x) = Y (x, 0), wS (x, t) = Y (x, t)|x∈S and K(T ) = T 2 . Figure 1 shows two optimal controls opt1 and opt2 obtained from different initial approximations that do not coincide with the control K(T ) = T 2 . Figure 1: Illustration of non-uniqueness of the solution of inverse problem 5 In [Zubov2016], the necessary conditions on the field are presented under which non-uniqueness of the solution of the inverse problem arises. To single out a unique solution of the optimal control problem, we suggest specifying a point T∗ at which the thermal conductivity coefficient is known, i.e., K∗ = K(T∗ ). If the approximate function K(T ) passes through the given point T∗ at each step of the minimization process, then the solution of the inverse problem is unique. Figure 2 illustrates the convergence of the approximate solutions of the inverse problem (the thermal conductivity coefficient) to the exact solution K(T ) = T 2 . The line marked by 0 corresponds to the initial solution. The further approximations to the solution of the identification problem (to the function K(T )) are shown in several iterative steps and are marked by increasing numbers. The point with the coordinates (T∗ , T∗2 ), which is marked by a circle, belongs to all approximate solutions. It is important to note that if we use functional (4) with β(x, t) > 0, then a unique solution is obtained (even if µ(x, t) ≡ 0). Figure 2: Identification of the thermal conductivity coefficient in the case of setting a single point 2) A large number of calculations of the optimal control problem (1)-(4) in the two-dimensional case were performed (see [Albu2018], [Albu2018-2]). For example, the experimental field was determined as a result of solving the primal problem (1)-(3), for which Θ = 1, C(x) ≡ 1, the coefficient of thermal conductivity { 0.1 · (T − 3)(T − 6)(T − 7) + 3.4, T ≥ 3, K∗ (T ) = 1.2 · (T − 3) + 3.4, T < 3, and with such conditions on the parabolic boundary of the domain that are there the trace of the function g(x1 , x2 , t) = 9/(x1 + 2x2 + 5t + 1), (0 ≤ x1 ≤ 1), (0 ≤ x2 ≤ 1). It was assumed that µ(x1 , x2 , t) ≡ 1 and β(x1 , x2 , t) = 0. Analysis of the constructed temperature field allowed us to determine the range of temperature changes: a = 1, b = 9. It was found that the gradient of the cost functional is distributed over temperature strongly nonuniformly, which noticeably deteriorates the convergence of the iterative process. This difficulty in solving the problem can be effectively overcome by applying an approach based on a sequential increase in the number M of partitions of the interval [a, b]. It is desirable to begin the process with M = 1. After an optimal solution has been obtained, it is used as an initial approximation for the variant with M = 2. The optimal solution found for M = 2 is used as an initial approximation for M = 4, etc. Figure 3 shows the initial control (ini) and the optimal controls obtained for M = 4, M = 8 and M = 128(opt). It can be seen that a practically acceptable solution is already obtained for M = 8. As for the solution for M = 128(opt), here the cost functional decreased by 15 orders of magnitude, the maximum gradient modulus by 11 orders of magnitude, and the deviation of the obtained thermal conductivity coefficient K(T ) from its analytical value K∗ (T ) did not exceed 10−7 . It should be noted that in this example the solution of the inverse problem is unique. Calculations performed for different parameters of the optimization process (initial approximation, size of the computational grid, number of temperature interval, parameters of the gradient descent algorithm) always led to the same solution. 6 Figure 3: Recovery of the thermal conductivity coefficient in two- dimensional case The above studies are made on the assumption that all the initial data are accurate. In this case, the parameter ε in the cost functional (4) was assumed to be zero. 3) A large series of calculations was performed to evaluate the effect of the error in the experimental data on the accuracy of the obtained solution of the inverse problem. The perturbed experimental data were constructed by superimposing some random field on the unperturbed data. As shown by a numerical study of the stability of the obtained solutions, the perturbation of the desired thermal conductivity coefficient has the same order as the perturbations that caused it in the experimental data. 5 Conclusion An algorithm for determining the thermal conductivity coefficient from a heat flux given on the boundary of the object was proposed. A problem of special interest is to compare the thermal conductivity coefficients obtained in two different formulations of the inverse problem, namely, µ(x, t) > 0, β(x, t) ≡ 0 (field functional) and µ(x, t) ≡ 0, β(x, t) > 0 (flux functional). The solution of the problem may be nonunique in the former case, whereas an analysis of numerous problems in the latter case has never revealed the nonuniqueness of the solution. Solving the same identification problem, but with different cost functionals showed that the thermal conductivity coefficient is determined more accurately in the case of the flux functional (although much more experimental data is used in the case of the field functional). An additional argument for using the flux functional in determining the thermal conductivity is that the heat flux on the boundary of the object is easier to measure than the temperature in the object. References [Dul’nev1974] G.N. Dul’nev and Yu.P. Zarichnyak. Thremal Conductivity of Mixtures of Composite Materials. Energiya, Leningrad, 1974. [in Russian] [Zverev2008] V.G. Zverev, V.D. Gol’din, and V.A. Nazarenko. Radiation-conduction heat transfer in fibrous heat-resistant insulation under thermal effect. High Temp., 46:108–114, 2008. [Alifanov2009] O.M. Alifanov and V.V. Cherepanov. Mathematical simulation of high-porosity fibrous materials and determination of their physical properties. High Temp., 47:438–447, 2009. [Kozdoba1982] L.A. Kozdoba and P.G. Krukovskii. Methods for Solving Inverse Thermal Transfer Problems. Naukova Dumka, Kiev, 1982. [in Russian] [Alifanov1988] O.M. Alifanov. Inverse Problems of Heat Transfer. Mashinostroenie, Moscow, 1988. [in Russian] [Vabishchevich1990] P.N. Vabishchevich and A.Yu. Denisenko. Numerical methods for solving inverse coefficient problems. Method of Mathematical Simulation and Computational Diagnostics, Mosk. Gos. Univ., Moscow, 35–45, 1990. 7 [Samarskii1997] A.A. Samarskii and P.N. Vabishchevich. Difference metods for solving inverse problems of mathematical physics. Fundamentals of Mathematical Simulation, Nauka, Moscow, 5–97, 1997. [Samarskii2003] A.A. Samarskii and P.N. Vabishchevich. Computational Heat Transfer. Editorial URSS, Moscow, 2003. [in Russian] [Marchuk1992] G.I. Marchuk. Adjoint Equation and the Analysis of Complex Systems. Nauka, Moscow, 1992. [in Russian] [Czel2012] Balazs Czel, Gyula Grof. Simultaneous Identification of Temperature-Dependent Thermal Properties via Enhanced Genetic Algorithm. Int J Thermophys, 33:1023–1041, 1979. [Zubov2016] V.I. Zubov. Application of fast automatic differentiation for solving the inverse coefficient problem for the heat equation. Comput. Math. Math. Phys., 56(10):1743–1757, 2016. [Albu2018] A.F. Albu and V.I. Zubov. Identification of the thermal conductivity coefficient using a given surface heat flux. Comput. Math. Math. Phys., 58(12):2031–2042, 2018. [Sirazetdinov1997] G.I. Marchuk. Optimization of Systems with Distributed Parameters. Nauka, Moscow, 1997. [in Russian] [Kraiko1979] A.N. Kraiko. Variational Problems in Fluid Dynamics. Nauka, Moscow, 1979. [in Russian] [Vasil’ev2002] F.P. Vasil’ev. Variational Problems in Fluid Dynamics. Nauka, Moscow, 1979. [in Russian] [Peaceman1955] D.W. Peaceman and H.H. Rachford. The Numerical Solution of Parabolic and Elliptic Differ- ential Equations. Journal of the Society for Industrial and Applied Mathematics, 3(1):28–41, 1955. [Samarskii2001] A.A. Samarskii. Theory of Finite Difference Schemes. Marcel Dekker, New York, 2001. [Douglas1956] D.W. Peaceman and H.H. Rachford. On the numerical solution of heat conduction problems in two and three space variables. Trans. Amer. Math. Soc., 82:421–439, 1956. [Thibault1985] Jules Thibault. Comparison of Nine Three-Dimensional Numerical Methods for the Solution of the Heat Diffusion Equation. Numerical Heat Transfer Fundamentals, 8(3):281–298, 1985. [Gao1996] Ch. Gao , Y. Wang. A general formulation of Peaceman and Rachford ADI method for the N- dimensional heat diffusion equation. Int. Comm. Heat Mass Transfer, 23(6):845–854, 1996. [Albu2018-2] A.F. Albu and V.I. Zubov. Identification of thermal conductivity coefficient using a given temper- ature field. Comput. Math. Math. Phys., 58(10):1585–1599, 2018. 8