ESTIMATING EIGENVALUES OF AN ANISOTROPIC THERMAL TENSOR FROM TRANSIENT THERMAL PROBE MEASUREMENTS

We propose a new method for estimating the eigenvalues of the thermal tensor of an anisotropically heat-conducting material, from transient thermal probe measurements of a heated thin cylinder. We assume the principal axes of the thermal tensor to have been identified, and that the cylinder is oriented parallel to one of these axes (but we outline what is needed to overcome this limitation). The method involves estimating the first two Dirichlet eigenvalues (exponential decay rates) from transient thermal probe data. These implicitly determine the thermal diffusion coefficients (thermal tensor eigenvalues) in the directions of the other two axes. The process is repeated two more times with cylinders parallel to each of the remaining axes. The method is tested by simulating a temperature probe time-series (obtained by solving the anisotropic heat equation numerically) and comparing the computed thermal tensor eigenvalues with their true values. The results are generally accurate to less than 1% error.

1. Introduction. This paper proposes a new method for estimating the eigenvalues of the thermal tensor of an anisotropically conducting material (for example, a nanocomposite) from transient thermal decay data.
Techniques for estimating thermal diffusion coefficients have been largely of three types: • Those making use of a steady state distribution of temperature (see, for example, [7]). • Those making use of a transient heat flow (see, for example, [2,3]).
• Exact analytic formulations, e.g., "final overdetermination" ( [5, p. 218]). The first category uses a "divided bar" technique: two samples ("bars") are exposed to the same heat source. One is the material with unknown thermal diffusion coefficient and the other ("reference bar") has known parameters. Comparison at equilibrium yields an estimate for the unknown coefficient. See [7].

