A Generalization of Beck's Method for Inverse Heat Conduction Problems

A finite dimensional variational approach for Inverse Heat Conduction problems is studied. By means of the gradients and the Hessian of the minimizing functional, the unknown heat flux in form of a step function is characterized as the solution of a linear system of equations. In special cases, one obtains the well-known Beck method and a generalization thereof. The latter allows the determination of an optimal number of future times.


Introduction
Inverse Heat Conduction Problems (abbr.: IHCP) occur whenever surface temperatures and surface heat fluxes should be determined at inaccessible portions of the surface from corresponding measurements at accessible parts. These problems are known to be severely illposed. There are various approaches to regularize and stabilize the illposedness (see [4,11,20] and the references therein). An IHCP may be considered as a non-characteristic Cauchy problem for parabolic equations, the solvability is known only for special cases (see [25,12]). Stability estimates of Hölder type for non-characteristic Cauchy problems for the heat equation are given in [2,3,4,6,7,22,27] and for parabolic equations with coefficients depending smoothly on t in [8,23,26,28,29,30,35]. Some related results can also be found in [31,32].
In a series of papers we have investigated such problems from the analytical and numerical point of view (see [8] - [19], [37] - [40]). A powerful tool for illposed problem is the variational approach where the problem is considered as a control problem with an appropriate defect functional to be minimized. The Conjugate Gradient Method (CGM) can then be used to solve the minimization problem where the iteration index plays the role of a regularization parameter. The latter means that the iteration has to be stopped by an a-posteriori stopping criterion which, under certain assumption, yields a regularization method in the sense of Tikkonov ( [24,33,34]). The CGM we used in [14,17,18,40] has the interesting aspect that it is applicable to multidimensional IHCPs and, moreover, not only the surface heat flux can be determined but also the initial temperature or a coefficient in the parabolic equation or a source term.
In the engineering community the Beck method is well-known and successfully used since 40 years. However, the analysis of this method with respect to stability and regularization is not fully understood (see [1,18,37,38]). This method requires the knowledge of the initial temperature and proceeds sequentially. The latter means that the equations are only solved for a few "future" times while the heat flux (or temperature) is only taken up to the next time step, and the "sequential" procedure is started again by shifting it one time step further. The success of the method is based on the observation that the minimum of the defect functional is exactly obtained in every sequential step. This is true under the assumption that the desired heat flux (or temperature) is constant for all future times in every sequential step. Experimentally is has been observed that Beck's method is more stable the more future time steps are used.
In this paper we present the Inverse Heat Conduction Problem in a general setting and outline a variational approach in order to determine the desired temperature and heat flux at an inaccessible part of the boundary of the spatial domain from (Cauchy-) data of accessible parts. For the ease of presentation, we restrict ourselves to one-dimensional domains. For the minimizing functional in our variational approach, the gradient and second Fréchet derivative can be given by solutions of parabolic equations which are, in a certain sense, adjoint to the original problem. It is interesting to note that not only the boundary flux but also the initial temperature can be determined. In section 2 this variational approach will be specialized to finite dimensions when the desired heat flux is a step function. Here, we assume that the initial temperature is given. With the sensitivity function (see (16)), the finite dimensional solution fulfills a system of equation which is ill-conditioned. In section 3, a special case is considered where the desired heat flux is assumed to be constant over a certain number of time intervals. These leads to the well-known Beck method. We explain this method from the point of view of a variational approach. In the final section a generalization of Beck's method is studied where two fluxes are used -instead of one. We have introduced this generalization already in [18]. The most interesting aspect of this method is the fact that an optimal number of time intervals can be determined for which the pair of heat fluxes have to be calculated and for which temperature data are used. At the end we present these numbers for a model problem.

