Fourth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 03, 1999, pp 119–125. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

Determination of the source/sink term in a heat equation ^{∗}

Ping Wang & Kewang Zheng

Abstract

In this work, we consider the problem of determining an unknown parameter in a heat equation with ill-posed nature. Applying Tikhonov regularization, we obtain a stable approximation to the unknown parameter from over-specified data. We also present numerical computations that verify the accuracy of our approximation.

1 Introduction

Cannon and Zachmann [4] considered the question of determining an unknown source in the heat equation from over-specified data. More precisely, find the source f(t) in the heat equation

u_{t}(x, t) = u_{xx}(x, t) + f(t), 0 < x, 0 < t < T ,
u(x, 0) = g(x), 0 < x ,
u(0, t) = ϕ(t), 0 < t < T ,  (1)
u_{x}(0, t) = ψ(t), 0 < t < T ,
g(0) = ϕ(0) = 0 ,

were u(x, t) is the unkown temperature, and ϕ(t), ψ(t), g(x) are the known data. Assuming that ϕ and ψ are smooth functions, Cannon and Zachmann were able to determine the source or the sink term f(t) explicitly or implicitly in several cases. In this short note, we will study (1) with non-smooth data applying the regularization approach used in [9].

Assume that the pair (u, f) is a classical solution of (1). Then

u(x, t) = ∫_{0}^{t} f(τ)d τ − 2 ∫_{0}^{t} K(x, t − τ)ψ(τ)d τ  (2)

+ ∫_{0}^{∞} g(ξ)(K(x − ξ, t) + K(x + ξ, t))d ξ,

∗Mathematics Subject Classifications: 35K05, 35R25.

Key words: Heat equation, inverse problem, regularization.

\begin{matrix}c \\ ◯ \end{matrix}2000 Southwest Texas State University and University of North Texas.

Published July 10, 2000.

where K(x, t) = \frac{1}{\sqrt{4πt}}e^{−x^{2}/(4t)} is the heat kernel.

Therefore, f(t) can be expressed as the solutions to the integral equation

Af = F ,  (3)

where

Af(t) = ∫_{0}^{t} f(τ) dτ ,  (4)

F (t) = ϕ(t) + \frac{1}{\sqrt{π}} ∫_{0}^{t} \frac{ψ(τ)}{\sqrt{t − τ}} dτ − \frac{1}{\sqrt{πt}} ∫_{0}^{∞} g(x)e^{− \frac{x^{2}}{4t}} dx .

Now we see that problem (1) is equivalent to (3). Therefore, we will focus our attention on this equation (3).

2 Ill-posedness and Regularization

For practical purposes, it is more interesting to assume that the data functions are non-smooth. Suppose that ϕ, ψ ∈ L^{2}[0, T] and g ∈ L^{2}[0, ∞] with compact support in [0, T]. Now the integral operator A defined in (3) from space C[0, T] to space L^{2}[0, T] is not surjective and the inverse operator A^{−1} defined on the range of A is not continuous. This means that the problem of solving equation (3) for f in C from data F ∈ L^{2} is ill-posed in the sense of Hadamard [7]

In what follows, we will apply a regularization technique to construct a regularizing operator for equation (3) and then define an approximation to the unknown term f. For this, we first introduce the Tikhonov functional

M^{α}[f, F] = ∥Af − F ∥_{L^{2}[0,T]}^{2} + α∥f∥_{W_{2}^{1}[0,T]}^{2},  (5)

where α is a positive parameter.

Theorem 2.1 For every function F ∈ L^{2}[0, T] and every positive number α, there exists a unique function f_{α} ∈ W_{2}^{1}[0, T] that minimizes the functional M^{α}.

Proof: Consider the first variation of the functional M^{α}. A straightforward calculation shows that the minimizer f_{α} is the solution of the following Euler differential integral equation

α[f′′(t) − f(t)] = ∫_{t}^{T} (∫_{0}^{τ} f(ξ)dξ − F (τ))dτ,  (6)

subject to boundary conditions f′(0) = f′(T) = 0. It is also easy to show that the solution of (6) in W_{2}^{1} is unique. ♢