STEVE ROSENCRANS, XUEFENG WANG AND SHAN ZHAO
Our approach is of the transient category. But previous procedures in that category have been confined to unbounded geometries for which there are explicit formulas (see [4]) for the transient temperature flow. An example is [2], which uses the half-space x 3 ≥ 0 in 3D. These explicit formulas, when compared with the measured temperature, enable diffusion coefficients to be estimated. Our experimental domain will be a bounded region. Unbounded domains (such as a half-space) can of course only be approximated in a real experiment, and this truncation would be an additional source of error in such treatments. Our proposed method is new and remains to be tested, but we do provide evidence (in Section 6) that suggests it would provide accurate recovery of diffusion coefficients.
Consider a bounded 3D domain Ω composed of the material. Suppose homogeneous Dirichlet/Neumann boundary conditions and that at time t = 0 the temperature is assigned everywhere on Ω, say u = u 0 (x 1 , x 2 , x 3 ). The temperature u at later times is the solution to the anisotropic heat equation in which the symmetric positive-definite matrix A = (a ij ) (assumed constant) is the so-called thermal tensor. After a rotation to its principal coordinate system, the matrix becomes diagonal, with entries equal to the eigenvalues of A. We assume that the three (mutually orthogonal) principal axes of the thermal tensor A have been identified ( [2] also begins with this assumption), and in effect where the eigenvalues a > 0, b > 0, and c > 0 are unknown. The heat equation in these principal coordinates is For for some coefficients {c n }, and distinct exponential decay rates (distinct Dirichlet eigenvalues, see Section 2) The eigenspace of λ 1 is simple. In this generality the eigenspace of λ 2 may not be simple, but in our main application it is. Provided the sampling point (x 0 , y 0 , z 0 ) ∈ Ω and initial data are acceptable (this will be clarified in Section 2), knowledge of u(x 0 , y 0 , z 0 , t) for t ≥ 0 implicitly determines the decay rates λ 1 , λ 2 , . . .. At any acceptable sampling point the time-series of temperature measurements (0 < t 1 < t 2 < . . . < t N ) can be used to estimate (several terms of) the λ-sequence, e.g., the first three terms. 1 In the original x 1 , x 2 , x 3 coordinates we would have u( The λ's are the same as in (3) because they are rotational invariants.
If λ 1 , λ 2 , λ 3 could be accurately estimated, the recovery of the thermal tensor eigenvalues a, b, c would be a simple matter (see Appendix A). However, the curvefitting on the basis of the time-series data (4) is an ill-posed problem, see [8]. The accurate estimation of three λ's seems to be out of the question. We take Ω to be a thin cylinder with axis parallel to one of the principal axes of the thermal tensor. We will then (in Section 2) outline a procedure requiring in effect only the estimation of λ 1 and λ 2 , but for two different initial data, so that the estimation of only one dominant exponential decay coefficient (at a time) is required. This permits much greater precision in our estimates. In Section 3 we carry out the estimation of λ 1 and λ 2 from the time-series data u A and u B specified in Section 2.
In Section 4 we show how the estimates of λ 1 and λ 2 lead to estimates of the eigenvalues a and b of the thermal tensor A, the first two coefficients in equation (2). The method requires knowledge of the functions →λ i ( ), 0 < ≤ 1, i = 1, 2 whereλ 1 ( ) andλ 2 ( ) are the principal and second Dirichlet eigenvalue of the elliptic operator − 2 ∂ 2 ∂x 2 − ∂ 2 ∂y 2 on the unit disk. These are tabulated in Appendix B. Section 5 shows how to estimate c, the last of the three thermal tensor eigenvalues. In Section 6 we verify the accuracy of our method by simulating experimental temperature probe data for a 2D disk with given diagonal thermal tensor by numerically solving (using Matlab's Partial Differential Equation Toolbox) the anisotropic heat equation (2) with appropriate initial data and Dirichlet boundary conditions (and, optionally, adding small random noise to the data and solution). This gives us a time-series to which we apply the methods outlined above, obtaining estimates a est and b est . These generally have error less than 1%. Matlab scripts for these calculations are included in Appendices C, D, E, and F.
2. Thin cylinder. Initial data. Now we specify the domain Ω to be a cylinder whose axis is the z-axis: with Dirichlet conditions (u = 0) on the vertical sides and Neumann conditions (∂u/∂z = 0) on the top and the bottom. 2 Since the domain is the cartesian product of the 2D disk of radius R and the 1D interval 0 < z < h, the eigenfunctions of Ω are of the form V · W , where (after a rotation about the z-axis) V is a (normalized) Dirichlet eigenfunction of the diagonal operator −a ∂ 2 ∂x 2 − b ∂ 2 ∂y 2 on the same disk (assume a ≤ b) with eigenvalues λ 1 , λ 2 , . . . and W is a Neumann eigenfunction of the interval 0 < z < h, with eigenvalues cj 2 π 2 /h 2 . The eigenvalues of the full problem are of the form Our intuition tells us that if the cylinder is thin enough (h small enough) the situation is closer to 2D. Indeed, the following are eigenvalues not necessarily in numerical order: Assume from now on that h is small enough, then and consequently the first two (Dirichlet/Neumann) eigenvalues of the thin cylinder are the same as the first two Dirichlet eigenvalues of the 2D disk. 3 We emphasize that λ 1 , λ 2 , etc., are invariant under rotations about the z-axis, just as is the cylinder itself. The first two corresponding normalized eigenfunctions of the cylinder With initial data independent of z (which we assume), c 1 and c 2 in the expansion (3) are independent of z and are given by the 2D formulas (in this situation the eigenspace of λ 2 is simple) The principal eigenfunction V 1 is positive in the interior of the disk, whereas (because b > a, which we assume) the y-axis is a nodal line for V 2 , and in fact it can be shown that V 2 is an odd function of x and V 1 is an even function of both x and y.
We use two initial data (both for the initial-value problem to be used to simulate experimental data, and also for the experiment itself).
Data u 0A : initial data even in x and hence orthogonal to . Data u 0B : initial data odd in x and hence orthogonal to V 1 , u 0B = x · u 0A . The solution to the heat equation u B is of the form d 2 e −λ2t +d 3 e −λ3t +O(e −λ4t ) ∼ d 2 e −λ2t (t → ∞).
Thus to estimate λ 1 and λ 2 it is only necessary to estimate the leading exponential decay coefficient for each of the two time-series u A and u B . We have not yet discussed the choice of sampling point (for the temperature probe).
• It can be on the top or bottom of the cylinder.
• It should not be near the lateral boundary, for then all terms will be very small. • In the case of u B , it should not be near the y-axis. This would make d 2 small, where u B ≈ d 2 e −λ2t + d 3 e −λ3t + . . . and this would delay the dominance of the first term. • Apart from this, the method is not sensitive to the sampling point. The knowledge of all three principal axes was only used in the specification of odd and even data, just above. If only one principal axis were known (call it the z-axis), we would not have knowledge of data satisfying the above orthogonality conditions and the estimation of λ 1 and λ 2 from thermal probe data would be much less accurate. After estimating λ 1 and λ 2 , we would pass to Section 4. 4 3. Estimation of Dirichlet eigenvalues λ 1 and λ 2 from the two time-series of temperature probe data, u A and u B . Following a suggestion of Mac Hyman (personal communication), we consider the ratio u(t + δ)/u(t), δ ∼ O(1) not yet determined. From the above it easily follows that which suggests simple estimates for the leading exponential decay coefficients λ 1 and λ 2 , from the u A and u B time-series respectively.
But the time-series of temperature measurements is unlikely to be accurate for such large times; the temperature has decayed considerably. The same seems to be true of the simulated time-series discussed in Section 6: evaluating the right sides of the last two equations (15) and (16) for large t does not yield correct results. We proceed therefore to use these estimates for smaller t, with no guarantee that they will work. We take u A (t + δ) and u A (t) to be time-adjacent points in the u A time-series, and we do the same for u B , i.e., δ is the time between two successive measurements. It may change from one pair of measurements to another. In the case of the simulated experimental data, there are generally 81 observations, hence as many as 80 separate estimates for each of λ 1 and λ 2 .
In every calculation we have made (i.e., Section 6 with various values of the thermal tensor eigenvalues a and b), we have found that the 80 estimates of the points (λ 1 , λ 2 ) display a pronounced clustering in the plane: many outliers but also many points very close together. There is always one and only one cluster. See Figures 1, 2, 3.
Our algorithm (Matlab script presented in Appendices D and F) for eliminating outliers in the (λ 1 , λ 2 ) plane is: • Drop all data points whose λ 1 or λ 2 z-score 5 is greater than 1 in absolute value. • Repeat the process until only one or two data points remain.
-If one point remains it is taken as best estimate.
-If two points remain the mean of the two is taken as best estimate. We compared the above procedure with • Exponential curve-fitting using Matlab's Curve Fitting Toolbox (fitting 'exp2').
• Introduction of a single time-series u C , arising from initial data not orthogonal to either the principal or second eigenfunction and application of the above Toolbox to find the best fit to u C of the form • Per Sundqvist's script 'exp2fit' [9]. None of these three approached the accuracy obtained by using data u 0A and u 0B and (13) and (14). 4. Recovery of thermal tensor eigenvalues a and b from the Dirichlet eigenvalues λ 1 and λ 2 . Recall that the z-axis is a principal axis for the thermal tensor A, hence A is similar to is a 2×2 positive-definite matrix, whose eigenvalues we have called a and b, with = a/b ∈ (0, 1). A rotation to principal coordinates, (x 1 , x 2 ) → (x, y) and a rescaling of the independent variables yields the following: The last characterization (Dirichlet eigenvalues of minus Laplacian on an ellipse) will not be used, but is included for completeness.
Let λ 1 , λ 2 andλ 1 ( ),λ 2 ( ) be the principal and second Dirichlet eigenvalues for the respective situations. We have computedλ 1 ( ),λ 2 ( ) in steps of ∆ = .025. These results are included in Appendix B. 6 Assuming λ 1 , λ 2 to have been estimated, the relation λ = bλ( )/R 2 impliesλ We solve this equation for . For this step table-lookup and interpolation (Matlab's interp1 ) are used. Finally, requiring another table-lookup) and This completes the recovery of the thermal tensor eigenvalues a and b from thermal decay measurements.

Remark 1.
We have used a simple re-scaling. The argument does not assume any eigenspace dimensions and does not even use the assumption that λ 1 and λ 2 are the first and second distinct Dirichlet eigenvalues. They could be, for example, the fifth and eighth distinct Dirichlet eigenvalues. (Thenλ 1 ( ) andλ 2 ( ) would be the fifth and eighth distinct Dirichlet eigenvalues in the scaled situation.) The Matlab script for these steps is included in Appendix E. It uses the tabular dataλ 1 ( ),λ 2 ( ), ratio12 in steps of ∆ = .025. (These arrays are denoted lam1shan, lam2shan and epsshan. The array ratio12 is lam1shan./lam2shan.) See the table in Appendix B.
The two table-lookups (requiring the publication and dissemination of tables!) and the interpolation can be avoided (at the expense of a slight increase in relative error) by explicit curve-fitting of the two implicitly defined functions λ 1 /λ 2 → 1/λ 1 ( ) and λ 1 /λ 2 → 2 /λ 1 ( ), where is the solution of Equation (17). Matlab's Curve-Fitting Toolbox gives explicit approximations for a and b in terms of λ 1 and γ :≡ λ 1 /λ 2 : 5. Recovery of the third thermal tensor eigenvalue c. We now repeat the process with a new thin cylinder whose axis is parallel to, say, the x-axis. Our procedure will estimate the thermal tensor eigenvalues b and c. Finally, take a last thin cylinder whose axis is parallel to the y-axis, and thus estimate a and c. Altogether, we have two estimates for each of the three thermal tensor eigenvalues. They should agree.
6. Simulated temperature probe data; numerical results. Matlab's PDE Toolbox was used to create the time-series The unit disk and its boundary are selected in the graphical user interface (gui). The mesh is initialized and refined twice. At this point the arrays b (boundary), p (nodes), e (edges) and t (triangles) are exported to the Matlab workspace, to allow the PDE to be solved using command-line functions rather than the gui. In the m-file q1 quoted in Appendix C, the command load(DISK_BPET.mat) loads all these arrays.
The command-line function parabolic outputs the solution at a finite sequence of times tlist, which we have taken to be 0, h, 2h, 3h, . . . , 2 plus a random vector 7 of the same length with mean h/2. We allow the solution and data to be perturbed by a small normal relative noise. These random effects are meant to simulate real experimental effects.
To select h appropriately we must consider the order of magnitude of the thermal tensor eigenvalues a and b. Recall b > a. If = a/b << 1, the problem is a singular perturbation, and can't generally be handled by our methods. (But even, for example, b = 1, a = .1 is successfully computed.) If a and b are roughly of order 1, we take h = .025. Much smaller or much larger values for h would degrade the results. But the scaling of (2) shows that if a and b are (roughly) of order 10, we should replace h → h/10. If a and b are roughly of order .1, then we should set h → 10h ≈ .25, etc. The tlist variable has to be altered accordingly. So in a real experiment, while a and b are unknown, we do need to have some rough a priori estimate of their order of magnitude.