Inverse Heat Conduction Problems
Let ν ≥ µ be given positive constants, and a, c ∈ L ∞ (Q T ), f ∈ L 2 (Q T ) be given functions with ν ≥ a ≥ µ, Q T = (0, 1) × (0, T ], T > 0. Linear IHCPs lead us to the following problem: Find u, and in particular au x | x=1 , from the following non-characteristic Cauchy problem for linear parabolic equations Here ϕ, g ∈ L 2 (0, T ), and the solution can be understood in a weak sense namely as an element of cf. [18], 2.2 . Let us remark that the initial condition need not be given.
A solvability criterion for the problem (1) was first given by Holmgren ([25]) in 1904. Namely, for a necessary and sufficient condition for the existence of a solution of the noncharacteristic Cauchy problem (2) is that the function ψ(t) defined by  Several generalizations of this result for special cases of the problem (1)are given in [12,21]. However, solvability criteria for the general problem with nonsmooth coefficients are not known up to now.
Up to now there exist various methods for solving IHCPs. We mention here: 1) methods of establishing stability estimates and then using Tikhonov's regularization, 2) the sequential function specification method of Beck, 3) variational methods and iterative methods, 4) mollification methods, 5) perturbation methods, 6) the method of quasireversibility, and 7) function-theoretic methods. Comprehensive surveys of methods for solving IHCPs are given in [11,38]. As it has been noted, most of the available methods need a supplementary condition: either the initial condition, or a boundary condition.
We consider the problem (1) by variational methods; since the initial condition u| t=0 and the heat flux au x | x=1 are not given, we consider them as a control. We, therefore, reformulate (1) into the following optimal control problem: Minimize the functional is the control and u satisfies the following problem: The solution of (4) will be understood in a weak sense, namely in we can write the functional J(v) in (3) in operator form We note that A maps v : . It is quite interesting and important to investigate the properties of A, namely its domain, image, singular values and singular functions. One knows that the variational problem has still the ill-posed nature of the inverse problem (1) (cf. [18], 4.2). Moreover, (3), (4) has a solution for bounded convex sets V 0 , V 1 and every minimizing sequence converges weakly to the set of all its minimum points (see [18], 4.3).
To conclude this preliminary section, we provide results concerning the differentiability of J. We note that the mapping A defined in (5) is affine. Indeed by superposition one can write the solution of the direct problem (4) as denoting the three terms on the right-hand side by u 1 , u 2 , u 3 . The linear part of A is then given by i.e. the solution of (4) for g = 0 and F = 0. Hence, the shift part w of A = A L + w is given by w = u 3 = u(g, 0, 0, F ), i.e. the solution of (4) with v 0 = 0, v 1 = 0.
For an affine mapping A : B −→ H, B a Banachspace and H a Hilbert space, A = A L +w, with a bounded linear mapping A L ∈ L(B, H), the first and second (Fréchet-)derivative, resp. of a functional J(v) = 1 2 Av − ϕ 2 H are given by Here, A * L ∈ L(H, B ) is the adjoint operator to A L .
Specializing this to IHCPs, with B = H = L 2 (0, T ) one sees that J and J can be respresented as solutions of associated parabolic equations (see [18], 4.4), where ψ(x, t) = ψ(x, t; v) is the solution of the following problem adjoint to (4), where u = u(x, t; v) is the solution of Problem (4); in the representation of J , Ψ = Ψ(x, t; ∆v) is the solution of the problem with the solution ∆u = ∆u( One observes, that (10) and (11) for any Q ∈ L 2 (0, 1) × L 2 (0, T ) where A L is the linear part of A (see (7)).

Finite Dimensional Approximations
In the following we consider the case when the initial condition v 0 is given. Subdividing the time interval [0, T ] into N equidistant time steps, t i = iT /N, i = 0, . . . , N , we are looking for an approximation of the heat flux q := au x | x=0 at x = 0 by a step function Here, χ (t i ,t i+1 ] denotes the characteristic function of the interval (t i , t i+1 ]. By superposition, using (6) the affine mapping A in (5) can be represented as for t < t ≤ t +1 , = 0, 1, . . . , N − 1, The function s = s(x, x) is also called "sensitivity function". It is not a classical solution of the parabolic problem (15) but must be understood in the V 1,0 sense. For the sensitivity function, the heat flux at x = 0 is discontinuous in t = 0. One can define s(x, t) = 0, t < 0. For a given q, one can choose q i = q(t i+1/2 ) with t i+1/2 = 1 2 (t i+1 +t i ) in (15). This mapping may be called a "Neumann-to-Dirichlet" mapping. The linear part of A is obviously given by For an q as in (14), the minimizing functional (see (3)) in our variational approach is given by dt. (17) with the gradient ∇J( q) = ∂J/∂q k ( q) k=0,...,N −1 , where r = A q − ϕ denotes the residual (or defect).
The restriction of A and A L with values in R N will be denoted by A N , A N,L , where δ t s(x, τ ) := s(x, τ + ∆t) − s(x, τ ), (∆t := T /N ) denote forward differences w.r.t. the time variable. Instead of q we will use the notation q for the vector (q 0 , . . . , q N −1 ) determining the step function q. The integrals in the representation of J and its gradient will simply approximated by the trapezoidal rule which gives where r l = (A q − ϕ)(t l+1/2 ), l = 0, . . . , N − 1. From this, it is not difficult to determine the second partial derivatives of J N , We note that the Hessian H J of second derivatives is symmetric and has diagonal entries Furthermore, it is independent of q itself -or, in other words, H J is constant (in q). We finally like to mention, that, in the differences δ t s for j = N (cf. (20) and (21)) the term s(0, t −1/2 ) vanishes since s(·, t) = 0 for t < 0; hence δ t s(0, t −1/2 ) = s(0, t 1/2 ).
In order to minimize J N (q) we are looking for a "stationary point" of J N , i.e. a vector q = (q 0 , . . . , q N −1 ) which solves ∇J N (q) = 0. The crucial observation is now due to the fact that H J is independent of q itself, and Taylor's formula up to second derivatives yields the identity.
Thus, using an initial approximation q (0) , the optimal (stationary) q = q (0) + ∆q may be obtained by the increment provided H J is regular. This is nothing else than the first step of Newton's method. In the present situation we observe that the first step of Newton's method already yields the desired (stationary) solution. Solving (23) is obviously equivalent to the solution of Hence, the solution of (23) means the inversion of A * N,L A N,L which is known to be an extremely ill-conditioned system of linear equations. In the other words, the Hessian H J represents A * N,L A N,L and its condition number reflects the ill-posedness of the underlying problem.
For large conditon numbers of H J it is not recommended to solve equation (24) directly but, instead, to use the conjugate gradient method (cf. [18], 5.3). To check the condition numbers, of course one has to calculate the values and difference of the sensitivity function s. Since the latter has a discontinuity for its flux at t = 0, smaller time increments are recommended for a certain number of initial steps.

Constant Heat Flux Case (Beck's Method)
We consider the very special case that the q i 's are required to coincide for all intervals (t i−1 , t i ] , i = 1, . . . , N . In this situation, one wants to find one single value q such that where r j = u(0, t j−1/2 ; q) − ϕ(t j−1/2 ) , t ∈ (0, T ). In this special case, we have This follows from the representation and, thus, which yields (25). One step of Newton's method again yields the exact stationary solution, i. e. J N (q) = 0, which, in this situation, has the special form where q * denotes an initial approximation to q andr * j = u(0, t j−1/2 ; q * ) − ϕ(t j−1/2 ). Moreover, q provides a minimum of J N because J N > 0.
Comparing (28)  Here, A N,L q can be considered as a row vector (or a N × 1-matrix), A * N,L is corresponding adjoint 1 × N -matrix (or column vector) defined by with the discrete L 2 -scalar product The fundamental equation (23) -or (24) -then yields formula (28).
Formula (28) gives exactly the well-known Beck's method for N "future times" when the initial temperature distribution is given. Beck's method in its original form is used in a sequential way which means that, for a given temperature distribution u (k) 0 and heat flux q (k) (at x = 1) at time t k , the heat flux q (k+1) in the time interval (t k , t k+r ] is given by formula (28) with N = r future times and initial approximation (q * =)q (k) . We note also that u = u (k+1) depends on the sequential step. The new heat flux is then taken only in (t k , t k+1 ], and the (sequential) procedure is started again at t k+1 with initial function u (k+1) 0 := u (k+1) (·, t k+1 ; q (k+1) ). It is remarkable that this method yields a stationary solution in the constant heat flux case which implies that no further iteration is needed.
However, the initial function u (k+1) 0 for the next sequential step is not chosen in an optimal way but by some heuristic reasoning. (Instead of using an 'update' u (k+1) 0 = u (k+1) ·, t k+1 ; q (k+1) as the initial function for the next time step, one could use u The latter might be an explanation why Beck's method is difficult to analyse from the viewpoint of stability of a time stepping algorithm. On the other hand, the fact that it always meets a stationary point -at least in the special case considered -might give an explanation why Beck's method works so well and is successfully applied by many engineers over the last 40 years.
The coefficients k j = s(0, t j−1/2 )/D which multiply the residualsr * j in (28) to obtain the new q are called "gain coefficients" in the book of Beck et al. (see [1], chapt. 4.4.3). They provide a measure for stability of the time stepping algorithm to determine q = q (k+1) from q * = q (k) .

Two Fluxes Case
A generalization of Becks's method is obtained when one wants to find an optimal pair of heat fluxes (q,q) such that In this case the gradient vector ∇J N = (∂J N /∂q, ∂J N /∂q) has two components and the Hessian of second derivatives is a 2 × 2-matrix, It is not difficult to see that det H J > 0 ; here, the determinant of the Hessian is given by With an initial approximation q * = (q * ,q * ), a stationary solution q = q * + ∆q , ∆q = (∆q, ∆q) is given by the solution of (23). In the present situation, the first component of (23) solution may be written down explicitly in the following way, If we apply this method in a sequential way, then, with an initial approximation u (k) 0 and heat flux q (k) for time t k , we only use ∆q = ∆q (k+1) from formula (32) and obtain q (k+1) = q (k) +∆q (k+1) with N = r future times andr (k) j :=r * j := u (k+1) (0, t k+j−1/2 ; q (k) )− ϕ(t k+j−1/2 ) , j = 1, . . . , r. One may call this sequential procedure a "generalized Beck method".
For the model problem u t = u xx + F , the following figure presents the condition numbers of the corresponding 2 × 2-matrices H J for various ∆t and N (= r) in a logarithmic scale. For the mesh width in the spatial interval [0, 1] we have chosen ∆x = 1/8; we did not observe a great difference to corresponding condition numbers for ∆x = 1/10 or ∆x = 1/16. It can be seen that the inversion of H J is stable as long as the condition numbers remain below, say, 2 * 10 2 . As one expects, for a fixed r and for growing ∆t s, the condition numbers increase and the method becomes more stable.
The most interesting observation is that for any (fixed) ∆t one can find a minimum of the condition numbers for a specific r = r min given by the following Increasing r = 3, 4, . . . one observes that the condition numbers decrease up to r min ; and, thus, the method becomes more stable; if one increases r further the condition numbers grow. Another observation is that the r min 's are of moderate size. Consequently, using these r's in a sequential procedure the computational overhead is acceptable.
We suggest to use these optimal number of future times not only for the generalized Beck method but also for the original one. If more than one measurement point is available, one hat to sum over these x−positions in (32) with ths corresponding x−values for the sensitivity function.
In a forthcoming, more experimental work, we will compare the computational performance of the original and generalized Beck method.

Conclusions
We have studied a finite dimensional, variational approach for IHCP's where the initial temperature is given. If the sensitivity function is known with a sufficiently high accuracy, the minimizing functional, its gradient and matrix of second derivatives is explicitly given. By means of the Hessian, the solution -i.e. the temperature and flux at the inaccessible part of the spatial domain -can be obtained for all discrete times which is called a "whole domain method".
However, the corresponding system of equations is extremely ill-conditioned and has to be solved iteratively or by a sequential procedure. For the latter, one solves the underlying parabolic equation for a few "future time steps" assuming the desired heat flux to be constant during these steps. This is known as the "Beck method" which is explained here from the viewpoint of a variational approach. An important feature of this method is the fact that the optimal (constant) heat flux can be explicitly determined by only one step in Newton's method.
A generalization of Beck's method uses two fluxes for the future time steps and obtains the optimal pair of fluxes again by an explicit formula (after one Newton step). For this method one observes that for any (fixed) time step size there exists an optimal number of future times for the which the condition numbers of the corresponding Hessian -2 × 2matrix -are minimal. This answers the question how to chose the number of future times related to the size of the time step width at least for this variant of Beck's method.