Wavelet-Galerkin Method for Identifying an Unknown Source Term in a Heat Equation

We consider the problem of identification of the unknown source in a heat equation. The problem is ill posed in the sense that the solution if it exists does not depend continuously on the data. Meyer wavelets have the property that their Fourier transform has compact support. Therefore, by expanding the data and the solution in the basis of the Meyer wavelets, high-frequency components can be filtered away. Under the additional assumptions concerning the smoothness of the solution, we discuss the stability and convergence of a wavelet-Galerkin method for the source identification problem. Numerical examples are presented to verify the efficiency and accuracy of the method.


Introduction
Inverse source identification problems are important in many branches of engineering sciences; for example, an accurate estimation of pollutant source is crucial to environmental safeguard in cities with high population.This inverse problem has been investigated in some theoretical papers concerned with the conditional stability and the data compatibility of the solution, notably in 1-6 .The optimal error bound has been given in 7 .Several numerical methods 8-13 have been proposed for the inverse source identification problem.In the present paper, based on some ideas 14-16 , we propose a wavelet-Galerkin method to solve the inverse source problem.
We consider the following initial value problem for the nonhomogeneous heat equation: where u •, t ∈ L 2 R represents state variable and f x denotes the source sink term.

Mathematical Problems in Engineering
The problem 1.1 is a classical direct problem which has been extensively studied in the past decades.Unfortunately, in many practical situations, the characteristics of the source sink term are always unknown.Therefore, the problem is mathematically underdetermined, and an additional data must be supplied to fully determine the physical process.Our task is to determine the heat source on the usual initial conditions with the assistance of additionally supplied data.This is inversely determined and it is usually ill posed in the sense that small perturbation in the data may result, enormous deviation in the solution.
In this paper we consider the inverse problem of determining the function f x in 1.1 from the overspecified condition: This means that our purpose is to determine the pair of functions {u x, t , f x } from the following problem:

1.3
As shown in 17,18 , this problem has a unique solution, but the solution is very sensitive to small data perturbations; hence, it is ill posed.Since g can only be measured in practice, there would be measurement errors and we actually have the noisy data function g m ∈ L 2 R which satisfies where • denotes the L 2 R -norm, and the constant δ > 0 represents the noise level.The ill posedness of problem 1.3 can be seen by solving it in the frequency domain.Let denote the Fourier transform of the function v x .The problem 1.3 can now be formulated in frequency space as follows:

1.6
It is easy to know that the function f ξ in 1.6 can be given by , the existence of a solution depends on a rapid decay of g ξ at high frequencies.However, for the measurement data function g m is merely in L 2 R and in general does not possess such a decay property, high frequency components in the error are magnified and can destroy the solution.Therefore it is impossible to solve the problem using classical numerical methods and some special techniques are required to be employed.Since the convergence rates can only be given under a priori assumptions on the exact solution 19 , we will formulate such an a priori assumption in terms of the exact solution f x by considering where • p denotes the norm in the Sobolev space H p R defined by Meyer's wavelets are special because, unlike most other wavelets, they have compact support in the frequency domain but not in the time domain however, they decay very fast .The wavelet-Galerkin method for approximation solutions of the sideways heat equation has been used in 15, 16 , and so forth.It was demonstrated there that using this method the sideways heat equation can be solved efficiently and in a numerically stable way.
The purpose of this paper is to demonstrate that, using a wavelet-Galerkin approach, we can solve problem 1.3 efficiently.By using the method, we give an error estimate between the exact solution and its approximation, as well as the rule for choosing an appropriate wavelet subspace, depending on the noise level of data.
The outline of the paper is as follows.First in Section 2 we describe Meyer's wavelets and discuss the properties that make them useful for solving ill-posed problems.Then, in Section 3, we describe the wavelet-Galerkin method and give an error estimate which shows the continuous dependence of approximated solution on the data.

The Meyer Wavelets
Let ϕ x , ψ x be Meyer's scaling and wavelet function defined by their Fourier transform in 20 which satisfy π .

