Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations

Yikan Liu
Research Institute for Electronic Science, Hokkaido University

Abstract

This article is concerned with the derivation of numerical reconstruction schemes for the inverse moving source problem on determining source profiles in (time-fractional) evolution equations. As a continuation of the theoretical result on the uniqueness, we adopt a minimization procedure with regularization to construct iterative thresholding schemes for the reduced backward problems on recovering one or two unknown initial value(s). Moreover, an elliptic approach is proposed to solve a convection equation in the case of two profiles.

Keywords

Inverse moving source problem, (Time-fractional) evolution equation, Iterative thresholding scheme, Convection equation

\@hangfrom\@seccntformat

sectionIntroduction

\@xsect

2.3ex plus.2ex

Let $0<\alpha\leq 2$ , $T>0$ be constants and $\Omega\subset\mathbb{R}^{d}$ ( $d=1,2,\ldots$ ) be a bounded domain with a smooth boundary $\partial\Omega$ . In this article, we consider the following initial-boundary value problem for a (time-fractional) evolution equation

\begin{cases}(\partial_{0+}^{\alpha}-\triangle)u=F&\mbox{in }\Omega\times(0,T),\\ \partial_{t}^{k}u=0\ (k=0,\lceil\alpha\rceil-1)&\mbox{in }\Omega\times\{0\},\\ u=0&\mbox{on }\partial\Omega\times(0,T),\end{cases}

(1.1)

where the source term $F$ takes the form

F(\bm{x},t):=\begin{cases}f(\bm{x}-\bm{p}t),&0<\alpha\leq 1,\\ f(\bm{x}-\bm{p}t)+g(\bm{x}-\bm{q}t),&1<\alpha\leq 2.\end{cases}

(1.2)

In (1.1), by $\triangle:=\sum_{j=1}^{d}\frac{\partial^{2}}{\partial x_{j}^{2}}$ we denote the usual Laplacian in space, and $\lceil\,\cdot\,\rceil$ is the ceiling function. The notation $\partial_{0+}^{\alpha}$ stands for the forward Caputo derivative in time of order $\alpha$ , which will be defined in Section Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations. In (1.2), we assume that $\bm{p},\bm{q}\in\mathbb{R}^{d}$ are constant vectors, and assumptions on the space-dependent functions $f,g$ will be specified also in Section Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations.

As a continuation of the theoretical counterpart studied in [17], this article is concerned with the establishment of numerical reconstruction schemes for the following inverse moving source problem regarding (1.1)–(1.2).

Problem 1.1 (Inverse moving source problem)

Let $u$ be the solution to (1.1)–(1.2), and $\omega\subset\Omega$ be a suitably chosen nonempty subdomain of $\Omega$ . Given constant vectors $\bm{p},\bm{q}\in\mathbb{R}^{d}$ such that $\bm{p}\neq\bm{q},$ determine one source profile $f$ in the case of $0<\alpha\leq 1$ or two source profiles $f,g$ in the case of $1<\alpha\leq 2$ by the partial interior observation of $u$ in $\omega\times(0,T)$ .

Except for the representative evolution equations of parabolic and hyperbolic types, the governing equation in (1.1) is called a time-fractional diffusion equation if $0<\alpha<1$ and a time-fractional diffusion-wave one if $1<\alpha<2$ . In recent decades, time-fractional evolution equations have gathered increasing attention among applied mathematicians and, especially, general linear theories for (1.1) with $0<\alpha<1$ have been mostly established within the last decade (see, e.g., [24, 4, 13]). Meanwhile, the corresponding numerical methods and inverse problems have also been studied intensively, and we refer e.g. to [16, 11] and [12, 20, 14, 19] with the references therein, respectively. Compared with the case of $0<\alpha<1$ , initial-boundary value problems such as (1.1) with $1<\alpha<2$ have not been well investigated from both theoretical and numerical aspects. On this direction, we refer to [24] for the well-posedness of forward problems, and [15] as a recent progress on numerical approaches to related inverse problems.

Recently, it has been recognized that (time-fractional) evolution equations with orders $0<\alpha<2$ actually share several common properties, which can be applied to the qualitative analysis of some inverse problems. Besides the time-analyticity asserted in [24], the weak vanishing property was established in [7, 17] for $0<\alpha<1$ and $1<\alpha<2$ respectively, which resembles the classical unique continuation principle for parabolic equations. As direct applications, this property was employed to prove the uniqueness of an inverse $\bm{x}$ -source problem in [7] and that of Problem 1.1 in [17]. On the other hand, since hyperbolic equations exhibit essentially different properties from the cases of $0<\alpha<2$ , the uniqueness for Problem 1.1 with the exceptional case of $\alpha=2$ requires special treatments and extra assumptions. The interested readers are referred to [17] for theoretical details of Problem 1.1, whose main result will be recalled in Section Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations.

By a moving source we refer to an inhomogeneous term in (1.1) basically taking the form of $h(\bm{x}-\bm{\rho}(t))$ , where $h:\Omega\longrightarrow\mathbb{R}$ and $\bm{\rho}:[0,T]\longrightarrow\mathbb{R}^{d}$ stand for the profile and the orbit of a moving object, respectively. In contrast to conventional inverse source problems, literature on inverse moving source problems seems rather limited, among which the majority focused on hyperbolic equations due to practical significance. We refer e.g. to [22, 23] for the reconstruction of moving point sources in wave equations, and [5] for the unique determination of profiles or orbits in Maxwell equations. Very recently, the stability for identifying orbits was obtained in [6] for (time-fractional) evolution equations such as (1.1) with $0<\alpha\leq 2$ .

Parallel to [6], the uniqueness for Problem 1.1 was newly verified in [17]. Aiming at possible real-world applications, this article is concerned with the development of efficient reconstruction schemes for Problem 1.1, which consists of three key ingredients. The first one coincides with that of [17], that is, reducing the original problem to a backward problem in which the unknown profiles appear as initial conditions. The second one follows the same line as that in [18, 9, 10, 7], which utilized the iterative thresholding algorithm (see [2]) to deal with the corresponding minimization problem. Finally, for the simultaneous recovery of two profiles when $1<\alpha\leq 2$ , the last step is to solve a convection equation, which is reduced to solving a series of second order ordinary differential equations.

The remainder of this article is organized as follows. In Section Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations, we prepare necessary tools of fractional calculus and recall theoretical facts regarding (1.1)–(1.2) and Problem 1.1. The derivations of reconstruction schemes for one and two unknown profiles are explained in Sections Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations and Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations, respectively, and conclusions with future topics are given in Section Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations.

\@hangfrom\@seccntformat

sectionPreliminaries and revisit of theoretical results

\@xsect

2.3ex plus.2ex

We start with fixing notations and terminologies. Throughout this article, all vectors are by default column vectors, e.g., $\bm{x}=(x_{1},\ldots,x_{d})^{\mathrm{T}}\in\mathbb{R}^{d}$ , where $(\,\cdot\,)^{\mathrm{T}}$ denotes the transpose. The inner product in $\mathbb{R}^{d}$ is denoted by $\bm{x}\cdot\bm{y}$ , and the Euclidean distance $|\cdot|$ is induced as $|\bm{x}|=(\bm{x}\cdot\bm{x})^{1/2}$ . By $L^{2}(\Omega)$ we denote the space of square integrable functions in $\Omega$ , and let $H^{k}(\Omega)$ , $H_{0}^{k}(\Omega)$ , etc. and $W^{k,\gamma}(0,T;H_{0}^{1}(\Omega))$ , etc. ( $k=1,2,\ldots$ , $\gamma\in[1,\infty]$ ) be (vector-valued) Sobolev spaces (see e.g. [1, 3]).

For the definition of $\partial_{0+}^{\alpha}$ in (1.1), we recall the forward Riemann-Liouville integral operator for $\beta\in[0,1]$ :