Now for each α > 0, each F ∈ L^{2}[0, T], we define the operator f_{α} = R(F, α). For an approximate data function F_{δ} (δ measures the error in data), it is important to choose an appropriate parameter α(δ) so that the according minimizer f_{α}(δ) = R[F_{δ}, α(δ)] can be taken as a stable approximate solution of (3). The following theorem shows how to choose the regularizing parameter α.

Theorem 2.2 Let f_{T} ∈ C^{1}[0, T] be the exact solution corresponding to the exact data F_{T}, and F_{δ} be approximate datum. Then for every positive ϵ there exists δ(ϵ) such that the inequality

∥F_{δ} − F_{T} ∥_{L^{2}[0,T]} ≤ δ ≤ δ(ϵ)

implies

∥f_{α(δ)} − f_{T} ∥_{C[0,T]} < ϵ,  (7)

where f_{α(δ)} = R(F_{δ}, α(δ)) with α = α(δ) = δ^{λ} and 0 < λ ≤ 2.

Proof: By the definition of f_{α(δ)}, we know

M^{α(δ)}[f_{α(δ)}, F_{δ}] ≤ M^{α(δ)}[f_{T} , F_{δ}].

That is

∥Af_{α(δ)} − F_{δ}∥_{L^{2}}^{2} + α(δ)∥f_{α(δ)}∥_{W_{2}^{1}}^{2} ≤ ∥Af_{T} − F_{δ}∥_{L^{2}}^{2} + α(δ)∥f_{T} ∥_{W_{2}^{1}}^{2}
≤ δ^{2} + δ^{λ}∥f_{T} ∥_{W_{2}^{1}}^{2}
≤ δ^{λ}d^{2}, d = (1 + ∥f_{T} ∥_{W_{2}^{1}}^{2})^{1/2}.

Hence, ∥f_{α(δ)}∥_{W_{2}^{1}} ≤ d and ∥Af_{α(δ)} − F_{δ}∥_{L^{2}} ≤ dδ^{λ/2}. It is easy to see that both f_{α(δ)} and f_{T} belong to the set E = {f : ∥f∥_{W_{2}^{1}[0,T]} ≤ d}, which is a compact subset of space C[0, T]. The continuity of A^{−1} on A E implies that

∥f_{α(δ)} − f_{T} ∥_{C[0,T]} ≤ ∥A^{−1}∥ · ∥Af_{α(δ)} − Af_{T} ∥_{L^{2}}
≤ ∥A^{−1}∥(∥Af_{α(δ)} − F_{δ}∥_{L^{2}} + ∥Af_{T} − F_{δ}∥_{L^{2}})
≤ ∥A^{−1}∥(dδ^{λ/2} + δ)
≤ δ^{λ/2}∥A^{−1}∥(1 + d) .

By setting

δ(ϵ) = [\frac{ϵ}{∥A^{−1}∥(1 + d)}]^{2/λ}

we obtain (7) and the proof is complete. ♢

Next, we show that F depends continuously on the initial data ϕ, ψ, g.

Theorem 2.3 Suppose that exact data F_{T} , ϕ_{T} , ψ_{T} , g_{T} satisfy (4), and that the appximate data F_{δ}, ϕ_{δ}, ψ_{δ}, g_{δ} also satisfy (4). Then inequalities ∥ϕ_{T} − ϕ_{δ}∥_{L^{2}} ≤ δ, ∥ψ_{T} − ψ_{δ}∥_{L^{2}} ≤ δ and ∥g_{T} − g_{δ}∥_{L^{2}} ≤ δ imply

∥F_{T} − F_{δ}∥_{L^{2}} ≤ Dδ, D = [6(1 + \frac{2T}{π} + \sqrt{\frac{T}{2π}})]^{1/2}

Proof: Applying Cauchy’s inequality, we have

∥F_{δ} − F_{T} ∥_{L^{2}}^{2} ≤ 3(∫_{0}^{T} [ϕ_{T} − ϕ_{δ}]^{2}dt + \frac{1}{π} ∫_{0}^{T} [∫_{0}^{t} \frac{ψ_{δ}(τ) − ψ_{T} (τ)}{\sqrt{t − τ}}]^{2}dτ

+\frac{1}{π} ∫_{0}^{T} \frac{1}{t}[∫_{0}^{∞}(g_{δ}(x) − g_{T} (x))e^{−x^{2}/(4t)}dx]^{2}dt)

