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

Authors: Fangfang Dou
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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:(1.1)ut-uxx=f(x),x∈R,0<t,u(x,0)=0,x∈R,
where u(·,t)∈L2(ℝ) represents state variable and f(x) denotes the source (sink) term.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 under-determined, 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 functionf(x) in (1.1) from the overspecified condition:(1.2)u(x,1)=g(x),x∈R.
This means that our purpose is to determine the pair of functions {u(x,t),f(x)} from the following problem:(1.3)ut-uxx=f(x),x∈R,0<t,u(x,0)=0,x∈R,u(x,1)=g(x),x∈R.
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 gm∈L2(ℝ) which satisfies(1.4)‖gm(⋅)-g(⋅)‖=‖gm(⋅)-u(⋅,1)‖≤δ,
where ∥·∥ denotes the L2(ℝ)-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(1.5)v̂(ξ)=(Fv)(ξ):=12π∫-∞∞e-iξxv(x)dx,ξ∈R
denote the Fourier transform of the function v(x). The problem (1.3) can now be formulated in frequency space as follows:(1.6)ût(ξ,t)+ξ2û(ξ,t)=f̂(ξ),ξ∈R,0<t,û(ξ,0)=0,ξ∈R,û(ξ,1)=ĝ(ξ),ξ∈R.It is easy to know that the functionf̂(ξ) in (1.6) can be given by(1.7)f̂(ξ)=ξ21-e-ξ2ĝ(ξ).
On account of ξ2/(1-e-ξ2)~𝒪(ξ2) as |ξ|→∞ and f∈L2(ℝ), the existence of a solution depends on a rapid decay of ĝ(ξ) at high frequencies. However, for the measurement data function gm is merely in L2(ℝ) 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(1.8)‖f‖p≤E,p>0,
where ∥·∥p denotes the norm in the Sobolev space Hp(ℝ) defined by(1.9)‖f‖p=(∫-∞∞(1+ξ2)p|f̂(ξ)|2dξ)1/2.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 Section2 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.


## 2. The Meyer Wavelets

Letφ(x), ψ(x) be Meyer’s scaling and wavelet function defined by their Fourier transform in [20] which satisfy(2.1)supp⁡φ̂=[-43π,43π],supp⁡ψ̂=[-83π,-23π]∪[23π,83π].
The function family(2.2)ψjk(x)=2j/2ψ(2jx-k),j,k∈Z
constitute an orthonormal basis of L2(ℝ) and(2.3)supp⁡ψ̂jk=[-83π2j,-23π2j]∪[23π2j,83π2j].
The multiresolution analysis (MRA) {Vj} of the Meyer wavelet is generated by(2.4)Vj={φjk,k∈Z}¯,φjk:=2j/2φ(2jx-k),j,k∈Z,(2.5)supp⁡(φ̂jk)={ξ;|ξ|≤43π2j}.
The functions {ψjk}k∈ℤ constitute the orthonormal complement Wj of Vj in Vj+1; that is, Vj+1=Vj⊕Wj. Let Pj and Qj denote the orthogonal projections of L2(ℝ) onto Vj and Wj, respectively. Then the orthogonal projections of any function g∈L2(ℝ) on space Vj and Wj can be expressed by(2.6)Pjg:=∑k∈Z(f,φjk)φjk,Qjg:=∑k∈Z(f,ψjk)ψjk,
respectively.It is easy to see from (2.3) and (2.5) that(2.7)PJĝ(ξ)=0for|ξ|≥43π2J,(2.8)Qjĝ(ξ)=0forj>J,|ξ|<43π2J.From (2.7) we see that the projection PJ can be considered as a low-pass filter: frequencies higher than 4π2J/3 will be filtered away.