2.1
The function family Mathematical Problems in Engineering constitute an orthonormal basis of L 2 R and supp The multiresolution analysis MRA {V j } of the Meyer wavelet is generated by The functions {ψ jk } k∈Z constitute the orthonormal complement W j of V j in V j 1 ; that is, V j 1 V j ⊕ W j .Let P j and Q j denote the orthogonal projections of L 2 R onto V j and W j , respectively.Then the orthogonal projections of any function g ∈ L 2 R on space V j and W j can be expressed by respectively.
It is easy to see from 2.3 and 2.5 that From 2.7 we see that the projection P J can be considered as a low-pass filter: frequencies higher than 4π2 J /3 will be filtered away.

A Galerkin Method in V J
We now introduce the Galerkin method for approximation of the solutions of problem 1.3 based on the separation of variables and using wavelets approximation in the space variable.
Consider the weak formulation of the differential equation with test functions from V J : for all k ∈ Z.The problem 3.1 can be rewritten in the equivalent form: find {u J , f J } ∈ V J such that

3.2
Then with the Ansatz we get the infinite-dimensional system of ordinary differential equations for the vector of coefficients c and μ, where γ { g, ϕ Jν } ν∈Z ; that is, and the matrix D J is given by D J kν ϕ Jν , ϕ Jk .

3.6
The matrix D J is the second-order differentiation operator in V J , and the following boundness guarantees the well-posedness of the Galerkin equation 3.1 .
Proposition 3.1.The infinite matrix D J is symmetric and has Toeplitz structure.Its norm satisfies Moreover, if r is a continuous function, then The proof is similar to 15 and is given in the appendix.
Let us denote

3.9
We are interested in the norm estimation of the distance between the Galerkin solution f δ J of problem 3.1 for the noisy data g m and the unknown solution f of problem 1.3 for the exact data g.Let P J f denote the projection of f on the space V J ; by the triangle inequality we know where the first term of the right-hand side of 3.10 describes the approximation of the exact solution in the scaling subspaces V J , the second one represents the norm of the "error" function and the last one corresponds to the stability of the Galerkin method.Now we consider the three terms,respectively.First, let us consider the problem of stability of the Galerkin solution with respect to perturbations of the data function g.Theorem 3.2.Let g m be the measured data which satisfies 1.4 .Let f J and f δ J be the solution of the Galerkin problem 3.2 with data g and g m , respectively.Then Proof.Let f J and f δ J be given by where μ and μ m are the solution of the Galerkin equations 3.4 with data γ and γ m with the obvious definition of γ .Then, by the Parseval formula,

Mathematical Problems in Engineering 7
Due to Proposition 3.1, we have

3.15
Combining it with 1.4 , we get 3.12 .
Before investigating the relation between the exact solution and the corresponding Galerkin solution with the exact data g, we list two useful lemma and corollary whose proofs are similar to Lemma 3.3 in 16 and will be given in the appendix.Lemma 3.3.If the exact solution f x satisfies a priori condition 1.8 for p > 0, then there holds

3.16
Corollary 3.4.Suppose condition of Lemma 3.3 holds, u x, t is the exact temperature distribution of problem 1.1 and 1.2 , then there holds 3.17

3.18
Proof.Due to 2.7 and 2.8 , we know

3.19
Taking into account the assumption 1.8 , we can estimate the second term of the right-hand side of 3.10 as follows:

3.20
Combining 3.20 with Lemma 3.3, we complete the proof.

Mathematical Problems in Engineering
It remains to reckon with the second term in the right-hand side of 3.10 .
Theorem 3.6.Let f and f J be the solution of problems 1.3 and 3.2 , respectively, for exact data g.If 1.8 is satisfied, then there holds

3.21
Proof.Since the Galerkin equation 3.1 for {u J , f J } satisfies 3.2 , and the pair of functions {P J u, P J f} satisfy P J u t − P J P J u xx − P J I − P J u xx P J f, P J u •, 0 0, 3.23 if we denote v P J u − u J , w P J f − f J , then the error functions {v x, t , w x } satisfy the equation

3.24
Let y • , z • ∈ L 2 0, 1 be the representations in the wavelet basis {ϕ Jk } of the functions v x, t and w x P J I − P J u xx , respectively; that is,