J_{0+}^{\beta}h(t):=\left\{\!\begin{aligned} &h(t),&\quad&\beta=0,\\ &\frac{1}{\Gamma(\beta)}\int_{0}^{t}\frac{h(\tau)}{(t-\tau)^{1-\beta}}\,\mathrm{d}\tau,&\quad&0<\beta\leq 1,\end{aligned}\right.\quad h\in C[0,\infty),

where $\Gamma(\,\cdot\,)$ is the Gamma function. Then for $\beta>0$ , the forward Caputo derivative $\partial_{0+}^{\beta}$ and the forward Riemann-Liouville derivative $D_{0+}^{\beta}$ are formally defined as

\partial_{0+}^{\beta}=J_{0+}^{\lceil\beta\rceil-\beta}\circ\frac{\mathrm{d}^{\lceil\beta\rceil}}{\mathrm{d}t^{\lceil\beta\rceil}},\quad D_{0+}^{\beta}=\frac{\mathrm{d}^{\lceil\beta\rceil}}{\mathrm{d}t^{\lceil\beta\rceil}}\circ J_{0+}^{\lceil\beta\rceil-\beta},

where $\circ$ denotes the composition. Meanwhile, we also introduce the backward Riemann-Liouville integral operator as

J_{T-}^{\beta}h(t):=\left\{\!\begin{aligned} &h(t),&\quad&\beta=0,\\ &\frac{1}{\Gamma(\beta)}\int_{t}^{T}\frac{h(\tau)}{(\tau-t)^{1-\beta}}\,\mathrm{d}\tau,&\quad&0<\beta\leq 1,\end{aligned}\right.\quad h\in C(-\infty,T],

and we further define the backward Caputo and Riemann-Liouville derivatives for $\beta>0$ as

\partial_{T-}^{\beta}=J_{T-}^{\lceil\beta\rceil-\beta}\circ\frac{\mathrm{d}^{\lceil\beta\rceil}}{\mathrm{d}t^{\lceil\beta\rceil}},\quad D_{T-}^{\beta}=\frac{\mathrm{d}^{\lceil\beta\rceil}}{\mathrm{d}t^{\lceil\beta\rceil}}\circ J_{T-}^{\lceil\beta\rceil-\beta}.

For later use, we derive the useful fractional version of integration by parts.

Lemma 2.1

Let $h_{1},h_{2}\in C^{\lceil\alpha\rceil}[0,T]$ . If $0<\alpha\leq 1,$ then

\int_{0}^{T}(\partial_{0+}^{\alpha}h_{1})\,h_{2}\,\mathrm{d}t+h_{1}(0)(J_{T-}^{1-\alpha}h_{2})(0)=-\int_{0}^{T}h_{1}\,(D_{T-}^{\alpha}h_{2})\,\mathrm{d}t+h_{1}(T)(J_{T-}^{1-\alpha}h_{2})(T).

(2.1)

If $1<\alpha\leq 2,$ then

	$\displaystyle\int_{0}^{T}(\partial_{0+}^{\alpha}h_{1})\,h_{2}\,\mathrm{d}t+h_{1}^{\prime}(0)(J_{T-}^{2-\alpha}h_{2})(0)=-\int_{0}^{T}h_{1}^{\prime}\,(D_{T-}^{\alpha-1}h_{2})\,\mathrm{d}t+h_{1}^{\prime}(T)(J_{T-}^{2-\alpha}h_{2})(T),$		(2.2)
	$\displaystyle\int_{0}^{T}h_{1}^{\prime}\,(D_{T-}^{\alpha-1}h_{2})\,\mathrm{d}t+h_{1}(0)(D_{T-}^{\alpha-1}h_{2})(0)=-\int_{0}^{T}h_{1}\,(D_{T-}^{\alpha}h_{2})\,\mathrm{d}t+h_{1}(T)(D_{T-}^{\alpha-1}h_{2})(T).$		(2.3)

Proof.

We recall the following formula connecting the forward and backward Riemann-Liouville integral operators (see [7, Lemma 4.1])

\int_{0}^{T}(J_{0+}^{\beta}h_{1})\,h_{2}\,\mathrm{d}t=\int_{0}^{T}h_{1}\,(J_{T-}^{\beta}h_{2})\,\mathrm{d}t,\quad 0<\beta\leq 1,\ h_{1},h_{2}\in C[0,T],

which can be easily verified by direct calculation. Then for $0<\alpha\leq 1$ , we employ the above formula and the usual integration by parts to derive

	$\displaystyle\int_{0}^{T}(\partial_{0+}^{\alpha}h_{1})\,h_{2}\,\mathrm{d}t$	$\displaystyle=\int_{0}^{T}(J_{0+}^{1-\alpha}h_{1}^{\prime})\,h_{2}\,\mathrm{d}t=\int_{0}^{T}h_{1}^{\prime}\,(J_{T-}^{1-\alpha}h_{2})\,\mathrm{d}t$
		$\displaystyle=\left[h_{1}\,(J_{T-}^{1-\alpha}h_{2})\right]_{0}^{T}-\int_{0}^{T}h_{1}\,(J_{T-}^{1-\alpha}h_{2})^{\prime}\,\mathrm{d}t$
		$\displaystyle=h_{1}(T)(J_{T-}^{1-\alpha}h_{2})(T)-h_{1}(0)(J_{T-}^{1-\alpha}h_{2})(0)-\int_{0}^{T}h_{1}\,(D_{T-}^{\alpha}h_{2})\,\mathrm{d}t$

for $h_{1},h_{2}\in C^{1}[0,T]$ , which gives (2.1). The derivations of (2.2)–(2.3) are similar and we omit the proofs here. ∎

Next, following the same line as [17], we specify the key assumption on the source term (1.2) and the observation subdomain $\omega$ as

\bigcup_{0\leq t\leq T}\{(\mathrm{supp}\,f+\bm{p}t)\cup(\mathrm{supp}\,g+\bm{q}t)\}\subset\subset\Omega,\quad f,g\in H_{0}^{\lceil\alpha\rceil}(\Omega),\quad\partial\omega\supset\partial\Omega.

(2.4)

In other words, it is assumed that the source term $F$ defined in (1.2) is compactly supported in $\Omega\times[0,T]$ and is sufficiently smooth. Especially, we can suppose $f,g\in H_{0}^{2}(\Omega)$ for $1<\alpha\leq 2$ since $f,g$ have compact supports in $\Omega$ . Here for simplicity, $\omega$ is also assumed to surround the whole boundary of $\Omega$ in the case of $\alpha=1$ . For readers’ better understanding, we illustrate the geometrical aspect of (2.4) in Figure 1.

Figure 1: An illustration of the geometrical assumptions on the moving source (1.2) and the subdomain

\omega

Based on the above assumptions, now we state the well-posedness results concerning (1.1)–(1.2), which are slightly modified from [17, Lemma 2.2] to fit into the framework of Sections Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations–Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations.

Lemma 2.2

Let $u$ be the solution to (1.1)–(1.2), where the profiles $f,g$ satisfy the key assumption (2.4). Then the following statements hold true.

(a) If $0<\alpha\leq 1,$ then $u\in L^{\infty}(0,T;H^{2}(\Omega)\cap H_{0}^{1}(\Omega))$ and $\partial_{t}u\in L^{\gamma}(0,T;H_{0}^{1}(\Omega)),$ where $\gamma\in(1,\frac{1}{1-\alpha})$ is arbitrarily fixed and we interpret $\frac{1}{1-\alpha}=\infty$ for $\alpha=1$ .

(b) If $1<\alpha\leq 2,$ then $u\in L^{\infty}(0,T;H^{3}(\Omega)\cap H_{0}^{1}(\Omega)),$ $\partial_{t}u\in L^{\infty}(0,T;H^{2}(\Omega)\cap H_{0}^{1}(\Omega))$ and $\partial_{t}^{2}u\in L^{\gamma}(0,T;H_{0}^{1}(\Omega)),$ where $\gamma\in(1,\frac{1}{2-\alpha})$ is arbitrarily fixed and we interpret $\frac{1}{2-\alpha}=\infty$ for $\alpha=2$ .

Compared with [17, Lemma 2.2], the time regularity of $\partial_{t}^{\lceil\alpha\rceil}u$ is improved in Lemma 2.2. Such an improvement is straightforward and we omit the detail here.

Finally, we close this section by recalling the uniqueness result for Problem 1.1 obtained in [17, Theorem 2.4]. Thanks to the linearity of Problem 1.1, it suffices to assume $u=0$ in $\omega\times(0,T)$ .

Lemma 2.3

Let $u$ be the solution to (1.1)–(1.2), where the profiles $f,g$ and the subdomain $\omega$ satisfy the key assumption (2.4).

(a) If $0<\alpha\leq 1,$ then $u=0$ in $\omega\times(0,T)$ implies $f\equiv 0$ in $\Omega$ .

(b) In the case of $1<\alpha\leq 2$ , we additionally require $T>2\,\mathrm{diam}(\Omega)$ if $\alpha=2$ . Then $u=0$ in $\omega\times(0,T)$ implies $f=g\equiv 0$ in $\Omega$ .

\@hangfrom\@seccntformat

sectionIterative thresholding scheme for one unknown profile

\@xsect

2.3ex plus.2ex

This section is devoted to the construction of a numerical inversion scheme for Problem 1.1 on identifying a single profile via a minimization procedure. In the basic formulation (1.2) and Problem 1.1, we assumed two profiles $f,g$ for $1<\alpha\leq 2$ , which definitely includes the case and the determination of a single profile. Nevertheless, from a numerical viewpoint, the recovery of a single profile for $1<\alpha\leq 2$ also deserves consideration at least as a prototype of that of two profiles. Therefore, independent of the theoretical formulation, in this section we will consider the initial-boundary value problem

\begin{cases}(\partial_{0+}^{\alpha}-\triangle)u(\bm{x},t)=f(\bm{x}-\bm{p}t),&(\bm{x},t)\in\Omega\times(0,T),\\ \partial_{t}^{k}u(\bm{x},0)=0\ (k=0,\lceil\alpha\rceil-1),&\bm{x}\in\Omega,\\ u(\bm{x},t)=0,&(\bm{x},t)\in\partial\Omega\times(0,T),\end{cases}

(3.1)

where $0<\alpha\leq 2$ and $f\in H_{0}^{1}(\Omega)$ . To emphasize the dependence, we will denote the solution to (3.1) by $u(f)$ .

Let $f_{\mathrm{true}}\in H_{0}^{1}(\Omega)$ be the true profile and $u^{\delta}$ be the observation data. For technical convenience, we assume that there exists a $\gamma\in(1,\frac{1}{\lceil\alpha\rceil-\alpha})$ such that

u^{\delta}\in W^{1,\gamma}(0,T;H^{1}(\omega)),\quad\|u(f_{\mathrm{true}})-u^{\delta}\|_{W^{1,\gamma}(0,T;H^{1}(\omega)))}\leq\delta,

(3.2)

which means that the noise in the observation data is bounded by $\delta$ in $W^{1,\gamma}(0,T;H^{1}(\omega))$ . Although (3.2) looks restrictive for observation data, it is essential for the establishment of the inversion scheme in the sequel. In view of Lemma 2.2, it turns out that (3.2) is tolerable. On the other hand, it suffices to add a preconditioning procedure to mollify the noisy data by some stable numerical differentiation method (see e.g. [25]).

\@hangfrom\@seccntformat

subsectionCase of $0<\alpha\leq 1$

\@xsect

1.5ex plus.2ex

First we consider the case of $0<\alpha\leq 1$ . Similarly to the proof of [17, Theorem 2.4], we introduce an auxiliary function $v(f):=J_{0+}^{1-\alpha}(\partial_{t}+\bm{p}\cdot\nabla)u(f)$ , which satisfies the following initial-boundary value problem for a (time-fractional) diffusion equation with a Caputo derivative in time:

\begin{cases}(\partial_{0+}^{\alpha}-\triangle)v=0&\mbox{in }\Omega\times(0,T),\\ v=f&\mbox{in }\Omega\times\{0\},\\ v=J_{0+}^{1-\alpha}(\bm{p}\cdot\nabla u(f))&\mbox{on }\partial\Omega\times(0,T).\end{cases}

Here the boundary condition of $v(f)$ follows from that of $u(f)$ and the assumption $\partial\omega\supset\partial\Omega$ . Especially, on $\partial\Omega\times(0,T)$ , we know $v(f_{\mathrm{true}})=J_{0+}^{1-\alpha}(\bm{p}\cdot\nabla u(f_{\mathrm{true}}))$ is approximated by $J_{0+}^{1-\alpha}(\bm{p}\cdot\nabla u^{\delta})$ . Then for any $f\in H_{0}^{1}(\Omega)$ , we can freeze the boundary condition of $v(f)$ as $J_{0+}^{1-\alpha}(\bm{p}\cdot\nabla u^{\delta})|_{\partial\Omega\times(0,T)}$ , i.e., $v(f)$ satisfies

\begin{cases}(\partial_{0+}^{\alpha}-\triangle)v=0&\mbox{in }\Omega\times(0,T),\\ v=f&\mbox{in }\Omega\times\{0\},\\ v=J_{0+}^{1-\alpha}(\bm{p}\cdot\nabla u^{\delta})&\mbox{on }\partial\Omega\times(0,T).\end{cases}

(3.3)

According to Lemma 2.2(a) and Hölder’s inequality, it is readily seen that $v(f)\in L^{\infty}(0,T;H^{1}(\Omega))$ for any $f\in H_{0}^{1}(\Omega)$ .

Now we set $v^{\delta}:=J_{0+}^{1-\alpha}(\partial_{t}+\bm{p}\cdot\nabla)u^{\delta}\in L^{\infty}(0,T;L^{2}(\omega))$ and consider the following regularized minimization problem with a penalty term

\min_{f\in H_{0}^{1}(\Omega)}\Phi(f):=\left\|v(f)-v^{\delta}\right\|_{L^{2}(\omega\times(0,T))}^{2}+\kappa\|\nabla f\|_{L^{2}(\Omega)}^{2},

(3.4)

where $\kappa>0$ is the regularization parameter. Here we notice that Poincaré inequality guarantees the equivalence of $\|f\|_{H^{1}(\Omega)}$ and $\|\nabla f\|_{L^{2}(\Omega)}$ . Next, we calculate the Fréchet derivative of $\Phi(f)$ for $f\in H_{0}^{1}(\Omega)$ on the direction $\widetilde{f}\in H_{0}^{1}(\Omega)$ . Picking a small $\varepsilon>0$ , we calculate

	$\displaystyle\frac{\Phi(f+\varepsilon\,\widetilde{f}\,)-\Phi(f)}{\varepsilon}$	$\displaystyle=\int_{0}^{T}\!\!\!\int_{\omega}\frac{v(f+\varepsilon\,\widetilde{f}\,)-v(f)}{\varepsilon}\left(v(f+\varepsilon\,\widetilde{f}\,)+v(f)-2\,v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t$
		$\displaystyle\quad\,+\kappa\int_{\Omega}\nabla\widetilde{f}\cdot(2\nabla f+\varepsilon\nabla\widetilde{f}\,)\,\mathrm{d}\bm{x}.$

Taking difference between the initial-boundary value problems of $v(f+\varepsilon\,\widetilde{f}\,)$ and $v(f)$ , we see that a further auxiliary function $w(\widetilde{f}\,):=(v(f+\varepsilon\,\widetilde{f}\,)-v(f))/\varepsilon$ satisfies

\begin{cases}(\partial_{0+}^{\alpha}-\triangle)w=0&\mbox{in }\Omega\times(0,T),\\ w=\widetilde{f}&\mbox{in }\Omega\times\{0\},\\ w=0&\mbox{on }\partial\Omega\times(0,T).\end{cases}

(3.5)

By linearity, we have $w(\widetilde{f}\,)\in L^{\infty}(0,T;H_{0}^{1}(\Omega))$ . Then we deduce

$\displaystyle\frac{\Phi^{\prime}(f)\widetilde{f}}{2}$	$\displaystyle=\lim_{\varepsilon\to 0}\frac{\Phi(f+\varepsilon\,\widetilde{f}\,)-\Phi(f)}{2\varepsilon}$
	$\displaystyle=\frac{1}{2}\int_{0}^{T}\!\!\!\int_{\omega}w(\widetilde{f}\,)\left(\lim_{\varepsilon\to 0}v(f+\varepsilon\,\widetilde{f}\,)+v(f)-2\,v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t+\kappa\int_{\Omega}\nabla\widetilde{f}\cdot\nabla f\,\mathrm{d}\bm{x}$
	$\displaystyle=\int_{0}^{T}\!\!\!\int_{\omega}w(\widetilde{f}\,)\left(v(f)-v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t-\kappa\int_{\Omega}\widetilde{f}\,\triangle f\,\mathrm{d}\bm{x},$	(3.6)

where we observed $v(f+\varepsilon\,\widetilde{f}\,)\longrightarrow v(f)$ in $L^{1}(0,T;L^{2}(\omega))$ as $\varepsilon\to 0$ by the continuous dependence of the solution to (3.3) upon the initial value. Here we understand $\triangle f\in H^{-1}(\Omega)$ and that the inner product $\int_{\Omega}\widetilde{f}\,\triangle f\,\mathrm{d}\bm{x}$ stands for the duality pairing ${}_{H^{-1}}\langle\triangle f,\widetilde{f}\,\rangle_{H_{0}^{1}}$ .

Now it remains to derive the explicit form of $\Phi^{\prime}(f)$ . To this end, we introduce the backward problem

\begin{cases}(D_{T-}^{\alpha}+\triangle)y=-\chi_{\omega}\left(v(f)-v^{\delta}\right)&\mbox{in }\Omega\times(0,T),\\ J_{T-}^{1-\alpha}y=0&\mbox{in }\Omega\times\{T\},\\ y=0&\mbox{on }\partial\Omega\times(0,T),\end{cases}

(3.7)

where $\chi_{\omega}$ is the characteristic function of $\omega$ . As before, we denote the solution to (3.7) as $y(f)$ to emphasize its dependence on $f$ . For the theoretical analysis and the numerical solution of (3.7) in the case of $0<\alpha<1$ , one difficulty is the treatment of the terminal condition, which involves the backward Riemann-Liouville integral as $t\to T$ . To circumvent this, we further introduce $z(f):=J_{T-}^{1-\alpha}y(f)$ , which obviously satisfies

\begin{cases}(\partial_{T-}^{\alpha}-\triangle)z=-\chi_{\omega}\,J_{T-}^{1-\alpha}\left(v(f)-v^{\delta}\right)&\mbox{in }\Omega\times(0,T),\\ z=0&\mbox{in }\Omega\times\{T\},\\ z=0&\mbox{on }\partial\Omega\times(0,T).\end{cases}

(3.8)

Since $v(f)-v^{\delta}\in L^{\infty}(0,T;L^{2}(\omega))$ , the source term of (3.8) belongs to $L^{\gamma}(0,T;L^{2}(\Omega))$ for any $\gamma\in(1,\frac{1}{1-\alpha})$ . Similarly to the proof of Lemma 2.2(a), one can verify $\partial_{t}z(f)\in L^{\gamma}(0,T;H_{0}^{1}(\Omega))$ . Hence, differentiating $z(f)=J_{T-}^{1-\alpha}y(f)$ indicates $D_{T-}^{\alpha}y(f)=\partial_{t}z(f)\in L^{\gamma}(0,T;H_{0}^{1}(\Omega))$ for any $\gamma\in(1,\frac{1}{1-\alpha})$ .

To proceed, we shall turn to the variational method to define the weak solutions to (3.5) and (3.7). Motivated by the definition of weak solutions to traditional evolution equations e.g. in [3], we employ formula (2.1) in Lemma 2.1 to provide the following definition.

Definition 3.1

Let $0<\alpha\leq 1$ and $f,\widetilde{f}\in H_{0}^{1}(\Omega)$ .

(a) We say that $w(\widetilde{f}\,)\in L^{\infty}(0,T;H_{0}^{1}(\Omega))$ is a weak solution to (3.5) if

\int_{0}^{T}\!\!\!\int_{\Omega}\left(\nabla w(\widetilde{f}\,)\cdot\nabla y-w(\widetilde{f}\,)\,(D_{T-}^{\alpha}y)\right)\mathrm{d}\bm{x}\mathrm{d}t=\int_{\Omega}\widetilde{f}\,(J_{T-}^{1-\alpha}y)(\,\cdot\,,0)\,\mathrm{d}\bm{x}

holds for all test functions $y\in L^{1}(0,T;H_{0}^{1}(\Omega))$ satisfying $D_{T-}^{\alpha}y\in L^{1}(0,T;L^{2}(\Omega))$ and $J_{T-}^{1-\alpha}y=0$ in $\Omega\times\{T\}$ .

(b) For some fixed $\gamma\in(1,\frac{1}{1-\alpha}),$ we say that $y(f)\in L^{\gamma}(0,T;H_{0}^{1}(\Omega))$ satisfying $D_{T-}^{\alpha}y(f)\in L^{\gamma}(0,T;L^{2}(\Omega))$ is a weak solution to (3.7) if $J_{T-}^{1-\alpha}y(f)=0$ in $\Omega\times\{T\}$ and

\int_{0}^{T}\!\!\!\int_{\Omega}\left(\nabla w\cdot\nabla y(f)-w\,(D_{T-}^{\alpha}y(f))\right)\mathrm{d}\bm{x}\mathrm{d}t=\int_{0}^{T}\!\!\!\int_{\omega}w\left(v(f)-v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t

holds for all test functions $w\in L^{\infty}(0,T;H_{0}^{1}(\Omega))$ .

In the above definition, the weak solution to the backward problem (3.7) is slightly stronger than that to the forward problem (3.5) because $D_{T-}^{\alpha}y(f)$ makes sense in $L^{\gamma}(0,T;L^{2}(\Omega))$ . Now we can take $y=y(f)$ and $w=w(\widetilde{f}\,)$ as mutual test functions in Definition 3.1 to calculate

	$\displaystyle\int_{0}^{T}\!\!\!\int_{\omega}w(\widetilde{f}\,)\left(v(f)-v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t$	$\displaystyle=\int_{0}^{T}\!\!\!\int_{\Omega}\left(\nabla w(\widetilde{f}\,)\cdot\nabla y(f)-w(\widetilde{f}\,)\,(D_{T-}^{\alpha}y(f))\right)\mathrm{d}\bm{x}\mathrm{d}t$
		$\displaystyle=\int_{\Omega}\widetilde{f}\,(J_{T-}^{1-\alpha}y(f))(\,\cdot\,,0)\,\mathrm{d}\bm{x}=\int_{\Omega}\widetilde{f}\,z(f)(\,\cdot\,,0)\,\mathrm{d}\bm{x}.$

Plugging the above equality into (3.6), we obtain

\frac{\Phi^{\prime}(f)\widetilde{f}}{2}=\int_{\Omega}\widetilde{f}\,(z(f)(\,\cdot\,,0)-\kappa\,\triangle f)\,\mathrm{d}\bm{x}.

Since $\widetilde{f}\in H_{0}^{1}(\Omega)$ is chosen arbitrarily, it follows from the variational principle that the minimizer $f_{*}\in H_{0}^{1}(\Omega)$ of (3.4) is the weak solution to the following boundary value problem for an elliptic equation

\begin{cases}\kappa\,\triangle f_{*}=z(f_{*})(\,\cdot\,,0)&\mbox{in }\Omega,\\ f_{*}=0&\mbox{on }\partial\Omega.\end{cases}

(3.9)

Consequently, in the same manner as [7], we can propose the following iterative thresholding update

\left\{\!\begin{aligned} &\triangle f_{\ell+1}=\frac{1}{M+\kappa}z(f_{\ell})(\,\cdot\,,0)+\frac{M}{M+\kappa}\triangle f_{\ell}&\quad&\mbox{in }\Omega,\\ &f_{\ell+1}=0&\quad&\mbox{on }\partial\Omega,\end{aligned}\right.

(3.10)

where $M>0$ is a tuning parameter. Now we summarize the numerical reconstruction scheme for Problem 1.1 with $0<\alpha\leq 1$ as follows.

Algorithm 3.2

Choose a tolerance $\epsilon>0$ , a regularization parameter $\kappa>0$ and a tuning parameter $M>0$ . Give an initial guess $f_{0}\in H_{0}^{1}(\Omega)$ (e.g., $f_{0}\equiv 0$ ), and set $\ell=0$ .

1.

Compute $f_{\ell+1}$ by the iterative update (3.10).
2.

If $\|f_{\ell+1}-f_{\ell}\|_{H^{1}(\Omega)}/\|f_{\ell}\|_{H^{1}(\Omega)}<\epsilon$ , stop the iteration. Otherwise, update $\ell\leftarrow\ell+1$ and return to Step 1.

By choosing $M>0$ suitably large, the convergence of the sequence $\{f_{\ell}\}_{\ell=0}^{\infty}\subset H_{0}^{1}(\Omega)$ is guaranteed by [2]. In each step of the iteration, it suffices to solve $3$ differential equations, i.e., the modified forward problem (3.3) for $v(f_{\ell})$ , the backward problem (3.8) for $z(f_{\ell})$ , and the boundary value problem for the Poisson equation (3.10) for $f_{\ell+1}$ .

Remark 3.3

Similarly, if we penalize the $L^{2}$ -norm of $f$ instead of $\nabla f$ in the definition (3.4) of $\Phi(f)$ , then the resulting iterative thresholding update turns out to be

f_{\ell+1}=\frac{1}{M+\kappa}z(f_{\ell})(\,\cdot\,,0)+\frac{M}{M+\kappa}f_{\ell},

which differs from the scheme (3.10) only on the cost of solving a Poisson equation numerically. Nevertheless, since we assumed $f\in H_{0}^{1}(\Omega)$ , such sacrifice on the computational cost is unimportant in view of pursuing more accurate reconstruction of the profile in the $H^{1}$ sense.

\@hangfrom\@seccntformat

subsectionCase of $1<\alpha\leq 2$

\@xsect

1.5ex plus.2ex

Now we consider (3.1) in the case of $1<\alpha\leq 2$ with $f\in H_{0}^{1}(\Omega)$ . Analogously to the proof of Lemma 2.2(b), one can easily check

u(f)\in L^{\infty}(0,T;H^{2}(\Omega)\cap H_{0}^{1}(\Omega)),\quad\partial_{t}u(f)\in L^{\gamma}(0,T;H_{0}^{1}(\Omega)),

where $\gamma\in(1,\frac{1}{2-\alpha})$ . Again setting $v(f):=J_{0+}^{2-\alpha}(\partial_{t}+\bm{p}\cdot\nabla)u(f)$ , $v^{\delta}:=J_{0+}^{2-\alpha}(\partial_{t}+\bm{p}\cdot\nabla)u^{\delta}$ and freezing the boundary condition of $v(f)$ as $J_{0+}^{2-\alpha}(\bm{p}\cdot\nabla u^{\delta})$ , we know $v(f)\in W^{1,\infty}(0,T;H^{1}(\Omega))$ , $v^{\delta}\in L^{\infty}(0,T;L^{2}(\omega))$ and $v(f)$ satisfies

\begin{cases}(\partial_{0+}^{\alpha}-\triangle)v=0&\mbox{in }\Omega\times(0,T),\\ v=0,\ \partial_{t}v=f&\mbox{in }\Omega\times\{0\},\\ v=J_{0+}^{2-\alpha}(\bm{p}\cdot\nabla u^{\delta})&\mbox{on }\partial\Omega\times(0,T).\end{cases}

Treating the same minimization problem (3.4) and repeating the calculation of the Fréchet derivative, again we obtain (3.6), where $w(\widetilde{f})\in W^{1,\infty}(0,T;H_{0}^{1}(\Omega))$ satisfies

\begin{cases}(\partial_{0+}^{\alpha}-\triangle)w=0&\mbox{in }\Omega\times(0,T),\\ w=0,\ \partial_{t}w=\widetilde{f}&\mbox{in }\Omega\times\{0\},\\ w=0&\mbox{on }\partial\Omega\times(0,T).\end{cases}

(3.11)

In a parallel manner as before, we investigate the solution $y(f)$ to the backward problem

\begin{cases}(D_{T-}^{\alpha}-\triangle)y=\chi_{\omega}\left(v(f)-v^{\delta}\right)&\mbox{in }\Omega\times(0,T),\\ J_{T-}^{2-\alpha}y=D_{T-}^{\alpha-1}y=0&\mbox{in }\Omega\times\{T\},\\ y=0&\mbox{on }\partial\Omega\times(0,T)\end{cases}

(3.12)

and the further auxiliary function $z(f):=J_{T-}^{2-\alpha}y(f)$ satisfying

\begin{cases}(\partial_{T-}^{\alpha}-\triangle)z=\chi_{\omega}\,J_{T-}^{2-\alpha}\left(v(f)-v^{\delta}\right)&\mbox{in }\Omega\times(0,T),\\ z=\partial_{t}z=0&\mbox{in }\Omega\times\{T\},\\ z=0&\mbox{on }\partial\Omega\times(0,T).\end{cases}

By a similar argument as that in the previous subsection, one can verify $y(f),D_{T-}^{\alpha-1}y(f)\in L^{\gamma}(0,T;H_{0}^{1}(\Omega))$ with some $\gamma\in(1,\frac{1}{2-\alpha})$ . Next, we employ formulae (2.2) and (2.3) in Lemma 2.1(b) to define the weak solutions to (3.11) and (3.12), respectively.

Definition 3.4

Let $1<\alpha\leq 2$ and $f,\widetilde{f}\in H_{0}^{1}(\Omega)$ .

(a) We say that $w(\widetilde{f}\,)\in L^{\infty}(0,T;H_{0}^{1}(\Omega))\cap W^{1,\infty}(0,T;L^{2}(\Omega))$ is a weak solution to (3.11) if

\int_{0}^{T}\!\!\!\int_{\Omega}\left\{\nabla w(\widetilde{f}\,)\cdot\nabla y-\left(\partial_{t}w(\widetilde{f}\,)\right)(D_{T-}^{\alpha-1}y)\right\}\mathrm{d}\bm{x}\mathrm{d}t=\int_{\Omega}\widetilde{f}\,(J_{T-}^{2-\alpha}y)(\,\cdot\,,0)\,\mathrm{d}\bm{x}

holds for all test functions $y\in L^{1}(0,T;H_{0}^{1}(\Omega))$ satisfying $D_{T-}^{\alpha-1}y\in L^{1}(0,T;L^{2}(\Omega))$ and $J_{T-}^{2-\alpha}y=0$ in $\Omega\times\{T\}$ .

(b) For some fixed $\gamma\in(1,\frac{1}{2-\alpha}),$ we say that $y(f)\in L^{\gamma}(0,T;H_{0}^{1}(\Omega))$ satisfying $D_{T-}^{\alpha-1}y(f)\in L^{\gamma}(0,T;L^{2}(\Omega))$ is a weak solution to (3.12) if $J_{T-}^{2-\alpha}y(f)=D_{T-}^{\alpha-1}y(f)=0$ in $\Omega\times\{T\}$ and

\int_{0}^{T}\!\!\!\int_{\Omega}\left(\nabla w\cdot\nabla y(f)-(\partial_{t}w)\,(D_{T-}^{\alpha-1}y(f))\right)\mathrm{d}\bm{x}\mathrm{d}t=\int_{0}^{T}\!\!\!\int_{\omega}w\left(v(f)-v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t

holds for all test functions $w\in L^{\infty}(0,T;H_{0}^{1}(\Omega))\cap W^{1,\infty}(0,T;L^{2}(\Omega))$ satisfying $w(f)=0$ in $\Omega\times\{0\}$ .

Taking $w(\widetilde{f}\,)$ and $y(f)$ as mutual test functions in the above definition, again we obtain

	$\displaystyle\int_{0}^{T}\!\!\!\int_{\omega}w(\widetilde{f}\,)\left(v(f)-v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t$	$\displaystyle=\int_{0}^{T}\!\!\!\int_{\Omega}\left\{\nabla w(\widetilde{f}\,)\cdot\nabla y(f)-\left(\partial_{t}w(\widetilde{f}\,)\right)(D_{T-}^{\alpha-1}y(f))\right\}\mathrm{d}\bm{x}\mathrm{d}t$
		$\displaystyle=\int_{\Omega}\widetilde{f}\,(J_{T-}^{2-\alpha}y(f))(\,\cdot\,,0)\,\mathrm{d}\bm{x}=\int_{\Omega}\widetilde{f}\,z(f)(\,\cdot\,,0)\,\mathrm{d}\bm{x}.$

As a result, formally we arrive at the same Euler-Lagrange equation (3.9) and thus the thresholding iterative update (3.10).

\@hangfrom\@seccntformat

sectionReconstruction scheme for two unknown profiles

\@xsect

2.3ex plus.2ex

On the same direction as that of the previous section, in this section we attempt to develop a numerical reconstruction scheme for Problem 1.1 in the case of two unknown profiles. More precisely, let $1<\alpha\leq 2$ , $f,g\in H_{0}^{2}(\Omega)$ and denote by $u(f,g)$ the solution to (1.1)–(1.2). Similarly as before, by $f_{\mathrm{true}},g_{\mathrm{true}}\in H_{0}^{2}(\Omega)$ and $u^{\delta}$ we denote the true profiles and the observation data, respectively. As a generalization of (3.2), here we assume that there exists $\gamma\in(1,\frac{1}{2-\alpha})$ such that

u^{\delta}\in W^{2,\gamma}(0,T;H^{2}(\omega)),\quad\|u(f_{\mathrm{true}},g_{\mathrm{true}})-u^{\delta}\|_{W^{2,\gamma}(0,T;H^{2}(\omega))}\leq\delta.

Meanwhile, introducing

v(f,g):=J_{0+}^{2-\alpha}(\partial_{t}+\bm{p}\cdot\nabla)(\partial_{t}+\bm{q}\cdot\nabla)u(f,g),\quad v^{\delta}:=J_{0+}^{2-\alpha}(\partial_{t}+\bm{p}\cdot\nabla)(\partial_{t}+\bm{q}\cdot\nabla)u^{\delta}

and freezing the boundary condition of $v(f,g)$ as

v^{\delta}=(\bm{p}+\bm{q})\cdot\nabla(\partial_{0+}^{\alpha-1}u^{\delta})+\bm{p}\cdot\nabla(\bm{q}\cdot\nabla(J_{0+}^{2-\alpha}u^{\delta}))\quad\mbox{on }\partial\Omega\times(0,T),

we see that $v^{\delta}\in L^{\infty}(0,T;L^{2}(\omega))$ and it was shown in [17, §4.2] that $v(f,g)$ satisfies

\begin{cases}(\partial_{0+}^{\alpha}-\triangle)v=0&\mbox{in }\Omega\times(0,T),\\ v=f+g,\ \partial_{t}v=\bm{q}\cdot\nabla f+\bm{p}\cdot\nabla g&\mbox{in }\Omega\times\{0\},\\ v=v^{\delta}&\mbox{on }\partial\Omega\times(0,T).\end{cases}

(4.1)

Then it follows from Lemma 2.2(b) and Hölder’s inequality that $v(f,g)\in L^{\infty}(0,T;H^{1}(\Omega))$ for any $f,g\in H_{0}^{2}(\Omega)$ .

Since the unknown profiles $f,g$ appears in the initial values of (4.1), it is natural to divide the reconstruction into two steps, i.e., the simultaneous recovery of the initial values and then the determination of $f,g$ by the information of $f+g$ and $\bm{q}\cdot\nabla f+\bm{p}\cdot\nabla g$ . To this end, it is convenient to rewrite the initial condition of (4.1) as

v(\,\cdot\,,0)=a:=f+g\in H_{0}^{2}(\Omega),\quad\partial_{t}v(\,\cdot\,,0)=b:=\bm{q}\cdot\nabla f+\bm{p}\cdot\nabla g\in H_{0}^{1}(\Omega)

(4.2)

and denote the solution to (4.1) as $v(a,b)$ .

\@hangfrom\@seccntformat

subsectionIterative thresholding scheme for the initial values

\@xsect

1.5ex plus.2ex

As a natural generalization of the minimization problem (3.4), we consider the following multivariable minimization problem with two penalty terms

\min_{(a,b)\in H_{0}^{2}(\Omega)\times H_{0}^{1}(\Omega)}\Psi(a,b):=\left\|v(a,b)-v^{\delta}\right\|_{L^{2}(\omega\times(0,T))}^{2}+\kappa_{2}\|\triangle a\|_{L^{2}(\Omega)}^{2}+\kappa_{1}\|\nabla b\|_{L^{2}(\Omega)}^{2},

(4.3)

where $\kappa_{1},\kappa_{2}>0$ are regularization parameters. Here, owing to $a\in H_{0}^{2}(\Omega)$ and $b\in H_{0}^{1}(\Omega)$ , we know that $\|\triangle a\|_{L^{2}(\Omega)}$ and $\|\nabla b\|_{L^{2}(\Omega)}$ are equivalent to $\|a\|_{H^{2}(\Omega)}$ and $\|b\|_{H^{1}(\Omega)}$ , respectively.

To calculate the Fréchet derivative of $\Psi$ for $(a,b)\in H_{0}^{2}(\Omega)\times H_{0}^{1}(\Omega)$ on the direction $(\widetilde{a},\widetilde{b}\,)\in H_{0}^{2}(\Omega)\times H_{0}^{1}(\Omega)$ , we pick a small $\varepsilon>0$ to calculate

	$\displaystyle\frac{\Psi(a+\varepsilon\,\widetilde{a},b+\varepsilon\,\widetilde{b}\,)-\Psi(a,b)}{\varepsilon}$	$\displaystyle=\int_{0}^{T}\!\!\!\int_{\omega}w(\widetilde{a},\widetilde{b}\,)\left\{v(a+\varepsilon\,\widetilde{a},b+\varepsilon\,\widetilde{b}\,)+v(a,b)-2\,v^{\delta}\right\}\mathrm{d}\bm{x}\mathrm{d}t$
		$\displaystyle\quad\,+\kappa_{2}\int_{\Omega}\triangle\widetilde{a}\,\triangle(2a+\varepsilon\,\widetilde{a}\,)\,\mathrm{d}\bm{x}+\kappa_{1}\int_{\Omega}\nabla\widetilde{b}\cdot\nabla(2b+\varepsilon\,\widetilde{b}\,)\,\mathrm{d}\bm{x},$

where $w(\widetilde{a},\widetilde{b}\,):=(v(a+\varepsilon\,\widetilde{a},b+\varepsilon\,\widetilde{b}\,)-v(a,b))/\varepsilon\in L^{\infty}(0,T;H_{0}^{1}(\Omega))$ satisfies

\begin{cases}(\partial_{0+}^{\alpha}-\triangle)w=0&\mbox{in }\Omega\times(0,T),\\ w=\widetilde{a},\ \partial_{t}w=\widetilde{b}&\mbox{in }\Omega\times\{0\},\\ w=0&\mbox{on }\partial\Omega\times(0,T).\end{cases}

(4.4)

Then by passing $\varepsilon\to 0$ and integration by parts, we deduce

$\displaystyle\frac{\nabla\Psi(a,b)\cdot(\widetilde{a},\widetilde{b}\,)}{2}$	$\displaystyle=\lim_{\varepsilon\to 0}\frac{\Psi(a+\varepsilon\,\widetilde{a},b+\varepsilon\,\widetilde{b}\,)-\Psi(a,b)}{2\varepsilon}$
	$\displaystyle=\int_{0}^{T}\!\!\!\int_{\omega}w(\widetilde{a},\widetilde{b}\,)\left(v(a,b)-v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t+\kappa_{2}\int_{\Omega}\triangle\widetilde{a}\,\triangle a\,\mathrm{d}\bm{x}+\kappa_{1}\int_{\Omega}\nabla\widetilde{b}\cdot\nabla b\,\mathrm{d}\bm{x}$
	$\displaystyle=\int_{0}^{T}\!\!\!\int_{\Omega}w(\widetilde{a},\widetilde{b}\,)\left\{\chi_{\omega}\left(v(a,b)-v^{\delta}\right)\right\}\mathrm{d}\bm{x}\mathrm{d}t+\int_{\Omega}\left(\kappa_{2}\widetilde{a}\,\triangle^{2}a-\kappa_{1}\widetilde{b}\,\triangle b\right)\mathrm{d}\bm{x}$
	$\displaystyle=\int_{\Omega}\left[\int_{0}^{T}w(\widetilde{a},\widetilde{b}\,)\left\{\chi_{\omega}\left(v(a,b)-v^{\delta}\right)\right\}\mathrm{d}t+\begin{pmatrix}\widetilde{a}\,\\ \widetilde{b}\end{pmatrix}\cdot\begin{pmatrix}\kappa_{2}\triangle^{2}a\\ -\kappa_{1}\triangle b\end{pmatrix}\right]\mathrm{d}\bm{x},$	(4.5)

where we employed $v(a+\varepsilon\,\widetilde{a},b+\beta\,\widetilde{b}\,)\longrightarrow v(a,b)$ in $L^{1}(0,T;L^{2}(\omega))$ as $\varepsilon\to 0$ , and interpret $\triangle^{2}a\in H^{-2}(\Omega)$ and $\triangle b\in H^{-1}(\Omega)$ for convenience.

In order to derive the Euler-Lagrange equation for the minimizer of (4.3), again we shall turn to the solution $y(a,b)$ to the backward problem

\begin{cases}(D_{T-}^{\alpha}-\triangle)y=\chi_{\omega}\left(v(a,b)-v^{\delta}\right)&\mbox{in }\Omega\times(0,T),\\ J_{T-}^{2-\alpha}y=D_{T-}^{\alpha-1}y=0&\mbox{in }\Omega\times\{T\},\\ y=0&\mbox{on }\partial\Omega\times(0,T)\end{cases}

(4.6)

as well as its Caputo counterpart $z(a,b):=J_{T-}^{2-\alpha}y(a,b)$ satisfying

\begin{cases}(\partial_{T-}^{\alpha}-\triangle)z=\chi_{\omega}\,J_{T-}^{2-\alpha}\left(v(a,b)-v^{\delta}\right)&\mbox{in }\Omega\times(0,T),\\ z=\partial_{t}z=0&\mbox{in }\Omega\times\{T\},\\ z=0&\mbox{on }\partial\Omega\times(0,T).\end{cases}

(4.7)

Similarly to the argument for (3.7) and (3.8), one can employ Lemma 2.2(b) to verify that

y(a,b)=\partial_{T-}^{\alpha-1}z(a,b)\in L^{\gamma}(0,T;H_{0}^{1}(\Omega)),\quad D_{T-}^{\alpha}y(a,b)=\partial_{t}^{2}z(a,b)\in L^{\gamma}(0,T,L^{2}(\Omega))

for any $\gamma\in(1,\frac{1}{2-\alpha})$ .

Now we are well prepared to provide the definitions of the weak solutions to (4.4) and (4.6).

Definition 4.1

Let $1<\alpha\leq 2$ and $a,b,\widetilde{a},\widetilde{b}\in H_{0}^{2}(\Omega)$ .

(a) We say that $w(\widetilde{a},\widetilde{b}\,)\in L^{\infty}(0,T;H_{0}^{1}(\Omega))$ is a weak solution to (4.4) if

\int_{0}^{T}\!\!\!\int_{\Omega}\left(\nabla w(\widetilde{a},\widetilde{b}\,)\cdot\nabla y+w(\widetilde{a},\widetilde{b}\,)\,(D_{T-}^{\alpha}y)\right)\mathrm{d}\bm{x}\mathrm{d}t=\int_{\Omega}\left(\widetilde{b}\,J_{T-}^{2-\alpha}y-\widetilde{a}\,D_{T-}^{\alpha-1}y\right)(\,\cdot\,,0)\,\mathrm{d}\bm{x}

holds for all test functions $y\in L^{1}(0,T;H_{0}^{1}(\Omega))$ satisfying $D_{T-}^{\alpha}y\in L^{1}(0,T;L^{2}(\Omega))$ and $J_{T-}^{2-\alpha}y=D_{T-}^{\alpha-1}y=0$ in $\Omega\times\{T\}$ .

(b) For some fixed $\gamma\in(1,\frac{1}{2-\alpha}),$ we say that $y(a,b)\in L^{\gamma}(0,T;H_{0}^{1}(\Omega))$ satisfying $D_{T-}^{\alpha}y(a,b)\in L^{\gamma}(0,T;L^{2}(\Omega))$ is a weak solution to (4.6) if $J_{T-}^{2-\alpha}y(a,b)=D_{T-}^{\alpha-1}y(a,b)=0$ in $\Omega\times\{T\}$ and

\int_{0}^{T}\!\!\!\int_{\Omega}\left(\nabla w\cdot\nabla y(a,b)+w\,(D_{T-}^{\alpha}y(a,b))\right)\mathrm{d}\bm{x}\mathrm{d}t=\int_{0}^{T}\!\!\!\int_{\omega}w\left(v(a,b)-v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t

holds for all test functions $w\in L^{\infty}(0,T;H_{0}^{1}(\Omega))$ .

In comparison with Definition 3.4, the above definition requires lower time regularity for the solution to the forward problem (4.4), whereas requires higher time regularity for that to the backward problem (4.6). As before, it is routine to take $y(a,b)$ and $w(\widetilde{a},\widetilde{b}\,)$ as mutual test functions to deduce

	$\displaystyle\int_{0}^{T}\!\!\!\int_{\omega}w(\widetilde{a},\widetilde{b}\,)\left(v(a,b)-v^{\delta}\right)\mathrm{d}\bm{x}\mathrm{d}t$	$\displaystyle=\int_{0}^{T}\!\!\!\int_{\Omega}\left(\nabla w(\widetilde{a},\widetilde{b}\,)\cdot\nabla y(a,b)+w(\widetilde{a},\widetilde{b}\,)\,(D_{T-}^{\alpha}y(a,b))\right)\mathrm{d}\bm{x}\mathrm{d}t$
		$\displaystyle=\int_{\Omega}\left(\widetilde{b}\,J_{T-}^{2-\alpha}y(a,b)-\widetilde{a}\,D_{T-}^{\alpha-1}y(a,b)\right)(\,\cdot\,,0)\,\mathrm{d}\bm{x}$
		$\displaystyle=\int_{\Omega}\begin{pmatrix}\widetilde{a}\,\\ \widetilde{b}\end{pmatrix}\cdot\begin{pmatrix}-\partial_{t}z(a,b)\\ z(a,b)\end{pmatrix}(\,\cdot\,,0)\,\mathrm{d}\bm{x}.$

Substituting the above equality into (4.5), we arrive at

\frac{\nabla\Psi(a,b)\cdot(\widetilde{a},\widetilde{b}\,)}{2}=\int_{\Omega}\begin{pmatrix}\widetilde{a}\,\\ \widetilde{b}\end{pmatrix}\cdot\begin{pmatrix}-\partial_{t}z(a,b)(\,\cdot\,,0)+\kappa_{2}\triangle^{2}a\\ z(a,b)(\,\cdot\,,0)-\kappa_{1}\triangle b\end{pmatrix}\mathrm{d}\bm{x}.

In other words, the minimizer $(a_{*},b_{*})\in H_{0}^{2}(\Omega)\times H_{0}^{1}(\Omega)$ of (4.3) is the weak solution to the following Euler-Lagrange equation

\begin{cases}\kappa_{2}\triangle^{2}a_{*}=\partial_{t}z(a_{*},b_{*})(\,\cdot\,,0),\\ \kappa_{1}\triangle b_{*}=z(a_{*},b_{*})(\,\cdot\,,0)\end{cases}\mbox{in }\Omega.

Especially, $a_{*}$ satisfies a bi-Laplace equation in the sense of $H^{-2}(\Omega)$ . Therefore, in a similar manner as before, we can design the following iterative thresholding update to produce a sequence $\{(a_{\ell},b_{\ell})\}\subset H_{0}^{2}(\Omega)\times H_{0}^{1}(\Omega)$ approximating $(a_{*},b_{*})$ :

\left\{\!\begin{aligned} &\triangle^{2}a_{\ell+1}=\frac{1}{M_{2}+\kappa_{2}}\partial_{t}z(a_{\ell},b_{\ell})(\,\cdot\,,0)+\frac{M_{2}}{M_{2}+\kappa_{2}}\triangle^{2}a_{\ell},\\ &\triangle b_{\ell+1}=\frac{1}{M_{1}+\kappa_{1}}z(a_{\ell},b_{\ell})(\,\cdot\,,0)+\frac{M_{1}}{M_{1}+\kappa_{1}}\triangle b_{\ell}\end{aligned}\right.\quad\mbox{in }\Omega,

(4.8)

where $M_{1},M_{2}>0$ are large parameters guaranteeing the convergence of (4.8). By terminating (4.8) appropriately, one can obtain a good approximation of $(a_{*},b_{*})$ . Parallel to the iterative update (3.10) for recovering one profile, the implementation of (4.8) only involves the solutions of $2$ evolution equations, that is, (4.1) and (4.7).

We close this subsection by summarizing the numerical reconstruction scheme for Problem 1.1 with two unknown profiles as follows.

Algorithm 4.2

Choose a tolerance $\epsilon>0$ , two regularization parameters $\kappa_{1},\kappa_{2}>0$ and two tuning parameters $M_{1},M_{2}>0$ . Give an initial guess $(a_{0},b_{0})\in H_{0}^{2}(\Omega)\times H_{0}^{1}(\Omega)$ (e.g., $a_{0}=b_{0}\equiv 0$ ), and set $\ell=0$ .

1.

Compute $(a_{\ell+1},b_{\ell+1})$ by the iterative update (4.8).
2.

If $\max\{\|a_{\ell+1}-a_{\ell}\|_{H^{2}(\Omega)}/\|a_{\ell}\|_{H^{2}(\Omega)},\|b_{\ell+1}-b_{\ell}\|_{H^{1}(\Omega)}/\|b_{\ell}\|_{H^{1}(\Omega)}\}<\epsilon$ , stop the iteration. Otherwise, update $\ell\leftarrow\ell+1$ and return to Step 1.

\@hangfrom\@seccntformat

subsectionElliptic approach to solving the convection equation

\@xsect

1.5ex plus.2ex

Supposing that the initial values of (4.1) are known, in this subsection we develop an efficient numerical scheme to reconstruct the two unknown profiles. By the relation (4.2), it suffices to solve the convection equation

\bm{r}\cdot\nabla f=c:=b+\bm{p}\cdot\nabla a\in H_{0}^{1}(\Omega),\quad f\in H_{0}^{2}(\Omega)

(4.9)

for $f$ and thus $g=a-f$ , where $\mathbb{R}^{d}\ni\bm{r}:=\bm{q}-\bm{p}\neq\bm{0}$ . Since the unknown function $f$ is sufficiently smooth and vanishes on $\partial\Omega$ , we attempt to transform the first order equation (4.9) to some second order elliptic equations for better stability.

To this end, we introduce the change of variables

\bm{\xi}=\bm{\xi}(\bm{x}):=\bm{Q}^{\mathrm{T}}\bm{x},\quad\mathrm{SO}(d)\ni\bm{Q}:=\begin{pmatrix}\displaystyle\frac{\bm{r}}{|\bm{r}|}&*\end{pmatrix},\quad\widehat{\Omega}:=\{\bm{Q}^{\mathrm{T}}\bm{x}\mid\bm{x}\in\Omega\},

(4.10)

where $*$ denotes some $d\times(d-1)$ matrix. In other words, $\bm{Q}$ is a $d\times d$ orthogonal matrix whose first column is normalized along $\bm{r}$ . Correspondingly, defining the auxiliary functions

\widehat{f}(\bm{\xi}):=f(\bm{x}(\bm{\xi}))=f(\bm{Q}\bm{\xi}),\quad\widehat{c}(\bm{\xi}):=c(\bm{x}(\bm{\xi}))=c(\bm{Q}\bm{\xi}),

we know $\widehat{f}\in H_{0}^{2}(\widehat{\Omega})$ and $\widehat{c}\in H_{0}^{1}(\widehat{\Omega})$ due to the invertibility and the smoothness of this change of variables. Especially, by (4.9) we calculate

\frac{\partial}{\partial\xi_{1}}\widehat{f}(\bm{\xi})=\frac{\partial}{\partial\xi_{1}}\left(\frac{\xi_{1}\bm{r}}{|\bm{r}|}\right)\cdot\nabla_{\bm{x}}f(\bm{Q}\bm{\xi})=\frac{\bm{r}}{|\bm{r}|}\cdot\nabla_{\bm{x}}f(\bm{Q}\bm{\xi})=\frac{c(\bm{Q}\bm{\xi})}{|\bm{r}|}=\frac{\widehat{c}(\bm{\xi})}{|\bm{r}|},\quad\bm{\xi}\in\widehat{\Omega}.

(4.11)

Noting that (4.11) only involves the first derivative in $\xi_{1}$ , we write $\bm{\xi}=(\xi_{1};\bm{\xi}^{\prime})$ , where $\bm{\xi}^{\prime}\in\mathbb{R}^{d-1}$ is regarded as a parameter. Therefore, we can understand the partial derivative $\frac{\partial}{\partial\xi_{1}}$ in (4.11) as an ordinary one, and differentiating both sides of (4.11) with respect to $\xi_{1}$ leads to a series of boundary value problems for second order ordinary differential equations with parameters $\bm{\xi}^{\prime}$ :

\left\{\!\begin{aligned} &\frac{\mathrm{d}^{2}\widehat{f}(\xi_{1};\bm{\xi}^{\prime})}{\mathrm{d}\xi_{1}^{2}}=\frac{1}{|\bm{r}|}\frac{\mathrm{d}\widehat{c}(\xi_{1};\bm{\xi}^{\prime})}{\mathrm{d}\xi_{1}},&\quad&(\xi_{1};\bm{\xi}^{\prime})\in\widehat{\Omega},\\ &\widehat{f}(\xi_{1};\bm{\xi}^{\prime})=0,&\quad&(\xi_{1},\bm{\xi}^{\prime})\in\partial\widehat{\Omega}.\end{aligned}\right.

(4.12)

Hence, to reconstruct $\widehat{f}$ in $\widehat{\Omega}$ , it suffices solve (4.12) line by line along $\xi_{1}$ -axis for any fixed $\bm{\xi}^{\prime}$ , which is expected to be numerically cheap. As long as (4.12) is solved for all $\bm{\xi}\in\widehat{\Omega}$ , one can recover $f$ in $\Omega$ by the relation $f(\bm{x})=\widehat{f}(\bm{Q}^{\mathrm{T}}\bm{x})$ and finally $g=a-f$ .

We summarize the above procedure by an illustrative example in Figure 2. First, we perform the coordinate transformation (4.10) so that the direction of the first coordinate $\xi_{1}$ coincides with that of $\bm{r}$ . Next, direct calculations yield that the convection equation (4.9) in the original coordinate system is reduced to a special one (4.11) in the new system, which is a series of first order ordinary differential equations with respect to $\xi_{1}$ . Finally, utilizing the homogeneous boundary condition in (4.9), we differentiate (4.11) in $\xi_{1}$ and arrive at a series of boundary value problems (4.12) of the second order. As is indicated by the shade in Figure 2, (4.12) can be solved along each segment in $\Omega$ on the direction of $\bm{r}$ .

Figure 2: A two-dimensional illustration of the idea of transforming the convection equation (4.9) into a series of second order ordinary differential equations (4.12).

Example 4.3

(1) For $d=1$ , it suffices to differentiate (4.9) to obtain a two-point boundary value problem

\begin{cases}r\,f^{\prime\prime}=c^{\prime}&\mbox{in }\Omega,\\ f=0&\mbox{on }\partial\Omega,\end{cases}

l which is automatically included in the framework of (4.12).

(2) For $d=2$ , consider $\Omega:=\{|\bm{x}|<1\}$ and $\bm{r}:=(1,1)^{\mathrm{T}}$ . Then by the above argument, one can choose $\displaystyle\bm{Q}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-1\\ 1&1\end{pmatrix}$ , and (4.12) becomes a series of two-point boundary value problems with the parameter $\xi_{2}\in(-1,1)$ :

\left\{\!\begin{aligned} &\frac{\mathrm{d}^{2}\widehat{f}(\xi_{1};\xi_{2})}{\mathrm{d}\xi_{1}^{2}}=\frac{1}{\sqrt{2}}\frac{\mathrm{d}\widehat{c}(\xi_{1};\xi_{2})}{\mathrm{d}\xi_{1}},&\quad&|\xi_{1}|<\sqrt{1-\xi_{2}^{2}}\,,\\ &\widehat{f}(\xi_{1};\xi_{2})=0,&\quad&|\xi_{1}|=\sqrt{1-\xi_{2}^{2}}\,.\end{aligned}\right.

\@hangfrom\@seccntformat

sectionConcluding remarks

\@xsect

2.3ex plus.2ex

Regardless of the vast difference mainly in quantitative analysis, many inverse problems for (time-fractional) evolution equations share certain similarity in methodology and qualitative properties. Accordingly, for some specified inverse problems, it seems reasonable and possible to develop universal numerical reconstruction scheme valid for any fractional orders $\alpha\in(0,2]$ . As a typical candidate, in [17] and this article, we investigated the inverse moving source problem on determining moving profile(s) in (1.1)–(1.2), which turns out to be novel for (time-fractional) evolution equations from both theoretical and numerical aspects. Guaranteed by the uniqueness claimed in [17], here we adopted a conventional minimization procedure to construct iterative thresholding schemes (Algorithms 3.2 and 4.2) for recovering one or two unknown profile(s), which follows the same line as that in [7]. However, unlike [7], here we actually proposed two independent numerical methods, namely,

1.

an iterative method for classical backward problems on determining initial values, and
2.

an elliptic method for convection equations with whole boundary conditions.

The first method seems novel for the simultaneous reconstruction of two initial values in the case of $1<\alpha<2$ , and the second one is expected to improve the numerical stability and efficiency. Especially, since usual convection equations only require the in-flow condition on a partial boundary, to the author’s knowledge there seems no literature on this direction.

Finally, we summarize several important future topics related to this article.

1.

Needless to say, the immediate task should be the implementation of the proposed algorithms. Examining the whole process, we shall apply some mollification method for differentiating the noisy data, and prepare a forward solver for $1<\alpha<2$ . Meanwhile, it would be interesting to observe the numerical performance as $\alpha\to 2$ . If the equation behaves closer to a wave equation in the sense of the finite wave propagation speed, it should require more observation time for stable recovery.
2.

Similarly to [17], the numerical schemes proposed in previous sections work for more general formulation than (1.1), e.g., an elliptic part with constant coefficients instead of $-\triangle$ . However, the key introduction [17, eq. (4.2)] of auxiliary functions fails whenever the governing operator involves variable coefficients. It would be interesting to find alternative ways to overcome this difficulty.

The assumption (1.2) on the translation of moving sources seems restrictive and mostly unrealistic. Inspired by [5], one can consider the replacement of (1.2) e.g. with

F(\bm{x},t)=\begin{cases}f(\bm{x}-\bm{\rho}(t))\,\theta(t),&0<\alpha\leq 1,\\ f(\bm{x}-\bm{\rho}(t))\,\theta(t)+g(\bm{x}-\bm{\zeta}(t))\,\eta(t),&1<\alpha\leq 2,\end{cases}

(5.1)

where $\bm{\rho},\bm{\zeta}$ denote the orbits of moving sources, and $\theta,\eta$ model the time evolution of source magnitude.

4.

As another type of inverse moving source problems, the identification of moving orbits was studied in [6], whereas its numerical reconstruction is absent except for hyperbolic equations. In the framework of (5.1), the orbit inverse problems deserve reconsideration especially from a numerical viewpoint.
5.

In [17] and this article, we considered the partial interior data and required the observation subdomain to surround the whole boundary. Similarly, one can deal with the case of full lateral Cauchy data and investigate the same problem both theoretically and numerically.
6.

On the direction of Subsection Numerical schemes for reconstructing profiles of moving sources in (time-fractional) evolution equations, it would be challenging to develop an elliptic approach to boundary value problems for convection equations with variable coefficients, i.e., $\bm{r}=\bm{r}(x)$ in (4.9). It is conjectured that if the flows determined by $\bm{r}$ do not intersect each other and possess ergodicity, then we can reduce (4.9) into a series of second order ordinary differential equations on flows. In such a way, the convection equation can be related with elliptic equations with Laplace-Beltrami operators on Riemannian manifolds.

Acknowledgement The author appreciates the valuable discussions with Gen Nakamura (Hokkaido University), Masahiro Yamamoto (The University of Tokyo) and Guanghui Hu (Beijing Computational Science Research Center).

References

[1] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975.
[2] I. Daubechies, M. Defrise, C. De Mol, An iterative thresholding algorithm for linear inverse problems, Commun. Pure Appl. Math., 57, 2004, 1413–1457.
[3] L.C. Evans, Partial Differential Equations, AMS, Providence, RI, 1998.
[4] R. Gorenflo, Y. Luchko, M. Yamamoto, Time-fractional diffusion equation in the fractional Sobolev spaces, Fract. Calc. Appl. Anal., 18, 2015, 799–820.
[5] G. Hu, Y. Kian, P. Li, Y. Zhao, Inverse moving source problems in electrodynamics, Inverse Problems, 35, 2019, 075001.
[6] G. Hu, Y. Liu, M. Yamamoto, Inverse moving source problem for fractional diffusion(-wave) equations: Determination of orbits, in: Inverse Problems and Related Topics, Springer Proceedings in Mathematics & Statistics 310, Springer, Singapore, 2020, pp. 81–100.
[7] D. Jiang, Z. Li, Y. Liu, M. Yamamoto, Weak unique continuation property and a related inverse source problem for time-fractional diffusion advection equations, Inverse Problems, 33, 2017, 055013.
[8] D. Jiang, Y. Liu, D. Wang, Numerical reconstruction of the spatial component in the source term of a time-fractional diffusion equation, Adv. Comput. Math., accepted, arXiv:1812.04235.
[9] D. Jiang, Y. Liu, M. Yamamoto, Inverse source problem for a wave equation with final observation data, in: Mathematical Analysis of Continuum Mechanics and Industrial Applications, Mathematics for Industry 26, Springer, Singapore, 2017, 153–164.
[10] D. Jiang, Y. Liu, M. Yamamoto, Inverse source problem for the hyperbolic equation with a time-dependent principal part, J. Differential Equations, 262, 2017, 653–681.
[11] B. Jin, R. Lazarov, Z, Zhou, Error estimates for a semidiscrete finite element method for fractional order parabolic equations, SIAM J. Numer. Anal., 51, 2013, 445–466.
[12] B. Jin, W. Rundell, A tutorial on inverse problems for anomalous diffusion processes, Inverse Problems, 31, 2015, 035003.
[13] A. Kubica, M. Yamamoto, Initial-boundary value problems for fractional diffusion equations with time-dependent coefficients, Fract. Calc. Appl. Anal., 21, 2018, 276–311.
[14] Z. Li, Y. Liu, Yamamoto M. Inverse problems of determining parameters of the fractional partial differential equations, in: Handbook of Fractional Calculus with Applications. Volume 2: Fractional Differential Equations, De Gruyter, Berlin, 2019, pp. 431–442.
[15] K. Liao, T. Wei, Identifying a fractional order and a space source term in a time-fractional diffusion-wave equation simultaneously, Inverse Problems, 35, 2019, 115002.
[16] Y. Lin, C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys., 225, 2007, 1533–1552.
[17] Y. Liu, G. Hu, M. Yamamoto, Inverse moving source problem for fractional diffusion(-wave) equations: Determination of profiles, submitted, arXiv:2002.00194.
[18] Y. Liu, D. Jiang, M. Yamamoto, Inverse source problem for a double hyperbolic equation describing the three-dimensional time cone model, SIAM J. Appl. Math., 75, 2015, 2610–2635.
[19] Li Z, Yamamoto M. Inverse problems of determining coefficients of the fractional partial differential equations, in: Handbook of Fractional Calculus with Applications. Volume 2: Fractional Differential Equations, De Gruyter, Berlin, 2019, pp. 443–464.
[20] Y. Liu, Z. Li, M. Yamamoto, Inverse problems of determining sources of the fractional partial differential equations, in: Handbook of Fractional Calculus with Applications. Volume 2: Fractional Differential Equations, De Gruyter, Berlin, 2019, pp. 431–442.
[21] Y. Liu, W. Rundell, M. Yamamoto, Strong maximum principle for fractional diffusion equations and an application to an inverse source problem, Fract. Calc. Appl. Anal., 19, 2016, 888–906.
[22] E. Nakaguchi, H. Inui, K. Ohnaka, An algebraic reconstruction of a moving point source for a scalar wave equation, Inverse Problems, 28, 2012, 065018.
[23] T. Ohe, Real-time reconstruction of moving point/dipole wave sources from boundary measurements, Inverse Probl. Sci. Eng., accepted.
[24] K. Sakamoto, M. Yamamoto, Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems, J. Math. Anal. Appl., 382, 2011, 426–447.
[25] Y.B. Wang, X.Z. Jia, J. Cheng, A numerical differentiation method and its application to reconstruction of discontinuity, Inverse Problems, 18, 2002, 1461–1476.

Research Institute for Electronic Science
Hokkaido University
N12W7, Kita-Ward, Sapporo 060-0812, JAPAN
E-mail address: [email protected]