## 3. A Galerkin Method inVJ

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 fromVJ:(3.1)(ut-uxx,φJk)=(f,φJk),(u(⋅,0),φJk)=0,(u(⋅,1),φJk)=(g,φJk),
for all k∈ℤ. The problem (3.1) can be rewritten in the equivalent form: find {uJ,fJ}∈VJ such that(3.2)(uJ)t-PJ(uJ)xx=fJ,uJ(⋅,0)=0,uJ(⋅,1)=PJg.
Then with the Ansatz(3.3)uJ(x,t)=∑ν∈Zc(ν)(t)φJν(x),fJ(x)=∑ν∈Zμ(ν)φJν(x),
we get the infinite-dimensional system of ordinary differential equations for the vector of coefficients c and μ,(3.4)ct-DJc=μ,0≤t≤1,c(0)=0,c(1)=γ,
where γ={(g,φJν)}ν∈ℤ; that is,(3.5)PJg=∑ν∈Zγ(ν)φJν,
and the matrix DJ is given by(3.6)(DJ)kν=(φJν′′,φJk).
The matrix DJ is the second-order differentiation operator in VJ, and the following boundness guarantees the well-posedness of the Galerkin equation (3.1).Proposition 3.1.
The infinite matrixDJ is symmetric and has Toeplitz structure. Its norm satisfies
(3.7)‖DJ‖≤2(π2J)2.
Moreover, if r is a continuous function, then
(3.8)‖r(DJ)‖≤max⁡|λ|≤2(π2J)2|r(λ)|.The proof is similar to [15] and is given in the appendix.Let us denote(3.9)IJ:=[-43π2J,-23π2J]∪[23π2J,43π2J],IJ*:=(-∞,-43π2J]∪[43π2J,∞).We are interested in the norm estimation of the distance between the Galerkin solutionfJδ of problem (3.1) for the noisy data gm and the unknown solution f of problem (1.3) for the exact data g. Let PJf denote the projection of f on the space VJ; by the triangle inequality we know (3.10)‖f-fJδ‖≤‖(I-PJ)f‖+‖PJf-fJ‖+‖fJ-fJδ‖,
where the first term of the right-hand side of (3.10) describes the approximation of the exact solution in the scaling subspaces VJ, the second one represents the norm of the “error” function(3.11)w=PJf-fJ,
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 functiong.Theorem 3.2.
Letgm be the measured data which satisfies (1.4). Let fJ and fJδ be the solution of the Galerkin problem (3.2) with data g and gm, respectively. Then
(3.12)‖fJ-fJδ‖≤(2(π2J)2+1)δ.Proof.
LetfJ and fJδ be given by
(3.13)fJ(x)=∑ν∈Zμ(ν)φJv(x),fJδ(x)=∑ν∈Zμm(ν)φJν(x),
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,
(3.14)‖fJ-fJδ‖=‖μ-μm‖=‖DJeDJ-1(γ-γm)‖≤‖DJeDJ-1‖‖γ-γm‖≤‖DJeDJ-1‖‖g-gm‖.
Due to Proposition 3.1, we have
(3.15)‖DJeDJ-1‖≤max⁡|λ|≤2(π2J)2|λeλ-1|=2(π2J)21-e-2(π2J)2≤2(π2J)2+1.
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 datag, 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 solutionf(x) satisfies a priori condition (1.8) for p>0, then there holds
(3.16)‖QJf(x)‖IJ≤(21-Jπ)1/2E(23π2J)-p.Corollary 3.4.
Suppose condition of Lemma3.3 holds, u(x,t) is the exact temperature distribution of problem (1.1) and (1.2), then there holds
(3.17)‖PJ((I-PJ)uxx)‖≤(21-Jπ)1/2E(23π2J)-p.Theorem 3.5.
If (1.8) is satisfied for a certain p>0, then
(3.18)‖(I-PJ)f‖≤((21-Jπ)1/2+2-p)E(23π2J)-p.Proof.
Due to (2.7) and (2.8), we know
(3.19)‖((I-PJ)f)̂(ξ)‖=‖QJf̂(ξ)‖IJ+‖f̂(ξ)‖IJ*.
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)‖f̂‖IJ*=(∫|ξ|>(4/3)π2J(1+ξ2)-p(1+ξ2)p|f̂|2dξ)1/2≤max⁡|ξ|>(4/3)π2Jξ-pE=(43π2J)-pE.
Combining (3.20) with Lemma 3.3, we complete the proof.It remains to reckon with the second term in the right-hand side of (3.10).Theorem 3.6.
Letf and fJ 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)‖PJf-fJ‖≤(21-Jπ)1/2E(23π2J)-p.Proof.
Since the Galerkin equation (3.1) for {uJ,fJ} satisfies (3.2),
(3.22)(uJ)t-PJ(uJ)xx=fJ,uJ(⋅,0)=0,uJ(⋅,1)=PJg,
and the pair of functions {PJu,PJf} satisfy
(3.23)(PJu)t-PJ(PJu)xx-PJ((I-PJ)u)xx=PJf,PJu(⋅,0)=0,PJu(⋅,1)=PJg;
if we denote v=PJu-uJ, w=PJf-fJ, then the error functions {v(x,t),w(x)} satisfy the equation
(3.24)vt-PJvxx-PJ((I-PJ)u)xx=w,v(⋅,0)=0,v(⋅,1)=0.
Let y(·),z(·)∈L2[0,1] be the representations in the wavelet basis {φJk} of the functions v(x,t) and w(x)+PJ((I-PJ)u)xx, respectively; that is,
(3.25)yk(t)=(v(⋅,t),φJk),zk=(w(x)+PJ((I-PJ)u)xx,φJk).
Then (3.24) is equivalent to the infinite system of first-order differential equations for y={yk}k∈ℤ and z={zk}k∈ℤ:
(3.26)yt-DJy=z,-∞<x<∞,y(0)=0,y(1)=0.
Then we have
(3.27)∫01eDJ(1-t)z(t)dt=0;
taking into account (3.25), we get
(3.28)(w,φJk)=DJeDJt-1∫01eDJt(PJ((I-PJ)u)xx,φJk(x))dt;
hence we have
(3.29)‖PJf-fJ‖=‖DJ1-e-DJ∫01e-DJt(PJ((I-PJ)u)xx,φJk)dt‖≤∫01‖DJeDJteDJt-1‖‖(PJ((I-PJ)u)xx,φJk)‖dt≤(21-Jπ)1/2E(23π2J)-pmax⁡|λ|≤2(π2J)2λeλ-1∫01eλtdt≤(21-Jπ)1/2E(23π2J)-p;
in the last inequality we have used Proposition 3.1 and Corollary 3.4.Theorems3.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.
Letf be the exact solution of (1.1) and (1.2) satisfying (1.8) for p>0, and let fJδ be the Galerkin solution of (3.1) for the measured data gm such that (1.4) holds. If J=J(δ,E) satisfies
(3.30)23π2J=(Eδ)1/(p+2),
then there holds
(3.31)‖f-fJδ‖≤CE2/(p+2)δp/(p+2)+δ=(C+o(1))E2/(p+2)δp/(p+2),forδ⟶0,
where C is a positive constant independent of E and δ.


