APPROXIMATION METHOD FOR A BACKWARD HEAT CONDUCTION PROBLEM
XIANG-TUAN XIONG, CHU-LI FU, ZHI QIAN, AND XIANG GAO Received 22 March 2005; Revised 11 December 2005; Accepted 12 February 2006
We introduce a central difference method for a backward heat conduction problem (BHCP). Error estimates for this method are provided together with a selection rule for the regularization parameter (the space step length). A numerical experiment is presented in order to illustrate the role of the regularization parameter.
Copyright © 2006 Hindawi Publishing Corporation. All rights reserved.
1. Introduction
The backward heat conduction problem (BHCP) is also referred to as the final bound- ary value problem. In general no solution which satisfies the heat conduction equation with final data and the boundary conditions exists. Even if a solution exists, it will not be continuously dependent on the final data. The BHCP is a typical example of an ill- posed problem which is unstable in numerical simulations and requires special regular- ization methods. In the context of approximation method for this problem, many ap- proaches have been investigated. Such authors as Latt`es and Lions [6], Showalter [10], Ames [1], Miller [9] have approximated the BHCP by quasi-reversibility methods. In [11], Schr¨oter and Tautenhahn established an optimal error estimate for a special BHCP.
Mera and Jourhmane used many numerical methods with regularization techniques to approximate the problem in [8,5], and so forth. A mollification method has been stud- ied by H`ao in [4]. Recently, Liu [7] used a group preserving scheme to solve the backward heat equation numerically. A difference approximation method for solving sideways heat equation was provided by Eld´en in [2], where he used a central difference to replace the time derivative in heat equation. In this paper, we consider a special BHCP [4]:
ut(x,t)=uxx(x,t), x∈R, 0< t <1,
u(x, 1)=ϕ(x), x∈R. (1.1)
We want to obtain the temperature distributionu(x,t) for 0< t <1. Since the dataϕ(·) are based on (physical) observations and are not known with complete accuracy, we assume
Hindawi Publishing Corporation
International Journal of Mathematics and Mathematical Sciences Volume 2006, Article ID 45489, Pages1–9
DOI10.1155/IJMMS/2006/45489
thatϕ(·) andϕδ(·) satisfy
ϕ(·)−ϕδ(·)≤δ, (1.2)
whereϕ(·) andϕδ(·) belong toL2(R),ϕδ(·) denotes the measured data andδ denotes the noise level. As usual, assume that there exists an a priori condition for our problem (1.1):
u(·, 0)≤E, (1.3)
whereEis a positive bound. Problem (1.1) has a unique solution according to [3]. In order to use the Fourier transform technique, we define the Fourier transform of function
f(x)(x∈R) as the following:
f(ξ) := 1
√2π ∞
−∞e−ixξf(x)dx, ξ∈R. (1.4) We consider the problem (1.1) inL2-space with respect to the variablex. Then taking Fourier transform with respect tox, the problem (1.1) can be reformulated in frequency space as follows:
ut(ξ,t)=(iξ)2u(ξ,t), u(ξ, 1) =ϕ(ξ ), ξ∈R. (1.5) The solution to (1.5) is given by
u(ξ,t)=eξ2(1−t)ϕ(ξ). (1.6) Then by inverse Fourier transform, the unique solution of (1.1) can be expressed:
u(x,t)=√1 2π
∞
−∞eixξeξ2(1−t)ϕ(ξ)dξ. (1.7) From (1.6), we can easily see that
u(ξ, 0)=eξ2ϕ(ξ ). (1.8)
Since in our problemu(ξ, t) is assumed to be inL2(R), we see that the exact data function
ϕ(ξ) must decay rapidly asξ→ ∞. In the same time, it is easy to see that a small per- turbation in the dataϕ(ξ) may cause dramatically large error in the solution u(ξ, t). In addition, the magnifying factor iseξ2; hence, it is a severely ill-posed problem. By using central difference with step lengthhto approximate the second derivativeuxx, we can get a regularized problem (with noisy data):
vt(x,t)=v(x+h,t)−2v(x,t) +v(x−h,t)
h2 , x∈R, 0< t <1, v(x, 1)=ϕδ(x), x∈R.
(1.9) Similarly we can easily get a solution in the frequency space for problem (1.9):
v(ξ,t)=e4 sin2(hξ/2)(1−t)/h2ϕδ(ξ). (1.10)
Then by inverse Fourier transform, v(x,t)= 1
√2π ∞
−∞eixξe4 sin2(hξ/2)(1−t)/h2ϕδ(ξ)dξ. (1.11) In the forthcoming section, we will see that (1.11) is a stable approximate solution for the problem (1.1).
2. Error estimates
In this section, we will prove error estimates between the exact solution (1.7) and the approximate solution (1.11). The following conclusion is the main result of this paper.
Theorem 2.1. Supposed thatu(x,t) is given by (1.7) with exact dataϕand thatv(x,t) is given by (1.11) with noisy dataϕδ. If there exists a boundu(·, 0) ≤Eand the data functions satisfyϕ−ϕδ ≤δ, and ifh=2(ln(E/δ))−1/2is chosen, then the following error estimates can be obtained.
(I) IfE/δ < e3/2, then for 0< t <1,
u(·,t)−v(·,t)≤2E1−tδt. (2.1) (II) IfE/δ≥e3/2, then
(II-a) for a fixedtsatisfying 1/(3/2) ln(E/δ)< t <1, u(·,t)−v(·,t)≤max
(1−t)
lnE
δ −1/2
3 2 et
3/2
,E1−tδt
+E1−tδt, (2.2) (II-b) for a fixedtsatisfying 0< t≤(1)/((3/2) ln(E/δ)),
u(·,t)−v(·,t)≤max 1−t
3
lnE δ
E1−tδt,E1−tδt
+E1−tδt. (2.3) Proof. By Parseval relation and (1.8), (1.3), and (1.2), we can get
u(·,t)−v(·,t)
=u(·,t)−v(·,t)=e(1−t)ξ2ϕ(ξ) −e(1−t)(4 sin2(ξh/2)/h2)ϕδ(ξ)
≤e(1−t)ξ2ϕ(ξ) −e(1−t)(4 sin2(ξh/2)/h2)ϕ(ξ)
+e(1−t)(4 sin2(ξh/2)/h2)ϕ(ξ) −e(1−t)(4 sin2(ξh/2)/h2)ϕδ(ξ)
=e−tξ2eξ2ϕ(ξ) −e(1−t)(4 sin2(ξh/2)/h2)e−ξ2eξ2ϕ(ξ) +e(1−t)(4 sin2(ξh/2)/h2)ϕ(ξ) −e(1−t)(4 sin2(ξh/2)/h2)ϕδ(ξ)
≤sup
ξ∈R
e−tξ2−e(1−t)(4 sin2(ξh/2)/h2)−ξ2 u(ξ, 0)
+ sup
ξ∈R
e(1−t)(4 sin2(ξh/2)/h2) δ≤sup
ξ∈RA(ξ)E+ sup
ξ∈RB(ξ)δ,
(2.4)
where
A(ξ) := e−tξ2−e(1−t)(4 sin2(ξh/2)/h2)−ξ2 , B(ξ) :=
e(1−t)(4 sin2(ξh/2)/h2) .
(2.5)
Firstly we can estimateB(ξ)δas
B(ξ)δ≤e(1−t)(4/h2)δ. (2.6) According to the selection ofhinTheorem 2.1, we have
B(ξ)δ≤e(1−t) ln(E/δ)δ=E1−tδt. (2.7) Now we will devote to estimatingA(ξ)E,
A(ξ)E= e−tξ2−e(1−t)(4 sin2(ξh/2)/h2)−ξ2 E
=e−tξ2 1−e(1−t)(4 sin2(ξh/2)/h2−ξ2) E.
(2.8)
Obviously,
4 sin2(ξh/2)
h2 ≤ξ2. (2.9)
Case 1. If|ξ| ≥2/h, hence
e(1−t)(4 sin2(ξh/2)/h2−ξ2)≤1, (2.10) furthermore for the selection ofh, we have
A(ξ)E≤e−tξ2E≤e−t(4/h2)E=E1−tδt. (2.11) Case 2. If|ξ| ≤2/h, from (2.8) and by the inequality 1−e−y≤y(y≥0), we have
A(ξ)≤e−tξ2
ξ2−4 sin2(ξh/2) h2
(1−t). (2.12)
Letr:=h|ξ|/2, and note that 0≤r≤1,ξ2=4r2/h2, we have A(ξ)≤e−t(4r2/h2) 1
h2
4r2−4 sin2r(1−t)
=e−t(4r2/h2) 4
h2(r+ sinr)(r−sinr)(1−t)
≤2e−tr2(4/h2)4
h2(r−sinr)(1−t).
(2.13)
By using the inequality sinr≥r−(r3/6) (r≥0), we get A(ξ)≤2e−tr2(4/h2) 4
h2 r3
6(1−t) :=Q(r). (2.14)
We can easily find that the uniqueness maximum of Q(r) attains at the point r0= (3/2t)(h2/4)=
(3/2t)(lnE/δ)−1, then A(ξ)≤Qmax=Q(r0)=(1−t)
lnE
δ
−1/2 3 2te
3/2
, forr0<1. (2.15) Ifr0>1, then maximum ofQ(r) isQ(1), then
A(ξ)≤Qmax=Q(1)=1−t 3
lnE
δ
E−tδt. (2.16)
Fromr0=
(3/2)(1/ln(E/δ)/t), note that 0< t <1, hence if (3/2)(1/ln(E/δ))>1, that is, lnE/δ <3/2, thenr0>1. We can conclude that the following hold.
(I) If (E/δ)< e3/2, thenr0>1 is always satisfied. Hence for all ξ∈R, combining (2.11) and (2.16) we have
u(·,t)−v(·,t)≤sup
ξ∈RA(ξ)E+ sup
ξ∈RB(ξ)δ
≤maxQ(1)E,E1−tδt+E1−tδt
=max
(1−t) 3
lnE
δ
E1−tδt,E1−tδt
+E1−tδt
≤max
(1−t) 3
3
2E1−tδt,E1−tδt
+E1−tδt
≤
max
(1−t) (2, 1)
+ 1
E1−tδt
≤2E1−tδt−→0, forδ−→0.
(2.17)
(II) IfE/δ≥e3/2, then
(II-a) if 1/(3/2) ln(E/δ)< t <1, thenr0≤1 holds, combining (2.11) and (2.15), we have
u(·,t)−v(·,t)
≤sup
ξ∈RA(ξ)E+ sup
ξ∈RB(ξ)δ≤maxQ(r0)E,E1−tδt+E1−tδt
=max
(1−t) 3
2et 3/2
lnE δ
−1/2
,E1−tδt
+E1−tδt−→0, for fixedt,δ−→0;
(2.18)
(II-b) if 0< t≤1/(3/2) ln(E/δ), thenr0>1 holds, combining (2.11) and (2.16), we have
u(·,t)−v(·,t)
≤sup
ξ∈RA(ξ)E+ sup
ξ∈RB(ξ)δ≤maxQ(1)E,E1−tδt+E1−tδt
=max
(1−t) 3
lnE
δ
E1−tδt,E1−tδt
+E1−tδt−→0, for fixedt δ−→0.
(2.19)
Thus, the proof of the theorem is completed.
It is easy to see that the space step lengthh is the regularization parameter of this problem. In the conclusion, we give a rule for choosing the regularization parameter, which is very important for the study of ill-posed problems.
3. A numerical example
In this section, by a numerical experiment, we will study how the regularization parame- terhinfluences the approximation.
It is easy to verify that the function u(x,t)=√ 1
1 + 4te−x2/(1+4t) (3.1) is the unique solution of the initial problem
ut=uxx, x∈R,t >0,
u(x, 0)=e−x2, x∈R. (3.2)
Hence,u(x,t) given by (3.1) is also the solution of the following backward heat equation for 0≤t <1:
ut=uxx, x∈R, 0≤t <1, u(x, 1)= 1
√5e−x2/5, x∈R. (3.3)
Lets=1−t,w(x,s)=u(x,t), we have
ws= −wxx, x∈R, 0< s <1, (3.4) w(x,s=0)=√1
5e−x2/5, x∈R. (3.5)
Ifwxxatxiis replaced by a second-order central difference, then (3.4) becomes ws(xi,s)= −1
h2
wxi+h,s−2wxi,s+wxi−h,s. (3.6)
First let us list some notations:xi=ih, fori= −n,...,0,...,n.wi=wi(s)=w(ih,s). Then (3.6) with the initial condition (3.5) can be discretized as
⎛
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎝ w−n
... w0
... wn
⎞
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎠
s
=
⎛
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎝ 2 h2 −
1
h2 0
−1 h2
2 h2 −
1 h2 . .. . .. . ..
0 −1
h2 2 h2
⎞
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎠
(2n+1)×(2n+1)
⎛
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎜⎝ w−n
... w0
... wn
⎞
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎟⎠
, (3.7)
⎛
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎝ w−n(0)
... w0(0)
... wn(0)
⎞
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎠
=
⎛
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎜⎜
⎝
√1 5e−x2−n/5
...
√1 5e−x20/5
...
√1 5e−xn2/5
⎞
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎟⎟
⎠
. (3.8)
This ordinary differential equation system can easily be solved. In our numerical imple- mentation we use a Runge-Kutta method. A standard ODE solver ode45 in Matlab is such a method.
Now we focus on the numerical experiment. The main aim is to investigate the role of regularization parameterh. We do the numerical experiment in the intervalx∈[−20, 20]
andt∈[0, 1]. This is reasonable in that the initial data at the pointsx= −20, 20 in (3.5) can be considered to be 0 in the computation by noting that the final valueu(x, 1)→0 in (3.3) whenx→ ±∞. To test the accuracy of an approximate solution, we compute the root mean square error (RMSE) by
E(u)= 1
2n+ 1 n i=−n
(vi−ui)2 (3.9)
at total 2n+ 1 test points at thex-axis, wherevianduiare, respectively, the approximate and exact temperature at a test point. Obviously, (3.9) is to approximateL2norm error.
For convenience, in our numerical experiment the noisy data is generated by
wiδ(0)=wi(0) +δ, i= −n,...,n, (3.10) wherewi(0) are the exact data given in (3.8). Thus, Since (3.9), it followsE(w(0))=δ.
FromFigure 3.1, one can see that for fixedδ=0.001, theL2-error attain a minimum at the “optimal”h=2(ln(E/δ))−1/2≈0.76 (here we takeu(x, 0)L2(R)=E≈1.1).
FromFigure 3.2, we can see that for a fixedh, theL2-error increase astapproaches 0.
0.8 1 0.6 0.4 0.2
0 t
1 0.8 0.6 0.4
h 0 0.05 0.1 0.15 0.2 0.25
L2normerror
Figure 3.1. L2-error for varioustandhwith fixedδ=0.001 (“optimal”h=10/13=. 0.77).
0.8 1 0.6 0.2 0.4
0 t
0.01 0.008 0.006 0.004 0.002 0
δ 0 0.02 0.04 0.06 0.08 0.1 0.12
L2normerror
Figure 3.2. L2-error for varioustandδ(for fixedh=10/13).
Acknowledgments
This work is supported by the National Natural Science Foundation of China (no.
10271050), the Natural Science Foundation of Gansu Province of China (no. 3ZS051- A25-015), and the Fundamental Research Fund for Physics and Mathematics of Lanzhou University (no. Lzu05-05). The authors would like to thank the referees for their valuable comments and suggestions.
References
[1] K. A. Ames, G. W. Clark, J. F. Epperson, and S. F. Oppenheimer, A comparison of regularizations for an ill-posed problem, Mathematics of Computation 67 (1998), no. 224, 1451–1471.
[2] L. Eld´en, Numerical solution of the sideways heat equation by difference approximation in time, Inverse Problems 11 (1995), no. 4, 913–923.
[3] L. C. Evans, Partial Differential Equations, Graduate Studies in Mathematics, vol. 19, American Mathematical Society, Rhode Island, 1998.
[4] D. N. H`ao, A mollification method for ill-posed problems, Numerische Mathematik 68 (1994), no. 4, 469–506.
[5] M. Jourhmane and N. S. Mera, An iterative algorithm for the backward heat conduction problem based on variable relaxtion factors, Inverse Problems in Engineering 10 (2002), no. 4, 293–308.
[6] R. Latt`es and J.-L. Lions, M´ethode de quasi-r´eversibilit´e et applications, Travaux et Recherches Math´ematiques, no. 15, Dunod, Paris, 1967, English translation R. Bellman, Elsevier, New York, 1969.
[7] C. S. Liu, Group preserving scheme for backward heat conduction problems, International Journal of Heat and Mass Transfer 47 (2004), no. 12-13, 2567–2576.
[8] N. S. Mera, L. Elliott, L. D. B. Ingham, and D. Lesnic, An iterative boundary element method for solving the one dimensional backward heat conduction problem, International Journal of Heat and Mass Transfer 44 (2001), no. 10, 1973–1981.
[9] K. Miller, Stabilized quasi-reversibility and other nearly-best-possible methods for non-well-posed problems, Symposium on Non-Well-Posed Problems and Logarithmic Convexity (Heriot-Watt University, Edinburgh, 1972), Lecture Notes in Mathematics, vol. 316, Springer, Berlin, 1973, pp. 161–176.
[10] R. E. Showalter, The final value problem for evolution equations, Journal of Mathematical Analysis and Applications 47 (1974), 563–572.
[11] U. Tautenhahn and T. Schr¨oter, On optimal regularization methods for the backward heat equa- tion, Zeitschrift f¨ur Analysis und ihre Anwendungen 15 (1996), no. 2, 475–493.
Xiang-Tuan Xiong: Department of Mathematics, Lanzhou University, Lanzhou 730000, China
E-mail address:[email protected]
Chu-Li Fu: Department of Mathematics, Lanzhou University, Lanzhou 730000, China E-mail address:[email protected]
Zhi Qian: Department of Mathematics, Lanzhou University, Lanzhou 730000, China E-mail address:[email protected]
Xiang Gao: Department of Mathematics, Lanzhou University, Lanzhou 730000, China E-mail address:[email protected]
Special Issue on
Modeling Experimental Nonlinear Dynamics and Chaotic Scenarios
Call for Papers
Thinking about nonlinearity in engineering areas, up to the 70s, was focused on intentionally built nonlinear parts in order to improve the operational characteristics of a device or system. Keying, saturation, hysteretic phenomena, and dead zones were added to existing devices increasing their behavior diversity and precision. In this context, an intrinsic nonlinearity was treated just as a linear approximation, around equilibrium points.
Inspired on the rediscovering of the richness of nonlinear and chaotic phenomena, engineers started using analytical tools from “Qualitative Theory of Differential Equations,”
allowing more precise analysis and synthesis, in order to produce new vital products and services. Bifurcation theory, dynamical systems and chaos started to be part of the mandatory set of tools for design engineers.
This proposed special edition of the Mathematical Prob- lems in Engineering aims to provide a picture of the impor- tance of the bifurcation theory, relating it with nonlinear and chaotic dynamics for natural and engineered systems.
Ideas of how this dynamics can be captured through precisely tailored real and numerical experiments and understanding by the combination of specific tools that associate dynamical system theory and geometric tools in a very clever, sophis- ticated, and at the same time simple and unique analytical environment are the subject of this issue, allowing new methods to design high-precision devices and equipment.
Authors should follow the Mathematical Problems in Engineering manuscript format described at http://www .hindawi.com/journals/mpe/. Prospective authors should submit an electronic copy of their complete manuscript through the journal Manuscript Tracking System athttp://
mts.hindawi.com/according to the following timetable:
Manuscript Due December 1, 2008 First Round of Reviews March 1, 2009 Publication Date June 1, 2009
Guest Editors
José Roberto Castilho Piqueira,Telecommunication and Control Engineering Department, Polytechnic School, The University of São Paulo, 05508-970 São Paulo, Brazil;
Elbert E. Neher Macau,Laboratório Associado de Matemática Aplicada e Computação (LAC), Instituto Nacional de Pesquisas Espaciais (INPE), São Josè dos Campos, 12227-010 São Paulo, Brazil ; [email protected] Celso Grebogi,Center for Applied Dynamics Research, King’s College, University of Aberdeen, Aberdeen AB24 3UE, UK; [email protected]
Hindawi Publishing Corporation http://www.hindawi.com