≤ 3(δ^{2} + \frac{2}{π} ∫_{0}^{T} [ϕ_{δ}(τ) − ϕ_{T} (τ)]^{2} ∫_{τ}^{T} \sqrt{\frac{t}{t − τ}} dt dτ

+\frac{δ^{2}}{π} ∫_{0}^{T} \frac{1}{t} ∫_{0}^{∞} e^{−x^{2}/(2t)} dx dt)

≤ 3δ^{2}(1 + \frac{4T}{π} + \sqrt{\frac{2T}{π}} )
< D^{2}δ^{2}.
♢

Combining Theorems 2.2, 2.3, we obtain the following stability theorem.

Theorem 2.4 Suppose f_{T} is the exact solution of (3) corresponding to data functions ϕ_{T} , ψ_{T} , g_{T}. For any ϵ > 0 and approximate data ϕ_{δ}, ψ_{δ}, g_{δ}, there exists a δ(ϵ) and an α(δ)such that inequalities ∥ϕ_{T} − ϕ_{δ}∥_{L^{2}} ≤ δ, ∥ψ_{T} − ψ_{δ}∥_{L^{2}} ≤ δ and ∥g_{T} − g_{δ}∥_{L^{2}} ≤ δ imply that

∥f_{α(δ)} − f_{T} ∥_{C[0,T]} < ϵ,  (8)

where f_{α(δ)} = R(F_{δ}, α(δ)).

The above result shows that, for carefully chosen α, f_{α(δ)}, the minimizer of functional (5), can be taken as a stable approximate solution of (1).

3 Numerical Verification

We will study a concrete overdetermined system in this section to numerically test the applicability of the regularization approach discussed in Section 2.

For T = 1, we take

ϕ_{T} (t) = \frac{t^{3}}{3} − \frac{t^{4}}{2} + 0.0002\sqrt{\frac{t}{π}} (1 + 4t(e^{− \frac{1}{4t}} − 1)),
ψ_{T} (t) = \frac{256t^{4.5}}{315\sqrt{π}} ,
g_{T} (x) = 0.0001x(1 − x^{2}).

Then F_{T} (t) = \frac{t^{3}}{3} − \frac{t^{4}}{2} + \frac{t^{5}}{5} . The corresponding exact solution of (3) is

f_{T} (t) = t^{2}(t − 1)^{2}.

It remains to be seen how well the equation (6) recovers the value of f_{T} with the following altered initial data

ϕ_{δ}(t) = ϕ_{T} (t) + δ sin(50πt),
ψ_{δ}(t) = ψ_{T} (t) + δ sin(50πt),
g_{δ}(x) = g_{T} (x) + δ(1 − x)(1 − (1 − x)^{2}).

First of all, (6) is replaced by its finite difference approximation on a uniform grid with step h = T/(n + 1). Thus we obtain the following system of linear equations in which the coefficient matrix is of five diagonal form:

A^{h}f^{h} = h^{3}DF^{h},

where

A^{h} = ( \begin{matrix}ã & b̃ & α \\ b̃ & a & b & α \\ α & b & a & b & α \\ & · & · & · & · & · \\ & & α & b & a & b & α \\ & & & α & b & a & b̃ \\ & & & & α & b & ˜ã \end{matrix} ) , D = ( \begin{matrix}1 \\ −1 & 1 \\ & −1 & 1 \\ & & & −1 & 1 \\ & & & & −1 & 1 \end{matrix} ) ,

a = h^{4} + 2α(h^{2} + 3), ã = h^{4} + α(h^{2} + 2), ˜ ã = h^{4} + α(2h^{2} + 3), b = −α(h^{2} + 4), b̃ = −α(h^{2} + 3). F^{h} = (F_{1}^{h}, ⋯, F_{n}^{h}), f^{h} = (f_{1}^{h}, ⋯, f_{n}^{h}) are difference functions (f_{0}^{h} = f_{1}^{h}, f_{n+1}^{h} = f_{n}^{h}).

The difference scheme for (4) is

F_{i}^{h} = ϕ_{i}^{h} + \frac{1}{\sqrt{π}} ∑_{j=1}^{i} b_{ij}ψ_{j}^{h} − ∑_{j=1}^{n} c_{ij}g_{i}^{h}, i = 1, 2, ⋯, n .

