Wavelet and Fourier Methods for Solving the Sideways Heat Equation.
We consider an inverse heat conduction problem, the  sideways heat equation, which is a model of a problem, where one wants to determine the temperature on both sides of a thick wall, but where one side is inaccessible to measurements. Mathematically it is formulated as a Cauchy problem for the heat equation in a quarter plane, with data given along the line x=1, where the solution is wanted for $0 \leq x < 1$.The problem is ill-posed, in the sense that the solution (if it exists) does not depend continuously on the data. We consider stabilizations based on  replacing the time derivative in the heat equation by  wavelet-based approximations or a Fourier-based approximation. The resulting problem is an initial value problem for an ordinary differential equation, which can be solved by standard numerical methods, e.g., a Runge--Kutta method.We discuss the numerical implementation of Fourier and wavelet methods for solving the sideways heat equation. Theory predicts that the Fourier method and a method based on Meyer wavelets will give equally good results. Our numerical experiments indicate that also a method based on Daubechies wavelets gives comparable accuracy. As test problems we take model equations with constant and variable coefficients. We also solve a problem  from an industrial application with actual measured data.
Introduction
In many industrial applications one wishes to determine the temperature
on the surface of a body, where the surface itself is inaccessible for measurements
[1]. It may also be the case that locating a measurement device
(e.g. a thermocouple) on the surface would disturb the measurements so
that an incorrect temperature is recorded. In such cases one is restricted
to internal measurements, and from these one wants to compute the surface
temperature.
In a one-dimensional setting, assuming that the body is large, this situation
can be modeled as the following ill-posed problem for the heat equation
in the quarter plane: Determine the temperature u(x; t) for
temperature measurements g(\Delta) := u(1; \Delta), when u(x; t) satisfies
Of course, since g is assumed to be measured, there will be measurement
errors, and we would actually have as data some function g m 2 L 2 (R), for
which
where the constant ffl ? 0 represents a bound on the measurement error. For
can solve a well-posed quarter plane problem using g m as data.
For we have the sideways heat equation 1
Note that, although we seek to recover u only for 0 - x ! 1, the problem
specification includes the heat equation for x ? 1 together with the
boundedness at infinity. Since we can obtain u for x ? 1, also u x (1; \Delta) is de-
termined. Thus we can consider (1.3) as a Cauchy problem with appropriate
Cauchy data [u; u x ] given on the line
Although in this paper we mostly discuss the heat equation in its simplest
form u our interest is in numerical methods that can be used for
more general problems, e.g. equations with non-constant coefficients,
Also referred to as the Inverse Heat Conduction Problem (IHCP) [1].
or nonlinear problems,
which occur in applications. For such problems one cannot use methods
based on reformulating the problem as an integral equation of the first kind,
since the kernel function k(t) is explicitly known
only in the constant coefficient case. Instead, we propose to keep the problem
in the differential equation form and solve it essentially as an initial value
problem in the space variable ("space-marching", [1, 5, 6, 11, 12, 13, 22, 23]).
In [11, 12] we have shown that it is possible to implement space-marching
methods efficiently, if the time-derivative is approximated by a bounded op-
erator. In this way we obtain a well-posed initial-boundary value problem
for an ordinary differential equation (ODE). The initial value problem can
be solved by a standard ODE solver, see Section 2.1.
In this paper, we study three different methods for approximating the
time derivative and their numerical implementation. First, the time derivative
is approximated by a matrix representing differentiation of trigonometric
interpolants, see e. g. [14, Sec. 1.4]. In Section 3 we first give some error
estimates (for the non-discrete case, see also [24]), and then we discuss the
numerical implementation of this method.
Meyer wavelets have the property that their Fourier transform has compact
support. This means that they can be used to prevent high frequency
noise from destroying the solution 2 . In a previous paper [25] we studied
the approximation of the time derivative in (1.3) using Meyer wavelets, and
gave almost optimal error estimates. Here we discuss the numerical implementation
of this method, and, in particular, we show how to compute the
representation of the time derivative in wavelet basis.
Daubechies wavelets themselves have compact support and, consequently,
they cannot have a compactly supported Fourier transform. However, the
Fourier transform decays so fast that in practice they can be used for solving
the sideways heat equation in much the same way as Meyer wavelets. We
discuss this in Section 4.5.
We emphasize that although all theoretical results are for the model
problem with constant coefficients, we only deal with methods that can
be used for the more general problems (1.4)-(1.5). Also, the quarter plane
assumption is made mainly for deriving theoretical results; it is not essential
in the computations (provided that u x (1; \Delta) can be obtained, either from
measurements of by assuming symmetry, cf. Section 5.2). A constructed
2 It is very common in ill-posed problems that the ill-posedness manifests itself in the
blow-up of high frequency perturbations of the data.
variable coefficient example is given in Section 5. There we also present a
problem from an industrial application, with actual measured data.
Ill-posedness and Stabilization
The problem of solving the sideways heat equation is ill-posed in the sense
that the solution, if it exists, does not depend continously on the data. The
ill-posedness can be seen by solving the problem in the Fourier domain. In
order to simplify the analysis we define all functions to be zero for t ! 0.
Let
g(t)e \Gammai-t dt;
be the Fourier transform of the exact data function 3 . The problem (1.1)
can now be formulated, in frequency space, as follows:
bounded .
The solution to this problem, in frequency space, is given by
where
i- denotes the principal value of the square root,
ae
In order to obtain this solution we have used the bound on the solution at
infinity. Since the real part of
i- is positive and our solution b u(x; -) is
assumed to be in L 2 (R), we see that the exact data function, b g(-), must
decay rapidly as - !1.
Now, assume that the measured data function satisfies g m
(R), is a small measurement error. If we try to solve the
problem using g m as data we get a solution
Since we can not expect the error b ffi(-) to have the same decay in frequency
as the exact data bg(-) the solution bv(x; -) will not, in general, be in L 2 (R).
3 In this paper all Fourier transforms are with respect to the time variable.
Thus, if we try to solve the problem (1.3) numerically, high frequency components
in the error, ffi , are magnified and can destroy the solution. However,
if we impose an apriori bound on the solution at and in addition,
allow for some imprecision in the matching of the data, i.e. we consider the
then we have stability in the folllowing sense: any two solutions of (2.4), u 1
and
For this reason we call (2.4) the stabilized problem. The inequality (2.5)
is sharp and therefore we can not expect to find a numerical method, for
approximating solutions of (2.4) that satisfy a better error estimate.
Different methods for approximating solutions of (2.4) exist. One difficulty
is that for (2.4) we do not have uniqueness, so in order to get a
procedure that can be implemented numerically, it is necessary to somehow
modify the problem. Often, the dependence on ffl and M is included by
choosing the value of some parameter in the numerical procedure.
2.1 Stabilization by Approximating the Time Derivative
In this subsection we consider the "initial-value" problem
' u
x
@
'' u
with initial-boundary values 4
The initial values for u x can be obtained (in principle and numerically) by
solving the well-posed quarter plane problem with data u(1;
x - 1.
4 For a discussion of the problem of setting numerical boundary values at
to [12].
We can write the solution of (2.6)-(2.8) formally as
' u
hm (t)
@
Loosely speaking, since the operator B is unbounded, with unbounded eigen-values
in the left half plane, the solution operator
unbounded frequency noise in the data can be blown up and destroy
the numerical solution. Even if the data are filtered [4, 23], so that high
frequency perturbations are removed, the problem is still ill-posed: rounding
errors introduced in the numerical solution will be magnified and will
make the accuracy of the numerical solution deteriorate as the initial value
problem is integrated.
In a series of papers [11, 12, 25], we have investigated methods for solving
numerically the sideways heat equation, where we have replaced the operator
@=@t by a bounded operator. Thus we discretized the problem in time, using
differences [11] or wavelets [25], so that we obtained an initial value problem
U x
x
'' U
U x
where the matrix D is a discretization of the time derivative,
and U are semi-discrete representations of the solution and its
derivative, and Gm and Hm are vectors. Thus, (2.9) can be considered as a
method of lines.
The observation made in [11, 12] is that when D represents a discrete
approximation of the time derivative, then D is a bounded operator, and the
problem of solving (2.9) is a well-posed initial value problem for an ODE. In
[11, 25] error estimates are given for difference and wavelet approximations.
The coarseness of the discretization is chosen depending on some knowledge
about ffl and M (actually their ratio). For solving (2.9) numerically, we can
use a standard ODE solver. In [12] we show that, in the case of a finite
difference approximation, it is often sufficient to use an explicit method, e.g.
a Runge-Kutta code.
In the rest of this paper, we will discuss error estimates and numerical
implementation of approximations of (2.4) by discretized problems of the
type (2.9), by a Fourier (spectral) method, and by Meyer and Daubechies
wavelets.
5 Cf. (2.2), and also [12], where a discretized equation is considered.
3 A Fourier Method
Here we consider how to stabilize the sideways heat equation by cutting off
high frequencies in the Fourier space. Error estimates are given in Section
3.1. Similar results can be found in [24]. In Section 3.2 we then arque that
in order to useful for more general equations (1.4)-(1.5), the Fourier method
should be implemented as in (2.9).
We start with the family of problems in Fourier space,
bounded .
parameterized by -. In Section 2 the solution was shown to be
e
Since the principal value of
i- has a positive real part, small errors in high
frequency components can blow up and completely destroy the solution for
natural way to stabilize the problem is to eliminate all high
frequencies from the solution and instead consider (3.1) only for
Then we get a regularized solution.
e
where - max is the characteristic function of the interval [\Gamma- In
the following sections we will derive an error estimate for the approximate
solution (3.3) and discuss how to compute it numerically.
3.1
In this section we derive a bound on the difference between the solutions
(3.2) and (3.3). We assume that we have an apriori bound on the solution,
. The relation between any two regularized solutions (3.3) is
given by the following lemma.
Lemma 3.1 Suppose that we have two regularized solutions v 1 and v 2 defined
by (3.3) with data g 1 and ffl. If we select
then we get the error bound
Proof: From the Parseval relation we get
\Gamma-
je
6 je
Using

From Lemma 3.1 we see that the solution defined by (3.3) depends continuously
on the data. Next we will investigate the difference between the
solutions (3.2) and (3.3) with the same exact data g(t).
Lemma 3.2 Let u and v be the solutions (3.2) and (3.3) with the same
exact data g, and let - Suppose that ku(0; \Delta)k 6 M .
Then
Proof: As in Lemma 3.1 we start with the Parseval relation, and using
the fact that the solutions coincide for - 2 [\Gamma-
Z
je
Z
je
Now we use the bound ku(0; \Delta)k 6 M , and as before we have -
which leads to the error bound

Now we are ready to formulate the main result of this section:
Theorem 3.3 Suppose that u is given by (3.2) with exact data g and that v
is given by (3.3) with measured data g m . If we have a bound ku(0; \Delta)k 6 M ,
and the measured function g m satisfies we choose
then we get the error bound
Proof: Let v 1 be the solution defined by (3.3) with exact data g. Then
by using the triangle inequality and the two previous lemmas we get

From Theorem 3.3 we find that (3.3) is an approximation of the exact
solution, u(x; t). The approximation error depends continously on the
measurement error and the error bound is optimal in the sense (2.5).
3.2 Numerical implementation of the Fourier-based method
Here we discuss how to compute the regularized solution (3.3) numerically.
Given a data vector Gm with measured samples from g m on a grid ft i g, a
simple method is to approximate the solution operator directly. Thus
we have the following algorithm:
1. b
Gm := FGm .
2. b
ae exp(
3. V (0; :) := F H b
Vm (0; :),
where F is the Fourier matrix. The product of F and a vector can be computed
using the Fast Fourier Transform (FFT) which leads to an efficient
way to compute the solution (3.3). When using the FFT algorithm we implicitly
assume that the vector Gm represents a periodic function. This is not
realistic in our application; and thus we need to modify the algorithm. We
refer to Appendix A for a discussion on how to make the problem 'periodic'.
The main disadvantage of this method is that it cannot be used for problems
with variable coefficients (1.4), (1.5). Note also that by approximating
the solution operator directly, we make explicit use of the assumption that
our solution domain is the whole quarterplane, t ? 0 and x ? 0.
An alternative, more widely applicable, approach is to approximate the
time derivative and use the ordinary differential equation formulation given
in Section 2.1, i.e. we take the matrix D in (2.9) to be
where   is a diagonal matrix which corresponds to differentiation of the
trigonometric interpolant, but where the frequency components with
are explicitly set to zero [14, Sec. 1.4]. This will filter the data and ensure
that we remove the influence from the high frequency part of the measurement
error in the solution. Thus in the ODE solver, multiplication by D F is
carried out as a FFT, followed by multiplication by   and finally an inverse
FFT. Multiplication with D F thus requires O(n log(n)) operations. We conclude
the section with the remark that this method can be considered as a
Galerkin method (cf. Section 4.2) with trigonometric interpolants as basis
and test functions.
4 Wavelet Methods
4.1 Multiresolution analysis
Wavelet bases are usually introduced using a multiresolution analysis (MRA)
(see Mallat [18], Meyer [21]). Consider a sequence of successive approximation
spaces,
and [
For a multiresolution analysis we also require that
This means that all spaces V j are scaled versions of the space V 0 . We also
require that there exists a scaling function, OE, such that the set fOE jk g j;k2Z ,
where OE jk is an orthogonal basis for V j . Since OE
we find that the function OE must satisfy a dilation relation
The wavelet function / is introduced as a generator of an orthogonal basis
of the orthogonal complement W j of V j in V j+1 (V
function / satisfies a relation similar to (4.2)
where the filter coefficients fg k g are uniquely determined by the fh k g. Any
function f 2 L 2 (R) can be written
l-j
where P j is the orthogonal projection onto V j ,
4.2 Meyer Wavelets and a Galerkin Approach
The Meyer wavelet / 2 C 1 (R)) is defined by its Fourier transform [7,
e \Gammai-t
where
sin
\Theta -j
cos
\Theta -j
0; otherwise,
and j is a C k or C 1 function satisfying
ae
It can be proved that the set of functions
for Zis an orthonormal basis of L 2 (R) [20, 7]. The corresponding
scaling function is defined by its Fourier transform,
cos
\Theta -j
0; otherwise.
The functions b
OE jk have compact support (see e.g. [25])
for any k 2 Z.
In [25] we presented a wavelet-Galerkin method, where starting from the
weak formulation of the differential equation,
(v
with test functions from V j , and with the Ansatz
c (-) (x)OE j- (t);
we get the infinite-dimensional system of ordinary differential equations for
the vector of coefficients c,
ae c
The initial values fl m are defined
and the matrix D j is given by
(D
The well-posedness of the Galerkin equation (4.8) follows from the fact that
the matrix D j is bounded with norm [25]
In [25] we proved the following stability estimate for the wavelet-Galerkin
method.
Theorem 4.1 Let g m be measured data, satisfying kg
that chosen so that
log M
Then the projection onto V j of the Galerkin solution v j+1 satisfies the error
estimate
Note that the result in Theorem 4.1 is suboptimal in two ways. Firstly,
we get slower rate of convergence when ffl tends to zero than in the optimal
estimate (2.5). Secondly, the error is not for v j itself, but only for the
projection of v j+1 onto V j .
4.3 Discrete Wavelet Transform
Discrete wavelet transforms can be defined in terms of discrete-time multiresolution
analyses (see e. g. [26, Section 3.3.3]). We will use DMT as a
short form of 'Discrete Meyer (wavelet) Transform'. Given a vector c 2 R n ,
is assumed to hold the coefficients
of some function v in terms of the basis of V J . This is the fine level
all contributions on yet finer levels are assumed to be equal
to zero. The DMT of c is equivalent to a matrix-vector multiplication [16,
p. 11]
where G J
n\Thetan is an orthogonal matrix. The subscript j indicates that
the part of the vector ~ c that contains the coarse level coefficients has 2 j
components; it holds the coefficients of the projection of v onto V j . We
illustrate how the transformation by G J
breaks up the vector into blocks of
coefficients corresponding to different coarseness levels in Figure 4.1.
~ v:

Figure

4.1: Schematic picture of the discrete wavelet transform ~
J \Gamma3 c.
The rightmost part of the vector represents wavelet coefficients on the resolution
level block coefficients on the level 2 J \Gamma2 , etc. The
leftmost part represents the coarse level coefficients.
Algorithms for implementing the discrete Meyer wavelet transform are
described in the thesis of Kolaczyk [16]. These algorithms are based on the
fast Fourier transform (FFT), and computing the DMT of a vector in R n
requires O(n log 2
n) operations [16, p. 66]. The algorithms presuppose that
the vector to be transformed represents a periodic function; we will return
to this question in Appendix A. A further illustration is given in Figure 4.2.
Here we have taken a data vector from R 256 , added a normally distributed
perturbation of variance 0:5 extended the vector smoothly to
size 512 in order to make it "periodic". Then we computed the DMT with
Thus the partitioning of the transformed vector is the same as in

Figure

4.1.
We see that the wavelet transform can be considered as a low pass filter:
the noise is removed and the leftmost part of the transformed vector is
a smoothed version of the data vector. The noise is represented in the
data vector

Figure

4.2: The upper graph shows a data vector and the lower graph its
DMT.
finer level coefficients. Note also the locality properties of the fine level
coefficients: only the first half of the data vector is noisy, and consequently,
each segment of fine level coefficients has (noticeable) contributions only in
its left half.
We remark that we use periodicity only for computational purposes,
since the codes (based on discrete Fourier transforms) are implemented for
this case. The vectors used are considered as finite portions of sequences in
4.4 Numerical Implementation
In the solution of the sideways heat equation in V j , we replace the infinite-dimensional
ODE (4.8), by the finite-dimensional
x
represents the approximation of the solution in
Gm and Hm are the projections of the data vectors on V j , and j is chosen
according to Theorem (4.1). For simplicity we suppress the dependence on
j in the notation for vectors. The matrix D d
j represents the differentiation
operator in V j , and since we are dealing with functions, for which only
a finite number of coefficients are non-zero, D d
j is a finite portion of the
infinite matrix D j (we use superscript d to indicate this).
In the numerical solution of (4.12) by an ODE solver, we need to evaluate
matrix-vector products D d
C. The representation of differentiation operators
in bases of compactly supported wavelets are described in the literature, see,
e.g. [2]. For such wavelets exact and explicit representations can be found.
Also finite difference operators can be represented easily [19]. In our context
of Meyer wavelets, which do not have compact support, the situation is
different. The proof of (4.10) in [25] actually gives a fast algorithm for this.
For \Gamma- t - define the function
Extend \Delta periodically, and expand it in Fourier series. In [25, Lemma A.2]
is is shown that
where d k is the element in diagonal k of D 0 . From the definition of D j it
is easily shown that D . Thus, we can compute approximations of
the elements of D j by first sampling equidistantly the function \Delta, and then
computing its discrete Fourier transform.
4.5 Daubechies Wavelets
The Daubechies wavelets have compact support [7] and therefore only a finite
number of filter coefficients (fh k g and fg k g in (4.2) and (4.3)) are nonzero.
The filter coefficients fg k g are uniquely determined by the coefficients fh k g
where L is the number of filter coefficients. In addition we want the wavelet
to have a number of vanishing moments, that is
Daubechies wavelets are defined in such a way that they have the largest

Figure

4.3: The Daubechies db4 scaling function (left) and wavelet (right).
possible number of vanishing moments given the number of non-zero filter
coefficients; in fact, we have
There is no explicit expression for the Daubechies wavelet, instead we
can compute it from the filter coefficients fh k g k2Z . For example the filter
coefficients associated with the Daubechies wavelet db4 are [7]
Since we are interested in using a basis of Daubechies wavelets for solving
the system of ODE:s (2.6) in a stable way, we need to find an approximation
of the derivative operator in the approximation spaces V j . Here we will use
the Galerkin approximation, D j , which is a Toeplitz matrix with elements
given by
(D
where the last equality follows from the definition of OE jk . This means that
it is sufficient to compute the derivative approximation on the space V 0 .
Since the scaling functions OE jk have compact support, the matrix D j will
be banded. For compactly supported wavelets it is possible to compute the
matrix D j explicitly. This result is due to Beylkin [2]. Since the matrix D 0
is constant along diagonals and only a few diagonals are non-zero, we can
insert the dilation relation (4.2) into (4.16) and we will get a small system of
linear equations to solve for the elements in D 0 . In the paper by Beylkin, the
elements of D 0 have been listed for different bases of Daubechies wavelets.
As in Section 4.2 we want to use a Galerkin approximation of the differential
equation. For the problem to be well-posed it is sufficient that kD j k
is bounded. In the case of periodized wavelets, the differentiation matrix
per
0 is a circulant and we have kD per
[27]. The same estimate can be
expected to hold for non-periodic wavelets.
The implementation details are similar to those in Section 4.4. We have
measured data g m (t) - u(1; t) on the finite interval
A ). Thus we have to approximate (2.6) with a finite-dimensional system
similar to (4.12) for the coarse level coefficients.
For compactly supported wavelets it is possible to give a fast implementation
of the wavelet transform by convolution with the filter coefficients
g. Since the filters are of short length we can compute the
wavelet transform of a vector in R n with O(n log(n)) operations[26].
When we solve numerically the discretized version of (2.6) we are interested
in computing the product of D j and a vector. Since D j is banded
and explicitly known this is simply a product with a sparse matrix and thus
requires
5 Numerical Experiments
5.1 Experiments
Here we present some numerical experiments to illustrate the properties of
the methods presented in the previous sections. We have solved discretized
versions of (2.6) using Matlab in IEEE double precision with unit roundoff
1:1\Delta10 \Gamma16 . The space marching was performed using a Runge-Kutta-Fehlberg
method (ode45 in Matlab) with automatic step size control, where the basic
method is of order 4 and the embedded method is of order 5. In all tests
the required accuracy in the R-K method was 10 \Gamma4 .
The tests were performed in the following way: First we selected a solution
computed data functions u(1;
and u x (1; by solving a well-posed quarter plane problem for the
heat equation using a finite difference scheme. Then we added a normally
distributed perturbation of variance 10 \Gamma3 to each data function, giving vectors
m and hm . Our error estimates use the signal-to-noise ratio M=ffl and
therefore we computed
kfk
From the perturbed data functions we to reconstructed u(0; t) and compared
the result with the known solution.
We conducted two tests:
Test 1: We solved the model problem (1.1) using a discontinous function,
f(t) as the exact solution.
Test 2: We solved the more general problem (1.4) using the coefficient
function
The results from these tests are given in Figure 5.1. In both cases the
length of the data vectors g m and hm were 1024. The regularization parameters
were selected according to the recipes given in Theorems 3.3 and
4.1. In the Meyer case we used projection on the space V 5 in both tests and
computed 64 coarse level coefficients. We have used the auxiliary function
in the definition of the scaling function
(4.5). In the Fourier case we used the 42 lowest frequency components when
calculating the time derivative. Since we have no stability theory for the
Daubechies method we chose to solve the problem in space V 5 in both cases.
Thus we computed 70 coefficients on the coarsest level. We used Daubechies
wavelets with filters of length 8 in the Galerkin formulation. In all cases the
number of steps in the ODE-solver were between 6 and 11. Before presenting
the results we recomputed our coarse level approximation on the finer
scale, using the inverse wavelet transform.

Figure

5.1: Solution of the sideways heat equation using the different meth-
ods. From top we have the results from the Fourier method, the Meyer
method and the Daubechies method. The results from Test 1 are to the left
and the results from Test 2 are to the right. The dashed line illustrates the
approximate solution, the solid line represents the exact solution and the
dash-dotted line represents the data function, g m .
5.2 An Industrial Application
In this section we present an example of an industrial problem where the
methods presented in the previous sections can be useful. The viability of
the methods is demonstrated by an experiment conducted in cooperation
with the Department of Mechanical Engineering, Link-oping University.
Consider a particle board, on which a thin lacquer coating is to be ap-
plied. In order to reduce the time for the lacquer coating to dry, the particle
board is initially heated. Since the temperature gradients on and close to
the surface of the board influence the drying time and the quality of the
lacquer coating, it is important to estimate the temperature and the temperature
gradients close to the surface. Often it is difficult or impossible
to measure the temperature directly on the surface of the board. Instead a
hole was drilled from the other side of the board and a thermocouple was
placed close to the surface, as seen in Figure 5.2. After the thermocouple
had been placed, the hole was filled in using the same material as in the
surrounding board. In the experiment a particle board, initially heated to

Figure

5.2: The cross-section, in principle, of the particle board used in the
experiment. The temperature g m is measured by a thermocouple, and we
seek to recover the temperature, f m , on the surface of the board [15].
suddenly placed to cool in air of room temperature. The thermocouple
placed inside the plate, at distance 2:9 mm from the surface,
gave the temperature history, g m (t), for From the
measured temperature we reconstructed the temperature history, f m (t), on
the surface of the board. The experiment presented here does not include
the lacquer coating. The heat equation in this case is given by
Here the constant - represents the physical properties of the problem. We
can formulate this as an initial value problem similar to (2.6). One difficulty
is that the surface temperature is not determined by the measured data g m
alone. We need also the temperature gradient at L. If the thickness of
the particle board is large in comparison to the distance between the surface
and the thermocouple, then it is reasonable to consider this to be a quarter
plane problem. In that case we can solve the heat equation for x ? L and
compute u x (L; \Delta). Another possibility is to assume that the temperature
function u(x; t) is symmetric with respect to the center of the board. In
this particular experiment this should be a very accurate assumption. In

Figure

5.3 we present results obtained using both these assumptions. The
computations were performed using the Meyer wavelet method presented in
Section 4.2. Initially the length of the vector g m was 8192; we solved only
for 128 coarse level coefficients. It is interesting to note that the solutions
almost coincide the first three minutes.
Temperature
Time min
Temperature
gradient
U x
(0,t)

Figure

5.3: The measured temperature vector, g m , sampled at 10Hz (solid)
and the corresponding surface temperature, f m from a quarter plane assumption
(dashed) and from a symmetry assumption (dash-dotted) In the
rightmost plot we see the temperature gradients at the thermocouple computed
from a quarter plane assumption (dashed) and from a symmetry assumption
(solid).
Time min
Temperature
Time min

Figure

5.4: Attempt to reconstruct the surface temperature, f m , using a
coarse level grid of size 512 (left) and 1024 (right), in the Wavelet-Galerkin
method. This means that we use less regularization than was the case in

Figure

5.3. Clearly these solutions are unphysical.
Since the temperature inside the particle board must be between the
initial temperature, 70 the temperature of the surrounding air, approximately
C, we can easily get a bound, M , on the solution. The noise
level, ffl, must be estimated from the measured data. It is often possible
to get a rough estimate by inspecting the data vector visually. In our error
estimates we use only the signal-to-noise ratio, ffl=M , and since we know only
rough estimates of M and ffl, we solved the problem several times using different
values for these parameters, i.e. different levels of regularization. By
experimenting, it is often easy to find an appropriate level of regularization.
In this particular experiment we know that the solution is monotonically de-
creasing, and this information can be used to rule out unphysical solutions,
see

Figure

5.4.
6 Concluding Remarks
We have considered three methods for solving the sideways heat equation,
based on approximating the time derivative by a bounded operator (matrix),
and then solving the problem in the space variable using a standard ODE
solver. For the Fourier and the Meyer wavelet methods, there is a stability
theory, which is used for choosing the level of approximation. In the case of
Daubechies wavelets the theory is incomplete.
From the theory for the Fourier and Meyer wavelet methods one can
expect that these methods give results of comparable accuracy. Considering
the fast decay in frequency of Daubechies wavelets, one can hope that this
method is about as accurate as the other two. Our numerical experience
confirms this.
The Fourier method can be implemented in such a way that two FFT's
are computed in each step of the ODE solver of Runge-Kutta type, and
these are computed for a vector of the same length as the measured data
vector. On the other hand, the two wavelet methods reduce the problem size,
determined by the noise level, and multiplications by small dense (Meyer)
or banded (Daubechies) matrices are performed in the ODE solver. In this
respect, both wavelet methods, and in particular the Daubechies wavelet
method, have an advantage over the Fourier method, if the data vectors are
very long 6 .
The quarter plane assumption is used mainly because it makes it easier
to obtain stability estimates. The numerical methods can be used also
for the case when the equation is defined only for a bounded interval in
space, provided that also measurements for u x are available. In principle
this requires measurements with two thermocouples. In certain symmetric
configurations it is sufficient to use only one thermocouples.
All three methods considered can easily be applied to problems with
variable coefficients, as was seen in Section 5.1. Numerical methods for
non-linear equations will be the studied in our future research.

Acknowledgement

The numerical experiments in Section 5.2 were based on measurements conducted
at the Department of Mechanical Engineering, Link-oping university.
We are grateful to Ricardo Garcia-Padr'on and Dan Loyd for performing
the experiment and for several discussions on real world heat conduction
problems.
6 In our implementation we have assumed that also intermediate results, for
are of interest. If only the final solution, at needed, then, in principle, the Fourier
method can be implemented with only two FFT's altogether, which could make it faster
than the wavelet methods.


Appendix


A Periodization
As we remarked in Section 4.3, the algorithms for computing the DMT for
Meyer wavelets are based on the Discrete Fourier Transform (DFT) and
use the FFT (Fast Fourier Transform) algorithm. In using the DFT, it is
assumed that the sequence to be transformed is periodic (see, e.g. [3]). In
our application it is not natural to presuppose that the data vectors are
periodic: We have a function u(x; t) (temperature) equal to zero for t ! 0,
and we consider it in the unit square [0; 1] \Theta [0; 1]. Further, we assume that
there is a temperature change at represented by
is diffused through some medium (in the interval 0 - x - 1), and recorded
at 1. Thus we cannot assume that
the functions that we are interested in, are equal to zero for t - 1.
Nevertheless, for numerical reasons, in order to avoid wrap-around effects
(cf. [3]) in the computation of the FFT, it is necessary to assume periodicity.
Therefore we extended the (noisy) data function g m (t) to the interval 1 -
first computing a smoothed version of the last 16 components of
the original data vector, using a cubic smoothing spline. Then the last 7
components of the extended data vector were prescribed to be equal to zero,
and a cubic spline function was constructed that interpolated 16 equidistant
data points from the smooting spline and the 7 zero data points. Finally
this interpolating spline was sampled at equidistant points, as many as in
the original data vector. Thus the size of the data vector was doubled. The
periodization of data vectors is illustrated in Figure 4.1.
Therefore, in effect we solved the sideways heat equation for
but the computed values for t ? 1 are not used. Of course, the computed
values for t less than but close to 1 are affected by the artificial data for t ? 1.
However, these values should not be trusted anyway, see the discussions in
[8, 9, 10].
We used the same periodization for all three methods considered in this
paper.



--R

Inverse Heat Conduction.
On the representation of operators in bases of compactly supported wavelets.
The DFT: An Owner's Manual for the Discrete Fourier Transform.
Determining surface temperatures from interior obser- vations
Space marching difference schemes in the nonlinear inverse heat conduction problem.
Slowly divergent space marching schemes in the inverse heat conduction problem.
Ten Lectures on Wavelets.
The numerical solution of a non-characteristic Cauchy problem for a parabolic equation
Hyperbolic approximations for a Cauchy problem for the heat equation.
Numerical solution of the sideways heat equation.
Numerical solution of the sideways heat equation by difference approximation in time.
Solving the sideways heat equation by a 'method of lines'.
A mollified space marching finite differences algorithm for the inverse heat conduction problem with slab symmetry.
Time Dependent Problems and Difference Methods.
Heat Transfer
Wavelet Methods for the Inversion of Certain Homogeneous Linear Operators in the Presence of Noisy Data.
Continuous data dependence
Multiresolution approximation signal decomposition.

Principe d'incertitude
Ondelettes, fonctions splines et analyses gradu'ees.
The mollification method and the numerical solution of the inverse heat conduction problem by finite differences.
A stable space marching finite differences algorithm for the inverse heat conduction problem with no initial filtering procedure.
Sideways heat equation and wavelets.
Solving the sideways heat equation by a Wavelet-Galerkin method
Prentice Hall PTR
Spectral analysis of the differential operator in wavelet bases.
--TR

--CTR
Fu Chuli , Qiu Chunyu, Wavelet and error estimation of surface heat flux, Journal of Computational and Applied Mathematics, v.150 n.1, p.143-155, January
Chu-Li Fu, Simplified Tikhonov and Fourier regularization methods on a general sideways parabolic equation, Journal of Computational and Applied Mathematics, v.167 n.2, p.449-463, 1 June 2004
Xiang-Tuan Xiong , Chu-Li Fu, Determining surface temperature and heat flux by a wavelet dual least squares method, Journal of Computational and Applied Mathematics, v.201 n.1, p.198-207, April, 2007
Wei Cheng , Chu-Li Fu , Zhi Qian, A modified Tikhonov regularization method for a spherically symmetric three-dimensional inverse heat conduction problem, Mathematics and Computers in Simulation, v.75 n.3-4, p.97-112, July, 2007
cauchy problem;ill-posed;heat conduction;wavelet;inverse problem;fourier analysis