3.25
Then 3.24 is equivalent to the infinite system of first-order differential equations for y {y k } k∈Z and z {z k } k∈Z :

3.26
Then we have 3.28 hence we have in the last inequality we have used Proposition 3.1 and Corollary 3.4.
Theorems 3.2, 3.5, and 3.6 give the estimates of the three terms appearing in the error bound inequality 3.10 ; by combining the three results we can give a H ölder-type error estimate for the wavelet-Galerkin method in the following theorem.Theorem 3.7.Let f be the exact solution of 1.1 and 1.2 satisfying 1.8 for p > 0, and let f δ J be the Galerkin solution of 3.1 for the measured data g m such that 1.4 holds.If J J δ, E satisfies

3.30
then there holds where C is a positive constant independent of E and δ.

Numerical Complement
In this section, we will describe a numerical complement of the proposed method.Note that the problem 1.3 is essentially local.That is a strong source f x at some position x 0 will influence the solution g x for x → x 0 but have limited impact further away.This sort of local property of the problem 1.3 allows us to truncate the problem to a finite internal of x and still obtain reasonable solutions.

Solve u x, 1 g x from the Direct Problem
We select u x, 0 0 and a known f x for 0 ≤ x ≤ 1, and suppose that u 0, t u 1, t 0 this is to agree with compatibility condition and computed data functions u x, t , hence u x, 1 g x by solving a well-posed initial-boundary value problem on the interval 0 ≤ x ≤ 1, using the Crank-Nicolson implicit scheme.It is described as follows: let Δt 1/M and Δx 1/N be the step lengths on time and space coordinates, M, N ∈ N, 0 t 0 < t 1 < • • • < t M 1, and 0 x 0 < x 1 < • • • < x N π denote equidistant partitions of the 0, 1 .We define u j i u x i , t j and f i f x i , and the finite difference approximation is 0, j 0, . . ., M.

4.1
Then we can easily obtain the data u x, t and u x, 1 g x .

Discrete Wavelet Transform
In the numerical solution of 3.4 by an ODE solver, we need to evaluate matrix-vector products D J c.The representation of differentiation operators in bases of compactly supported wavelets is described in the literature; see, for example, 21 .In our context of Meyer's wavelets, which do not have compact support, the situation is different.The proof of Proposition 3.1 actually gives a fast algorithm for this.From the definition of D J , it is easily shown that D J 2 2J D 0 .Thus, we can compute approximations of the elements of D J by first sampling the function Δ equidistantly and then computing its discrete Fourier transform.We will use DMT as a short form of the "discrete Meyer wavelet transform."Algorithms for discretely implementing the Meyer wavelet transform are described in 20 .These algorithms are based on the fast Fourier transform FFT , and computing the DMT of a vector in R requires O n log 2  2 n operations 20 .The algorithms presuppose the vector to be transformed represents a periodic function.So we need to make periodic the vector at first.A discussion on how to make a function "periodic" can be found in 14 .

Solve μ m from Problem 3.4
In the solution of problem 1.3 in V J , we replace the infinite-dimensional ODE 3.4 by the finite-dimensional where c c t ∈ R 2 J represents the approximation of the solution in V J and J is chosen according to Theorem 3.7.For simplicity we suppress the dependence in V J , and since we are dealing with functions, for which only a finite number of coefficients are nonzero, D d J is a finite portion of the infinite matrix D J we use superscript d to indicate this .We also denote Δt 1/M is the step size of time axis t.Define c k c kΔt , k 0, . . ., M.Then, using a modified Euler scheme, we have where By the initial condition,

4.6
We know that if BA −1 < 1, there holds where then we obtain 4.9

Numerical Examples
In this section some numerical examples are presented to demonstrate the usefulness of the approach.The tests were performed using Matlab and the wavelet package WaveLab 850.Suppose that the sequence {g x i } n i 1 represents samples from the function g x on an equidistant grid and n is even, then we add a random uniformly distributed perturbation to each data and obtain the perturbation data, g δ g μ randn size g , 5.1 where