where

b_{ij} = \{ \begin{matrix}2\sqrt{h}(\sqrt{i − j + 1} − \sqrt{i − j}) & j ≤ i \\ 0 & j > i \end{matrix} ,

c_{ij} = erf(\frac{j + 1}{2}\sqrt{\frac{h}{i}}) − erf(\frac{j}{2}\sqrt{\frac{h}{i}}), i, j = 1, 2, …, n,

and

ϕ_{i}^{h} = ϕ(ih), ψ_{j}^{h} = ψ(jh), g_{j}^{h} = g(jh).

The results of the numerical simulation are shown in the following table (δ = 0.0001 and λ = 1.2. f_{α1}, f_{α2}, f_{α3} are approximate solutions corresponding to n = 39, 79, 159 respectively).

One can see from the data in the Table 3 that the numbers generated through the computation show that the approximate solutions and the exact solution match better as n becomes larger. The numbers also show that the approximation when t is very close to 0 is not nearly as good as the approximation

\begin{matrix}t & f_{T} (t) & f_{α1}(t) & f_{α2}(t) & f_{α3}(t) \\ 0.025 & 0.000594 & 0.001808 & 0.003578 & 0.004813 \\ 0.05 & 0.002256 & 0.003663 & 0.005205 & 0.006199 \\ 0.1 & 0.008100 & 0.011018 & 0.011220 & 0.011334 \\ 0.2 & 0.025600 & 0.031311 & 0.028778 & 0.027295 \\ 0.3 & 0.044100 & 0.049948 & 0.046540 & 0.044596 \\ 0.4 & 0.057600 & 0.062285 & 0.059168 & 0.057404 \\ 0.5 & 0.062500 & 0.066572 & 0.063701 & 0.062082 \end{matrix}

Table 1: Exact and approximate solutions

elsewhere. We think it is because we know little about f(0) in advance. The only assumption on f at 0 is f′(0) = 0 (see (6)). Therefore we do not have much control of f when t is very very small. Overall, the table shows that, for large n, our regularization approach i s a reliable way of recovering unknown source or sink term in a heat equation from non-smooth overspecified data.

References

[1] Cannon, J. R., Determination of an unknown coefficient in a parabolic differential equation, Duke Math. J., 30, pp 313-323, (1963).

[2] Cannon, J. R., Determination of the unknown coefficient k(u) in the equation γ ·k(u)γu = 0 from overspecified boundary data, J. Math. Anal. Appl,, 18, (1967).

[3] Cannon, J. R. and Lin, Y.,. An inverse problem of finding a parameter in a semi-linear heat equation, J. Math. Anal. Appl., Vol. 145, No. 2, 1990.

[4] Cannon, J. R. and Zachmann, D., Parameter determination in parabolic partial differential equations from overspecified boundary data, Int. J. Engineering, 20 (6), pp. 779-788, (1982).

[5] Rundell, W., An inverse problem for a parabolic partial differential equation, Rocky Mountain J. Math., 13, pp. 667-688, (1983).

[6] Shcheblov, A. Yu., A method for the approximate solution of an inverse problem for the heat-conduction equation, Comp. Math. and Math. Phys., Vol. 32, No. 6, 1992.

[7] Tikhonov, A. N. and Arsenin, V. Y., Solution of Ill-posed Problems, John Wiley & Sons, (1977).

[8] Wang, P. and Zheng, K., Regularization of an Abel equation, Integr. Equ. Oper. Theory., 29, pp. 243 - 249, (1997).

[9] Zheng, K. and Wang, P., Numerical Approximation of an unknown boundary term in a heat equation, Neur. Parall. Sci. Comp., 2, pp. 451 - 458, (1994).

[10] Wang, P., Zheng, K., Reconstruction of heat sources in heat conduction equations, to appear in Comput. and Appl. Math.

[11] Wang, P., Zheng, K., Determination of conductivity in a heat equation, to appear in Int. J. of Math. and Math. Sci.

Ping Wang Department of Mathematics, Pennsylvania State University Schuylkill Haven, PA 17972, USA e-mail: pxw10@psu.edu

Kewang Zheng Department of Mathematics Hebei University of Science and Technology, China