## 4. 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 x0 will influence the solution g(x) for x→x0 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.


### 4.1. Solveu(x,1)=g(x) from the Direct Problem

We selectu(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∈ℕ, 0=t0<t1<⋯<tM=1, and 0=x0<x1<⋯<xN=π denote equidistant partitions of the [0,1]. We define uij=u(xi,tj) and fi=f(xi), and the finite difference approximation is
(4.1)uij+1-uijΔt-ui-1j+1-2uij+1+ui+1j+1+ui-1j-2uij+ui+1j2Δx2=fi,i=1,…,N-1,j=1,…,M-1,ui0=0,i=0,…,N,u0j=0,uNj=0,j=0,…,M.
Then we can easily obtain the data u(x,t) and u(x,1)=g(x).


### 4.2. Discrete Wavelet Transform

In the numerical solution of (3.4) by an ODE solver, we need to evaluate matrix-vector products DJc. 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 DJ, it is easily shown that DJ=22JD0. Thus, we can compute approximations of the elements of DJ 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 ℝ requires 𝒪(nlog22n) 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].


### 4.3. Solveμm from Problem (3.4)

In the solution of problem (1.3) in VJ, we replace the infinite-dimensional ODE (3.4) by the finite-dimensional(4.2)ct-DJdc=μm,0≤t≤T,c(0)=0,c(T)=γm,
where c=c(t)∈ℝ2J represents the approximation of the solution in VJ and J is chosen according to Theorem 3.7. For simplicity we suppress the dependence in VJ, and since we are dealing with functions, for which only a finite number of coefficients are nonzero, DJd is a finite portion of the infinite matrix DJ (we use superscript d to indicate this). We also denote Δt=1/M is the step size of time axis t. Define ck=c(kΔt), k=0,…,M. Then, using a modified Euler scheme, we have(4.3)ck+1-ckΔt-DJck+1+ck2=μm,c(0)=0,c(M)=γm;
that is,(4.4)ck+1=A-1(Bck+μ),
where(4.5)A:=IΔt-DJ2,B:=IΔt+DJ2.
By the initial condition,(4.6)c1=A-1μm,c2=A-1(Bc1+μm)=A-1(BA-1+I)μm,c3=A-1(Bc2+μm)=A-1[B(A-1(BA-1+I)μm)+μm]=A-1[BA-1(BA-1+I)μm+μm]=A-1∑i=02(BA-1)iμm,⋮γm=cM=A-1∑i=0M-1(BA-1)iμm.
We know that if ∥BA-1∥<1, there holds(4.7)∑i=0M-1(BA-1)i=(I-(BA-1))-1(I-(BA-1)M),
where ∥BA-1∥:=maxi|λi|, {λi}i=1M-1 are the eigenvalues of BA-1. Hence(4.8)γm=cM=A-1(I-BA-1)-1(I-(BA-1)M)μm,
then we obtain(4.9)μm=(I-(BA-1)M)-1(I-BA-1)Aγm.


## 4.1. Solveu(x,1)=g(x) from the Direct Problem