5.2
Then the total noise δ can be measured in the sense of root mean square error according to where "randn • " is a normally distributed random variable with zero mean and unit standard deviation and dictates the level of noise."randn size g " returns an array of random entries that is the same size as g.
The numerical examples were constructed in the following way.First we selected function f x , for 0 ≤ x ≤ 1, and computed u x, t , and hence u x, 1 g x , by solving a well-posed initial-boundary value problem on the domain x, t 0, 1 × 0, 1 , using the Crank-Nicolson implicit scheme see Section 4.1 .Then we added a normally distributed perturbation to data function g giving vectors g δ .From the perturbed data functions, we reconstructed f x and compared the result with the known solution.
Example 5.1.It is easy to verify that the function is the exact solution of problem 1.3 with data g x 2 − e −π 2 sin πx , 0 ≤ x ≤ 1.

5.5
Example 5.2.We examine the reconstruction of a Gaussian normal distribution Mathematical Problems in Engineering 13 where μ 0.5 is the mean and σ 0.1 is the standard deviation.Note that when σ is small expression, 5.6 mimics a Dirac delta distribution δ x − μ .Since the direct problem with f given by 5.6 does not have an analytical solution, the data g is obtained by solving the direct problem using finite difference.
Example 5.3.Consider a continuous piecewise smooth heat source; namely, Example 5.4.This example involves reconstructing a discontinuous heat source given by The results from these examples are given in Figures 1, 2, 3, and 4. In all cases, the length of the data vector g δ was 512.The regularization parameters were selected according to the recipe given in Theorem 3.7.In all cases the number of step length Δt in the ODE solver were 1/20; that is, M 20.Before presenting the results, we recomputed our coarse level approximation on the fine scale, using the inverse Meyer wavelet transform.
Figures 1-4 show that the proposed approach seems to be useful.Moreover, the smaller the error δ, the better the approximation result f δ J .The scheme works equally well for piecewise smooth and discontinuous heat sources.To illustrate this, the numerical results retrieved for Examples 5.3 and 5.4 are presented in Figures 3 and 4. From these figures, it can be seen that the numerical solutions are less accurate than that of Examples 5.1 and 5.2.It is not difficult to see that the well-known Gibbs phenomenon and the recovered data near the nonsmooth and discontinuities points are not accurate.Note that the same situation happened for iterative method 17, 18 .Taking into consideration the ill posedness of the problems, the results presented here are quite satisfactory.

A. Proof of Proposition 3.1
For the proof we use the following two lemmas.
Lemma A.1.The matrix D J is symmetric and has the Toeplitz structure.
Proof.It can be easily shown by integration by parts that D J is symmetric.Moreover, A.1 hence, D J is constant along diagonals; that is, the matrix D J has the Toeplitz structure.Denote D J k the element of the kth diagonal of the matrix D J , then Mathematical Problems in Engineering Lemma A.2.For −π2 J ≤ x ≤ π2 J , define the function extend it periodically, and expand it in the Fourier series Then for all k, δ k d k , where d k is the element in diagonal k of D J .Proof.The Fourier coefficients are given by where we have used the three terms in the definition A.3 of Δ J x .As the result of periodicity, we can rewrite the first term Rewriting δ k similarly, combining the expression for δ k , δ 0 k and δ − k , and noting that ϕ J x 0 for |x| ≥ 4/3 π2 J , we get From the definition of D J , we now see that d k δ k .
We can now prove Proposition 3.1.From 15 we know that, since First, due to Δ J −x Δ J x , Δ J x is an even function, we only need consider the interval 0, π2 J .Here, ϕ J x 2π2 J is identically zero, A.9 Finally we get A.10 Mathematical Problems in Engineering The estimate 3.7 for D J is proved.Since D J is a symmetric matrix, it can be written as

B.2
On the other hand, each coefficient f • , ψ Jk can be written as where Thus,

Theorem 3 . 5 .
If 1.8  is satisfied for a certain p > 0, then t P J I − P J u xx , ϕ Jk x dt;
dE λ , A.11 where E λ is a family of orthogonal projections; see Engl et al. 19 .It follows that if r is a continuous function, B. Proof of Lemma3.3 k∈Z f • , ψ Jk ψ Jk k∈Z f • , ψ Jk e −ikξ•2 −J ψ J .