We selectu(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∈ℕ, 0=t0<t1<⋯<tM=1, and 0=x0<x1<⋯<xN=π denote equidistant partitions of the [0,1]. We define uij=u(xi,tj) and fi=f(xi), and the finite difference approximation is
(4.1)uij+1-uijΔt-ui-1j+1-2uij+1+ui+1j+1+ui-1j-2uij+ui+1j2Δx2=fi,i=1,…,N-1,j=1,…,M-1,ui0=0,i=0,…,N,u0j=0,uNj=0,j=0,…,M.
Then we can easily obtain the data u(x,t) and u(x,1)=g(x).


## 4.2. Discrete Wavelet Transform

In the numerical solution of (3.4) by an ODE solver, we need to evaluate matrix-vector products DJc. 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 DJ, it is easily shown that DJ=22JD0. Thus, we can compute approximations of the elements of DJ 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 ℝ requires 𝒪(nlog22n) 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].


## 4.3. Solveμm from Problem (3.4)

In the solution of problem (1.3) in VJ, we replace the infinite-dimensional ODE (3.4) by the finite-dimensional(4.2)ct-DJdc=μm,0≤t≤T,c(0)=0,c(T)=γm,
where c=c(t)∈ℝ2J represents the approximation of the solution in VJ and J is chosen according to Theorem 3.7. For simplicity we suppress the dependence in VJ, and since we are dealing with functions, for which only a finite number of coefficients are nonzero, DJd is a finite portion of the infinite matrix DJ (we use superscript d to indicate this). We also denote Δt=1/M is the step size of time axis t. Define ck=c(kΔt), k=0,…,M. Then, using a modified Euler scheme, we have(4.3)ck+1-ckΔt-DJck+1+ck2=μm,c(0)=0,c(M)=γm;
that is,(4.4)ck+1=A-1(Bck+μ),
where(4.5)A:=IΔt-DJ2,B:=IΔt+DJ2.
By the initial condition,(4.6)c1=A-1μm,c2=A-1(Bc1+μm)=A-1(BA-1+I)μm,c3=A-1(Bc2+μm)=A-1[B(A-1(BA-1+I)μm)+μm]=A-1[BA-1(BA-1+I)μm+μm]=A-1∑i=02(BA-1)iμm,⋮γm=cM=A-1∑i=0M-1(BA-1)iμm.
We know that if ∥BA-1∥<1, there holds(4.7)∑i=0M-1(BA-1)i=(I-(BA-1))-1(I-(BA-1)M),
where ∥BA-1∥:=maxi|λi|, {λi}i=1M-1 are the eigenvalues of BA-1. Hence(4.8)γm=cM=A-1(I-BA-1)-1(I-(BA-1)M)μm,
then we obtain(4.9)μm=(I-(BA-1)M)-1(I-BA-1)Aγm.


## 5. 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(xi)}i=1n 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,(5.1)gδ=g+μrandn(size(g)),
where(5.2)g=(g(x1),…,g(xn))T,xi=(i-1)Δx,Δx=1n-1,i=1,…,n.
Then the total noise δ can be measured in the sense of root mean square error according to(5.3)δ:=‖gδ-g‖l2=1n∑i=1n(giδ-gi)2,
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 functionf(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(5.4)u(x,t)=(2-e-π2t)sin⁡(πx),0≤x≤1,0≤t≤1,f(x)=2π2sin⁡(πx),0≤x≤1,
is the exact solution of problem (1.3) with data
(5.5)g(x)=(2-e-π2)sin⁡(πx),0≤x≤1.Example 5.2.
We examine the reconstruction of a Gaussian normal distribution(5.6)f(x)=1σ2πe-((x-μ)2/2σ2),
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,(5.7)f(x)={6x-2,13≤x≤12,4-6x,12≤x≤23,0,else.Example 5.4.
This example involves reconstructing a discontinuous heat source given by(5.8)f(x)={1,13≤x≤23,0,else.The results from these examples are given in Figures1, 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.Exact solution (solid) and its approximation (dashed) for Example5.1. We select J=4, (a) ϵ=10-2; (b) ϵ=10-3.
(a)(b)Exact solution (solid) and its approximation (dashed) for Example5.2. We select J=3, (a) ϵ=10-2; (b) ϵ=10-3.
(a)(b)Exact solution (solid) and its approximation (dashed) for Example5.3. We select J=4, (a) ϵ=10-3; (b) ϵ=10-4.
(a)(b)Exact solution (solid) and its approximation (dashed) for Example5.4. We select J=4, (a) ϵ=10-3; (b) ϵ=10-4.
(a)(b)Figures1–4 show that the proposed approach seems to be useful. Moreover, the smaller the error δ, the better the approximation result fJδ. 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.



---

*Source: 904183-2012-01-29.xml*