On the identification of a nonlinear term in a reaction-diffusion equation

TL;DR

This paper investigates the inverse problem of recovering a nonlinear reaction term in a reaction-diffusion equation, including fractional models, using overposed data and fixed point and Newton scheme methods, supported by numerical reconstructions.

Contribution

It introduces a novel approach to recover the entire nonlinear reaction term in reaction-diffusion equations, including fractional models, using fixed point and Newton methods with numerical validation.

Findings

Established conditions for the fixed point map to be contractive.

Developed Newton schemes for the inverse problem.

Numerical reconstructions demonstrate the effectiveness of the methods.

Abstract

Reaction-diffusion equations are one of the most common partial differential equations used to model physical phenomenon. They arise as the combination of two physical processes: a driving force $f(u)$ that depends on the state variable $u$ and a diffusive mechanism that spreads this effect over a spatial domain. The canonical form is $u_t - \triangle u = f(u)$. Application areas include chemical processes, heat flow models and population dynamics. The direct or forwards problem for such equations is now very well-developed and understood. However, our interest lies in the inverse problem of recovering the reaction term $f(u)$ not just at the level of determining a few parameters in a known functional form, but recovering the complete functional form itself. To achieve this we set up the usual paradigm for the parabolic equation where $u$ is subject to both given initial and boundary data, then prescribe overposed data consisting of the solution at a later time $T$. For example, in the case of a population model this amounts to census data at a fixed time. Our approach will be two-fold.First we will transform the inverse problem into an equivalent nonlinear mapping from which we seek a fixed point. We will be able to prove important features of this map such as a self-mapping property and give conditions under which it is contractive. Second, we consider the direct map from $f$ through the partial differential operator to the overposed data. We will investigate Newton schemes for this case. In recent decades various anomalous processes have been used to generalize classical Brownian diffusion. Amongst the most popular is one that replaces the usual time derivative by a fractional one of order $\alpha \leq 1$. We will also include this model in our analysis. The final section of the paper shows numerical reconstructions that demonstrate the viability of the suggested approaches.

Equations

$$u_t(x,t)-\mathbb{L}u(x,t)=f(u(x,t))+r(x,t),(x,t)\in \Omega \times (0,T)$$

$${\partial }_{\nu }u(x,t)+\gamma (x)u(x,t)$$

$$u(x,0)$$

$$u(x,T)=g(x).$$

$$D_t^{\alpha }u(x,t)-\mathbb{L}u(x,t)=f(u(x,t))+r(x,t),(x,t)\in \Omega \times (0,T)$$

$${}_a{D}_t^{\alpha }u=\frac{1}{\Gamma (1-\alpha )}{\int }_a^t\frac{u^{\prime }(s)}{(t-s{)}^{\alpha }}ds$$

$$(Au)(x)=(-\mathbb{L}u)(x)x\in \Omega ,$$

$${\dot{H}}^s(\Omega )=\left\{v\in L^2(\Omega ):\sum _{j=1}^{\infty }{\lambda }_j^s|(v,{\phi }_j){|}^2<\infty \right\},$$

$$\Vert v{\Vert }_{{\dot{H}}^s(\Omega )}^2=\sum _{j=1}^{\infty }|(v,{\lambda }_j^{\frac{s}{2}}{\phi }_j){|}^2=\sum _{j=1}^{\infty }(v,A^{\frac{s}{2}}{\phi }_j{)}^2=\sum _{j=1}^{\infty }(A^{\frac{s}{2}}v,{\phi }_j{)}^2=\Vert A^{\frac{s}{2}}v{\Vert }_{L^2(\Omega )}^2.$$

$$D_t^{\alpha }u(x,t)-\mathbb{L}u(x,t)$$

$${\partial }_{\nu }u(x,t)+\gamma (x)u(x,t)$$

$$u(x,0)$$

$$u(x,t)$$

$$+\sum _{j=1}^{\infty }{\int }_0^t(t-s{)}^{\alpha -1}{E}_{\alpha ,\alpha }(-{\lambda }_j(t-s{)}^{\alpha })(r(x,s),{\phi }_j)ds{\phi }_j(x).$$

$$E_{\alpha ,\beta }(z)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\beta )}\alpha >0,\beta \in \mathbb{R},z\in \mathbb{C},$$

$$E_{\alpha ,\beta }(z)=\frac{1}{\alpha }{z}^{\frac{1-\beta }{\alpha }}{e}^{z^{\frac{1}{\alpha }}}-\sum _{k=1}^N\frac{1}{\Gamma (\beta -\alpha k)}\frac{1}{z^k}+O\left(\frac{1}{z^{N+1}}\right),$$

$$E_{\alpha ,\beta }(z)=-\sum _{k=1}^N\frac{1}{\Gamma (\beta -\alpha k)}\frac{1}{z^k}+O\left(\frac{1}{z^{N+1}}\right).$$

$$|E_{\alpha ,\beta }(z)|$$

$$u(t)=\mathbb{E}(t)u_0+{\int }_0^t\overline{\mathbb{E}}(t-s)r(\cdot ,s)ds,$$

$$\mathbb{E}(t)v=\sum _{j=1}^{\infty }{E}_{\alpha }(-{\lambda }_j{t}^{\alpha })(v,{\phi }_j){\phi }_j$$

$$\overline{\mathbb{E}}(t)v=\sum _{j=1}^{\infty }{t}^{\alpha -1}{E}_{\alpha ,\alpha }(-{\lambda }_j{t}^{\alpha })(v,{\phi }_j){\phi }_j.$$

$$\frac{d}{dt}\mathbb{E}(t)v=\bar{\mathbb{E}}(t)\mathbb{L}v.$$

$$\underset{t\to 0}{\lim }\Vert u(\cdot ,t)-u_0{\Vert }_{{\dot{H}}^{-\sigma }(\Omega )}=0$$

$$\Vert \mathbb{E}(t)v{\Vert }_{{\dot{H}}^p(\Omega )}\le c{t}^{\frac{q-p}{2}\alpha }\Vert v{\Vert }_{{\dot{H}}^q(\Omega )},$$

$$\Vert \mathbb{E}(t)v{\Vert }_{{\dot{H}}^p(\Omega )}^2$$

$$=t^{(q-p)\alpha }\sum _{j=1}^{\infty }{\lambda }_j^{p-q}{t}^{(p-q)\alpha }|E_{\alpha }(-{\lambda }_j{t}^{\alpha }){|}^2{\lambda }_j^q(v,{\phi }_j{)}^2$$

$$|E_{\alpha }(-{\lambda }_j{t}^{\alpha })|\le \frac{c}{1+{\lambda }_j{t}^{\alpha }},$$

$$\Vert \mathbb{E}(t)v{\Vert }_{{\dot{H}}^p(\Omega )}^2\le c{t}^{(q-p)\alpha }\underset{j}{\sup }\frac{{\lambda }_j^{p-q}{t}^{(p-q)\alpha }}{(1+{\lambda }_j{t}^{\alpha }{)}^2}\sum _{j=1}^{\infty }{\lambda }_j^q(v,{\phi }_j{)}^2\le c{t}^{(q-p)\alpha }\Vert v{\Vert }_{{\dot{H}}^q(\Omega )}^2,$$

$${\lambda }_j^{p-q}{t}^{(p-q)\alpha }\le \frac{p-q}{2}({\lambda }_j{t}^{\alpha }{)}^2+1-\frac{p-q}{2}\le (1+{\lambda }_j{t}^{\alpha }{)}^2.$$

$$\Vert \mathbb{E}(t)v{\Vert }_{{\dot{H}}^p(\Omega )}\le c{t}^{(q-p)/2}\Vert v{\Vert }_{{\dot{H}}^q(\Omega )},$$

$$\Vert g\ast h{\Vert }_{L^p(0,T)}\le \Vert g{\Vert }_{L^q(0,T)}\Vert h{\Vert }_{L^r(0,T)},$$

$$\Vert \bar{\mathbb{E}}(t)v{\Vert }_{{\dot{H}}^p(\Omega )}\le c{t}^{-1+(1+\frac{q-p}{2})\alpha }\Vert v{\Vert }_{{\dot{H}}^q(\Omega )},$$

Full text

However, our interest lies in the inverse problem of recovering the reaction term $f(u)$ not just at the level of determining a few parameters in a known functional form, but recovering the complete functional form itself. To achieve this we set up the usual paradigm for the parabolic equation where $u$ is subject to both given initial and boundary data, then prescribe overposed data consisting of the solution at a later time $T$. For example, in the case of a population model this amounts to census data at a fixed time. Our approach will be two-fold. First we will transform the inverse problem into an equivalent nonlinear mapping from which we seek a fixed point. We will be able to prove important features of this map such as a self-mapping property and give conditions under which it is contractive. Second, we consider the direct map from $f$ through the partial differential operator to the overposed data. We will investigate Newton schemes for this case.

Classical, Brownian motion diffusion is not the only version and in recent decades various anomalous processes have been used to generalize this case. Amongst the most popular is one that replaces the usual time derivative by a subdiffusion process based on a fractional derivative of order $\alpha\leq 1$. We will also include this model in our analysis. The final section of the paper will show numerical reconstructions that demonstrate the viability of the suggested approaches. This will also include the dependence of the inverse problem on both $T$ and $\alpha$.

1 Introduction

Before the thorough mathematical theory of existence, uniqueness and regularity results for reaction-diffusion equations initiated during the 1960’s, there had been an extensive literature in applications dating from the 1930’s. These include the work of Fisher in considering the Verhulst logistic equation together with a spatial diffusion, $u_{t}-ku_{xx}=au(1-u)$ and that of Kolmogorov, Petrovskii and Piskunov in similar models now collectively referred to as the Fisher-KPP theory; more complex polynomial terms occurring in combustion theory due to Zeldovich and Frank-Kamenetskii.

The class of reaction-diffusion equations describe the behaviour of a wide range of chemical systems where diffusion of material competes with production by a chemical reaction. They also describe systems where heat is produced depending on the local temperature and diffuses away by some mechanism. These systems can also be based on particles as in the nuclear diffusion equation.

For general references here we note the books [CantrellCosner:2003, Grindrod:1996, Murray:2002] as well as the chapter [Kuttler17]. These equations and their generalization can be described as follows.

Let $\Omega$ be a bounded, simply-connected region in $\mathbb{R}^{d}$ with smooth boundary $\partial\Omega$, and let $\mathbb{L}$ be a uniformly elliptic operator of second order with with $L^{\infty}$ coefficients.

[TABLE]

and subject to the initial and boundary conditions

[TABLE]

In (1) the usual, direct problem is to know both the diffusion operator $L$, the nonhomogeneous linear term $r$ and the reaction term $f(u)$. However, even in the earliest applications it was appreciated that the diffusion term $-k\triangle u$ might have the diffusion coefficient $k$ as an unknown and while the specific form of $f(u)$ was assumed known from the underlying model, specific parameters may not be. In our case we assume that $f$ is unknown (save that it be a function of $u$ alone). Thus we are interested in the inverse problem of recovering $f$ from additional data that would be overposed if in fact $f$ were known. There has been previous work along such lines. For example, in a series of papers, [PilantRundell:1986, PilantRundell:1987, PilantRundell:1988], Pilant and Rundell investigated the recovery of $f$ from overposed data consisting of a time trace of the solution $u$ at a given point $x_{0}\in\partial\Omega$, namely $u(x_{0},t)$ for all $t>0$. See also, [DuChateauRundell:1985]. We will take a very different situation here by assuming that one is able to measure in the spatial direction by taking a snapshot at a fixed later time $T$ and so our overposed data will be

[TABLE]

We also note that in the case of such data, if the reaction term is of the form $q(x)f(u)$ where $q$ is unknown but the actual form of the nonlinearity $f(u)$ is known, then it is possible to recover the spatial component in a unique way, [KaltenbacherRundell:2019b].

Exchanging the time derivative in (1) for one of fractional order adds further complexity but also of considerable physical interest as we are modifying in a fundamental way the diffusion process (Fick’s law) to one of a so-called anomalous diffusion mechanism that would lead to a non-local, time-fractional equation.

[TABLE]

and subject to the same initial and boundary conditions (2). Here we have also included the possibility of an additional (known) inhomogeneity, which can be made use of to drive the evolution process.

By $D_{t}^{\alpha}$ we mean the Djrbashian-Caputo derivative of order $\alpha$ which is defined to be

[TABLE]

and since this is a non-local operator, the designation of the starting point $t=a$ is essential; in our case we shall take $a=0$. Of course there are other possibilities for anomalous diffusion models, but the one chosen is both standard in most applications and has an existing rich mathematical theory.

Examples of the above include applications in many of the areas in which traditional reaction diffusion equations have featured over the last one hundred years and form a natural generalisation. The central feature of a reaction-diffusion model is the interplay between the two central processes that are in some sense in competition. One of the main themes of this paper is to investigate the role the different diffusion models play in our ability to recover the reaction term $f(u)$.

Briefly, the outline of the paper is as follows. In the next section we collect crucial main existence and regularity results for the direct problem. These are well-known for the parabolic case but we require an addition to the existing literature for the case of (4). Sections 3 and 4 are devoted to two iterative approaches for numerically solving the inverse problem, a fixed point scheme and a Newton type method. While these methods are formulated in the general, possibly higher dimensional setting, the analysis will be restricted to the spatially one-dimensional case, due to the need for certain Sobolev embeddings. In the last section we show some numerical examples illustrating the above and use this to bring out certain features of the problem. In particular we investigate the difference between the classical and fractional diffusion operators and the role that the value of the final time $T$ plays.

2 The direct problem

In this section we provide the results and tools for the direct problem, i.e., solution of the initial boundary value problem (2), (4), that will be needed later on in the reconstruction methods. These are, first of all, a representation of the solution to the inhomogeneous linear problem by means of Mittag-Leffler functions as well as eigenvalues and eigenfunctions of the operator $-\mathbb{L}$. Moreover, versions of the familiar Gronwall’s inequality for ordinary differential equations in function space are required for the fractional case and we develop these below, first for the linear case, then later for the semilinear. These will then be used to provide estimates needed for the semilinear fractional operator.

2.1 The Solution Representation and the Mittag-Leffler function

Throughout this paper we will assume that

$\mathbb{L}=\sum_{ij}\partial_{x_{i}}(a_{ij}\partial_{x_{j}})+c$

is a symmetric uniformly elliptic operator with coefficients $a_{ij}\in W^{1,\infty}(\Omega)$, $b_{i},c\in L^{\infty}(\Omega)$, and that $\Omega$ is a bounded and convex or $C^{1,1}$ domain.

We define an operator $A$ in $L^{2}(\Omega)$ by

[TABLE]

where $\mathbb{L}$ is equipped with impedance boundary conditions, with its domain $D(A)\subseteq H^{2}(\Omega)$. Then the fractional power $A^{\sigma}$ is defined for any $\sigma\in\mathbb{R}$. Since $\mathbb{L}$ is a symmetric uniformly elliptic operator, the spectrum of $A$ is entirely composed of eigenvalues and counting according to the multiplicities, we can set $0<\lambda_{1}\leq\lambda_{2}\ldots$. By $\varphi_{j}\in H^{2}(\Omega)$, we denote the $L^{2}(\Omega)$ orthonormal eigenfunctions corresponding to $\lambda_{j}$. Next we introduce a space $\dot{H}^{s}(\Omega)$ by

[TABLE]

which is a Hilbert space with the norm $\|v\|_{\dot{H}^{s}(\Omega)}^{2}=\sum_{j=1}^{\infty}\lambda_{j}^{s}|(v,\varphi_{j})|^{2}$. By definition, we have the following equivalent form:

[TABLE]

We have $\dot{H}^{s}(\Omega)\subseteq H^{s}(\Omega)$ for $s\in[0,2]$. Since $\dot{H}^{s}(\Omega)\subseteq L^{2}(\Omega)$, identifying the dual $(L^{2}(\Omega))^{\prime}$ with itself, we have $\dot{H}^{s}(\Omega)\subseteq L^{2}(\Omega)\subseteq(\dot{H}^{s}(\Omega))^{\prime}$. Henceforth, we set $\dot{H}^{-s}(\Omega)=(\dot{H}^{s}(\Omega))^{\prime}$, which consists of bounded linear functionals on $\dot{H}^{s}(\Omega)$.

The solution $u$ to the linear problem

[TABLE]

can be represented by, [JinRundell:2015]

[TABLE]

Here $E_{\alpha,\beta}(z)$ is the Mittag-Leffler function [Mittag-Leffler:1903] defined by

[TABLE]

This is an entire function of order $1/\alpha$ and type 1. For the general properties of this function we suggest the references [MainardiGorenflo:2000, SamkoKilbasMarichev:1993]. The definitive references are due to Djrbashian, [Dzharbashyan:1964, Djrbashian:1970, Djrbashian:1993] although not so immediately accessible. See also the tutorial on inverse problems for fractional operators, [JinRundell:2015]. In the case $\alpha=1$, $E_{\alpha,1}(z)=E_{\alpha,\alpha}(z)=e^{z}$ and the model (7) as well its solution formula (8) recovers the standard parabolic equation.

For our purposes, in fact for a large range of inverse problems involving a subdiffusion process, one of the key features of the Mittag-Leffler function is its decay for large, real and negative argument. The asymptotic behaviour as $z\to\infty$ in various sectors of the complex plane $\mathbb{C}$ was first derived by Djrbashian [Dzharbashyan:1964], and subsequently refined by many over the intervening years.

Lemma 2.1**.**

Let $\alpha\in(0,2)$, $\beta\in\mathbb{R}$, and $\mu\in(\alpha\pi/2,\min(\pi,\alpha\pi))$, and $N\in\mathbb{N}$. Then for $|\mathrm{arg}(z)|\leq\mu$ with $|z|\to\infty$, we obtain

[TABLE]

and for $\mu\leq|\mathrm{arg}(z)|\leq\pi$ with $|z|\to\infty$

[TABLE]

In particular, for $x$ real and positive, $E_{\alpha,\beta}(-x)$ decays to zero linearly as $x\to\infty$.

The following useful estimate is a direct corollary of Lemma 2.1.

Corollary 2.1**.**

Let $0<\alpha<2$, $\beta\in\mathbb{R}$, and $\pi\alpha/2<\mu<\min(\pi,\pi\alpha)$. Then the following estimates hold

[TABLE]

The solution representation (8) can be succinctly rewritten as

[TABLE]

where the solution operators $\mathbb{E}$ and $\overline{\mathbb{E}}$ are defined by, respectively

[TABLE]

and

[TABLE]

Since $\frac{d}{dt}E_{\alpha,1}(-\lambda_{j}t^{\alpha})=-\lambda_{j}t^{\alpha-1}E_{\alpha,\alpha}(-\lambda_{j}t^{\alpha})$ we have

[TABLE]

2.2 Existence, Uniqueness and Regularity

The differences between the parabolic ($\alpha=1$) and fractional ($\alpha<1$) cases is substantial in terms of the regularity of the direct solution even in the homogeneous case $r=0$. For the parabolic problem, if $u_{0}\in L^{2}(\Omega)$ then $u(t)$ is infinitely differentiable in $t$ and so we have an infinite pick up in regularity of the solution. For the fractional problem, this is not the case; one can only recover at most two weak derivatives, see [SakamotoYamamoto:2011a]. This difference is due to the exponential decay of $e^{-\lambda_{j}t}$, in the parabolic case but the only linear decay of $E_{\alpha,1}(-\lambda_{j}t^{\alpha})$ for $0<\alpha<1$.

There is also a distinction for the case of a linear nonhomogeneous term as shown in [SakamotoYamamoto:2011a] and Lemma 2.3 below which we will use to obtain a corresponding result for the nonlinear situation (4). As in the parabolic case, the main tool for this is Gronwall’s inequality and we will need a version of this for the fractional operator case.

Now we introduce the concept of weak solution for problem (7).

Definition 2.1**.**

We call $u$ a solution to problem (7) if (7) holds in $L^{2}(\Omega)$ and $u(\cdot,t)\in H_{0}^{1}(\Omega)$ for almost all $t\in(0,T)$ and $u\in C([0,T];\dot{H}^{-\sigma}(\Omega))$, with

[TABLE]

with some $\sigma\geq 0$, which may depend on $\alpha$.

In the homogeneous case, existence and uniqueness of a solution to (7) can be concluded from (12) and the following stability estimate on the solution operator $\mathbb{E}(t)$ defined in (13).

Lemma 2.2**.**

For the operator $\mathbb{E}(t)$ defined in (13), we have, for any $t>0$,

[TABLE]

where $0\leq q\leq 2$ and $q\leq p\leq q+2$, and the constant $c$ depends only on $\alpha$ and $p-q$.

Proof.

By the definition of the operator $\mathbb{E}(t)$, there holds

[TABLE]

By Corollary 2.1, we have

[TABLE]

where the constant $c$ depends only on $\alpha$. Consequently

[TABLE]

where the last inequality follows from the choice of $p$ and $q$ such that $0\leq p-q\leq 2$ and Young’s inequality $ab\leq\frac{1}{P}a^{P}+\frac{P-1}{P}b^{\frac{P}{P-1}}$, which yields

[TABLE]

∎

Remark 2.1**.**

The choice $p\leq q+2$ in Lemma 2.2 indicates that the operator $\mathbb{E}(t)$ has at best a smoothing property of order two in space, which contrasts sharply with the classical parabolic case, for which the following estimate holds [Thomee:2006]

[TABLE]

for $q\geq 0$, and any $p\geq q$. This is reflected in the exponential decay of $e^{-\lambda_{j}t}$, instead of the linear decay $E_{\alpha,1}(-\lambda_{j}t^{\alpha})$ for $0<\alpha<1$.

Using Young’s inequality for convolution

[TABLE]

for $g\in L^{q}(0,T)$, $h\in L^{r}(0,T)$, $p^{-1}+1=q^{-1}+r^{-1}$, $p,q,r\geq 1$, one can similarly prove the following result on the the solution operator $\bar{\mathbb{E}}(t)$.

Lemma 2.3**.**

For the operator $\bar{\mathbb{E}}(t)$ defined in (14), there holds, for any $t>0$,

[TABLE]

where $0\leq q\leq 2$ and $q\leq p\leq q+2$, and the constant $c$ depends only on $\alpha$ and $p-q$.

Gronwall’s inequalities are a basic tool for differential equations; here is the classical integral formulation that does not require interior differentiability of the function $u(t)$.

Theorem 2.1**.**

Let $\beta$ and $u$ be real-valued continuous functions defined on the interval $[a,b]$ while $\alpha$ is integrable on every subinterval of $(a,b)$. If $\beta(t)\geq 0$ and if

[TABLE]

then

[TABLE]

If, in addition, the function $\alpha$ is non-decreasing, then

[TABLE]

We will need the corresponding results for fractional order equations

Theorem 2.2**.**

Let $\beta>0$. Suppose that $a(t)$ is a nonnegative, locally integrable function on $(0,T)$ for some $T>0$, while $g(t)$ is a nonnegative, nondecreasing continuous function on $(0,T)$. If $u(t)$ is nonnegative and locally integrable on $(0,T)$ such that

[TABLE]

Then

[TABLE]

In addition, if $a(t)$ is nondecreasing on $[0,T)$ then

[TABLE]

Proof.

For locally integrable functions $\phi(t)$ on $[0,T)$ define the operator $B$ by $B\phi(t)=g(t)\int_{0}^{t}(t-\tau)^{\beta-1}\phi(\tau)\,d\tau$. Then for each integer $n$ by successively applying (19) we obtain

[TABLE]

We claim that

[TABLE]

The proof of this is by induction and the inequality (23) is clearly true for $n=1$. Assume now it is true for $n=j$, then we have

[TABLE]

where the second line follows from the first from the fact that $g(t)$ is nondecreasing and the third line follows by an obvious change of variables. Writing the integrand as a Beta function and then invoking the Beta-Gamma function duplication formula, yields the identity $\int_{s}^{t}(t-\tau)^{\beta-1}(\tau-s)^{j\beta-1}\,d\tau=\frac{\Gamma(j\beta)\Gamma(\beta)}{\Gamma((j+1)\beta)}(t-s)^{(j+1)\beta-1}$. This now shows (23).

With $M=\sup_{[0,T]}g(t)$ we obtain

[TABLE]

and it is easily seen that

[TABLE]

From the above we now obtain (20).

If $a(t)$ is nondecreasing on $[0,T)$ then we can write (20) as

[TABLE]

where we have used the fact that $\int_{0}^{t}(t-\tau)^{n\beta-1}\,d\tau=\int_{0}^{t}s^{n\beta-1}=\frac{\Gamma(n\beta)}{\Gamma(n\beta+1)}t^{n\beta}$. ∎

An important particular case occurs when $g$ is constant and the expression on the right hand side of the estimate can be written as a Mittag-Leffler function.

Lemma 2.4**.**

Suppose $b\geq 0$, $\beta>0$ and $a(t)$ is a nonnegative function locally integrable on $0\leq t<T$ (for some $T<+\infty$), and suppose $u(t)$ is nonnegative and locally integrable on $0\leq t<T$ with

[TABLE]

Then

[TABLE]

Moreover, the following convolution version of Gronwall’s inequality will be employed.

Lemma 2.5**.**

Suppose $c\geq 0$ and $b(t)$ is a nonnegative function locally integrable on $0\leq t<T$ (for some $T<+\infty$), and suppose $u(t)$ is nonnegative and locally integrable on $0\leq t<T$ with

[TABLE]

Then

[TABLE]

Proof.

The proof proceeds as in that of Theorem 2.2 but in this case it is easier since there is no singularity involved. Define $(Ba)(t)=c\int_{0}^{t}e^{-\lambda(t-s)}a(s)\,ds$, then the inequality for $u$ reads $u(t)\leq b(t)+(Bu)(t)$ and implies $u(t)\leq\sum_{k=0}^{N}(B^{k}b)(t)+(B^{N}u)(t)$, where (by induction) $(B^{k}a)(t)=c^{k}\int_{0}^{t}\frac{(t-s)^{k-1}}{(k-1)!}e^{-\lambda(t-s)}a(s)\,ds\,\to\,0$ as $k\to\infty$, hence

[TABLE]

∎

2.3 Semilinear problems

Now we briefly discuss the slightly more complex case, i.e., semilinear problems, using a fixed point argument.

Consider the following initial boundary value problem for a semilinear subdiffusion equation

[TABLE]

In the model, we assume

[TABLE]

The argument below is based on the operator theoretic approach in $L^{2}(\Omega)$, see also [JinLiZhou17, Theorem 3.1].

Theorem 2.3**.**

Let $u_{0}\in\dot{H}^{2}(\Omega)$, $f(u_{0})\in L^{2}(\Omega)$, and (29) hold. Then the solution $u$ to problem (28) belongs to $C([0,T];\dot{H}^{2}(\Omega))\cap C^{1}((0,T];L^{2}(\Omega))$ such that $D_{t}^{\alpha}u\in C((0,T];L^{2}(\Omega))$.

Proof.

For the unique existence of $u$, it suffices to discuss the integral equation cf. (12)

[TABLE]

First we estimate $u_{t}(t)$, which is given by

[TABLE]

where we have used (15). Consequently, for any $0<t\leq T$, by Lemma 2.3 with $p=q=0$, there holds

[TABLE]

Using Assumption (29) on $f$, we have for $0<t\leq T$

[TABLE]

By Lemma 2.4, we obtain

[TABLE]

and $u\in C^{1}((0,T];L^{2}(\Omega))$ and thus

[TABLE]

i.e., $D_{t}^{\alpha}u\in C([0,T];L^{2}(\Omega))$. Further, using again (15), for any $\sigma\in[0,2]$

[TABLE]

Clearly, $\|\mathrm{I}_{1}\|_{L^{2}(\Omega)}\leq c\|u_{0}\|_{\dot{H}^{\sigma}(\Omega)}$. By Lemma 2.2, for $\mathrm{I}_{2}$ and $\mathrm{I}_{3}$, we have

[TABLE]

By the mean value theorem, we have

[TABLE]

Consequently,

[TABLE]

Hence, by Lemma 2.3 with $p=\sigma$, $q=0$, we deduce

[TABLE]

Likewise, since $r(s)-r(t)=\int_{t}^{s}r_{t}(\tau)d\tau$,

[TABLE]

Therefore, we have $A^{\frac{\sigma}{2}}u\in C([0,T];L^{2}(\Omega))$, i.e., $u\in C([0,T];\dot{H}^{\sigma}(\Omega))$. This completes the proof of the theorem. ∎

3 An iteration scheme to recover $f$

A natural scheme to recover $f$ is to evaluate the equation (4) on the overposed boundary $t=T$. That is, we define a map $\mathbb{T}:f\to u(x,T;f)$ by

[TABLE]

where $g$ is the given data and define a sequence of approximations by the fixed point iteration $f_{k+1}=\mathbb{T}(f_{k})$. We start with the idealized situation of exact data and later on in Section 3.5 consider the realistic setting of noisy data.

Before we can utilize the map (30) we must obtain conditions on the data that guarantee it is well defined. Specifically, the range of $g(x)$ must contain all values of the solution $u(x,t;f)$ for $t\leq T$:

[TABLE]

The operator $\mathbb{T}$ is a concatenation $\mathbb{T}=\mathbb{P}\circ\mathbb{S}$ of the operators $\mathbb{P}:Z\to X$ and $\mathbb{S}:X\to Z$ defined by

[TABLE]

where we choose $\rho_{0}\geq\|f_{act}(u_{0})+r(\cdot,T){-\mathbb{L}u_{0}}\|_{\dot{H}^{\sigma}(\Omega)}$ to make sure that a fixed point of $\mathbb{T}$ solves the original inverse problem (35), (3).

We use a bounded interval $I=[g_{\min},g_{\max}]$ with

$g_{\min}=\min\{g(x):x\in\Omega\}$, $g_{\max}=\max\{g(x):x\in\Omega\}$,

in order to be able to make use of embedding theorems. Then we use the function space setting

[TABLE]

with $\sigma$ such that ${\dot{H}^{\sigma}(\Omega)}$ continuously embeds into $W^{1,\infty}(\Omega)$ and with the norm

[TABLE]

Moreover, we employ the projection $P(z)=\max\{\min\{z,g_{\max}\},g_{\min}\}$ on $I$ to define $u(x,t;f)$ as solution to

[TABLE]

Hence, based on the range condition (31), which we need to assume to hold for the actual nonlinearity $f_{act}$ only, and the fact that $u(x,T;f_{act})$ solves (4), we can replace the original model (4) by the equation containing the projection (35).

In order to establish well-definedness of the operator $\mathbb{P}$ we will assume throughout this section that

[TABLE]

for some $0<\underline{g}<\overline{g}$, so that there exist $c(g),C(g)>0$ such that

[TABLE]

Indeed, we have in case $\sigma=2$

[TABLE]

as well as in case $\sigma=1$

[TABLE]

from which the general case $q\in[1,2]$ follows by interpolation. Moreover,

[TABLE]

Lemma 3.1**.**

Under condition (36), the operator $\mathbb{P}:Z\to X$ is well-defined by (32) between the spaces defined by (33), and satisfies

[TABLE]

Proof.

Existence of a minimizer follows from the fact that

$\{\|f^{\sharp}(g)-y\|_{Z}\,:\,f^{\sharp}\in X\mbox{ and }\|f(u_{0})\|_{\dot{H}^{\sigma}(\Omega)}\leq\rho_{0}\}$ is a nonempty closed convex (hence weakly* closed) subset of the space $X$ which is the dual of a separable space. Moreover, from (37) and the triangle inequality as well as by minimality of $\mathbb{P}f$, comparing with $f^{\sharp}=0$, we get

[TABLE]

∎

Remark 3.1**.**

For $f\in X$, $r(\cdot,T)\in Z$, under the attainability condition

[TABLE]

any minimizer $f^{+}$ of (32) satisfies $f^{+}(g)=y$ and is therefore unique, due to the fact that $g(\Omega)=I$. The second condition in (40) can be shown to be satisfied by the strict monotonicity of $g$ according to (36), if $\Omega=(0,L)\subseteq\mathbb{R}^{1}$ and $u_{0}\in W^{1,\infty}(\Omega)\cap H^{\sigma}(\Omega)$ is strictly monotone. This can be seen by using the identity $f(u(\cdot,T;f))+\mathbb{L}(u(\cdot,T;f)-g)=D_{t}^{\alpha}u(\cdot,T;f)-r(\cdot,T)-\mathbb{L}g\in Z$ (this regularity of $D_{t}^{\alpha}u(\cdot,T;f)$ will be shown at the beginning of the next section) and using the assumed regularity (36) of $g$, which transfers to its inverse and gives $f^{+}=y\circ g^{-1}\in W^{1,\infty}(I)\cap H^{\sigma}(I)\subseteq X$, as an easy verification of the chain rule in $W^{1,\infty}(a,b)\cap H^{\sigma}(a,b)$ for $(a,b)=\Omega$ and $(a,b)=I$ shows.

In our implementation we use a regularized version of $\mathbb{P}$ by restricting the search for a minimizer of $\|f^{\sharp}(g)-y\|_{X}$ to a finite dimensional subset of smooth functions in $X$. This is actually an intrinsic part of the algorithm and to some extent explains its good contraction performance in spite of the fact that in an infinite dimensional function space setting, contractivity of $\,\mathbb{T}$ can only be proven in the conditional sense of Theorem 3.2 below.

Well-definedness of the operator $\mathbb{S}:X\to L^{2}(\Omega)$ follows from Theorem 2.3. The higher regularity $\mathbb{S}f\in{\dot{H}^{\sigma}(\Omega)}$ will be established in the following section.

3.1 Self-mapping property of the fixed point operator

To prove that $D_{t}^{\alpha}u(x,T;f)\in Z$ we make use of the the representation formulae (12), (13), (14)

[TABLE]

hence

[TABLE]

which by the triangle inequality and Young’s inequality for convolution (16) (with $p=\infty$, $q=Q$, $r=Q^{*}=\frac{Q}{Q-1}$) yields

[TABLE]

for $Q\leq\frac{2(2-\theta)}{2-2\theta+\sigma}$, that is, $Q^{*}=\frac{Q}{Q-1}\geq\frac{2(2-\theta)}{2-\sigma}$.

Here we have used the fact that

[TABLE]

for $\theta\in(0,2-1/\alpha)$

$\alpha\in(\frac{1}{2},1]$

and

[TABLE]

which leads to the requirement $\sigma\leq\frac{2}{Q}(2-Q+\theta(Q-1))=2-2\frac{Q-1}{Q}(2-\theta)$, (that is, $Q\leq\frac{2(2-\theta)}{2-2\theta+\sigma}$,) in order to cancel the powers of $\lambda_{j}$ in the term containing $f^{\prime}(Pu)$, so that no space derivatives are finally applied to $f^{\prime}(Pu)$.

We can therefore further estimate, noting that by the range condition (31) $u_{0}=Pu_{0}$

[TABLE]

If $Q^{*}<\infty$ then for $\eta(t)=\|u_{t}(t)\|_{\dot{H}^{\sigma}(\Omega)}^{Q^{*}}$ this means

[TABLE]

for $C_{0}=\|f(u_{0})+r(\cdot,0)-\mathbb{L}u_{0}\|_{{\dot{H}^{\sigma}(\Omega)}}$, $C_{1}=\frac{\tilde{C}_{\alpha}^{2/Q^{*}}}{\lambda_{1}^{\sigma}}\|f^{\prime}\|_{L^{\infty}(I)}$, $C_{r}(t)=\tilde{C}_{\alpha}^{2/Q^{*}}\|r_{t}\|_{L^{Q^{*}}(0,t;L^{2}(\Omega))}$ $\tilde{C}_{i}=2^{Q^{*}-1}C_{i}^{Q^{*}}$ for $i=0,1,r$.

Gronwall’s inequality (17) yields

[TABLE]

provided $\alpha$ is sufficiently close to one: $\alpha>1-\frac{1}{Q^{*}}=\frac{1}{Q}$. That is,

[TABLE]

(note that $C_{0}$, $\tilde{C}_{0}$ contain $\|f(u_{0})+r(\cdot,0)-\mathbb{L}u_{0}\|_{{\dot{H}^{\sigma}(\Omega)}}$ and $C_{1}$, $\tilde{C}_{1}$ contain $\|f^{\prime}\|_{L^{\infty}(I)}$)

Note that $Q^{*}=\infty$ does not work in the case $\alpha<1$, because of the singularity at $t=0$.

Also note that the condition on $\alpha$ arising here is compatible with the previous one $Q\leq\frac{2(2-\theta)}{2-2\theta+\sigma}$with $\theta\in(0,2-1/\alpha)$, since both together mean $\alpha>\frac{1}{Q}\geq\frac{2-2\theta+\sigma}{2(2-\theta)}=1-\frac{2-\sigma}{2(2-\theta)}>1-\frac{2-\sigma}{2}\alpha$. and can therefore be satisfied provided $\sigma<2$ and $\alpha$ is sufficiently close to one: $\alpha>1-\frac{2-\sigma}{4-\sigma}$. That is, since we need $\sigma>\frac{3}{2}$ to guarantee continuity of the embedding $\dot{H}^{\sigma}(\Omega)\to W^{1,\infty}(\Omega)$, this implies $\alpha>\frac{4}{5}$.

The function $\Phi$ appearing in the estimate (43) is an increasing function of the initial data $\|f(u_{0})+r(\cdot,0)-\mathbb{L}u_{0}\|_{{\dot{H}^{\sigma}(\Omega)}}$, of $\|f^{\prime}\|_{L^{\infty}(I)}$ and of $\|r_{t}\|_{L^{Q^{*}}(0,T;L^{2}(\Omega))}$; moreover for fixed $T$ and bounded $\|f^{\prime}\|_{L^{\infty}(I)}$, the value of $\Phi$ can be made arbitrarily small by making $\|f(u_{0})+r(\cdot,0)-\mathbb{L}u_{0}\|_{{\dot{H}^{\sigma}(\Omega)}}$ and $\|r_{t}\|_{L^{Q^{*}}(0,T;L^{2}(\Omega))}$ small.

We thus define

[TABLE]

and apply the above estimate (43) to get, for any $f\in B$,

[TABLE]

where the second inequality in (45) can be achieved by choosing $\rho$ sufficiently large $\rho>C(g)C_{\dot{H}^{\sigma},W^{1,\infty}}^{\Omega}\|\mathbb{L}g\|_{Z}$, and then choosing $\rho_{0}$, $\|r_{t}\|_{L^{Q^{*}}(0,t;L^{2}(\Omega))}$ small enough so that \Phi\bigl{(}\rho_{0},T,\rho,\|r_{t}\|_{L^{Q^{*}}(0,t;L^{2}(\Omega))}\bigl{)}\leq\frac{1}{C(g)C_{\dot{H}^{\sigma},W^{1,\infty}}^{\Omega}}\rho-\|\mathbb{L}g\|_{Z}. Moreover, by the constraint in the definition (32) of $\mathbb{P}$, we have

[TABLE]

i.e., altogether $\mathbb{T}(f)\in B$, provided $f\in B$.

3.2 Weak* compactness of $B$ and weak* continuity of $\mathbb{T}$

The set $B$ is by definition weakly* compact and convex in the Banach space $X$ with norm $\|f\|_{X}=\|f\|_{W^{1,\infty}(I)}+\|f(u_{0})\|_{{\dot{H}^{\sigma}(\Omega)}}$. For any sequence $(f_{n})_{n\in\mathbb{N}}\in B$ converging weakly* in $X$ to $f\in X$ we have that the sequence of images under $\mathbb{T}$, i.e., $\mathbb{T}(f_{n})$, is contained in the weakly* compact set $B$. Using this and compactness of the embeddings $X\to L^{\infty}(I)$ (where boundedness of the interval $I$ is crucial) we can extract a subsequence with indices $(n_{k})_{k\in\mathbb{N}}$ and an element $f^{+}\in B$ such that

[TABLE]

It remains to prove that $f^{+}=\mathbb{T}(f)$. For this purpose, we use the fact that the difference $\hat{u}_{n}:=u(x,t;f_{n})-u(x,t;f)$ of solution to (35) solves

[TABLE]

with homogeneous initial and impedance boundary conditions, where

$q_{n}(x,t)=\int_{0}^{1}f_{n}^{\prime}(P(u(x,t;f)+\sigma\hat{u}_{n}(x,t)))\,d\sigma P$, $\hat{f}_{n}=f_{n}-f$, and $u=u(x,t;f)$. From the representation (12), (13), (14)

[TABLE]

Young’s inequality (16) and (41), as well as the fact that $f_{n}\in B$, allows us to obtain the crude estimate

[TABLE]

which by Gronwall’s inequality (17) yields

[TABLE]

Using the fact that $\hat{u}_{n}^{j}$ as defined in (49) satisfies the fractional ODE $D_{t}^{\alpha}\hat{u}_{n}^{j}(t)+\lambda_{j}\hat{u}_{n}^{j}(t)=((q_{n}\hat{u}_{n}+\hat{f}_{n}(Pu))(\cdot,s),\varphi_{j})$ we obtain

[TABLE]

From Young’s inequality (16) and (41) we therefore obtain an estimate of $D_{t}^{\alpha}\hat{u}_{k}$ in a rather weak norm

[TABLE]

Thus,

under the attainability condition (40)

[TABLE]

by (47). On the other hand, (47) also implies

[TABLE]

Hence, the two limits need to coincide $\mathbb{T}(f)(g)=f^{+}(g)\in Z\subseteq C(\Omega)$, thus $\mathbb{T}(f)=f^{+}$. A subsequence-subsequence argument yields weak* convergence in $X$ of the whole sequence $\mathbb{T}(f_{n})$ to $\mathbb{T}(f)$.

Invoking Schauder’s Fixed Point Theorem in locally convex topological spaces; see [Fan1952], we have proven the following result.

Theorem 3.1**.**

*Let $\alpha\in(\frac{4}{5},1]$, $\Omega\subseteq\mathbb{R}^{1}$,

an open bounded interval

$\sigma\in(\frac{3}{2},2]$ and assume that $g$

is strictly monotone and

satisfies (36) with (33), that $\rho_{0}$ as well as $\|r_{t}\|_{L^{Q^{*}}(0,T;L^{2}(\Omega))}$ are sufficiently small.

Then for large enough $\rho>0$ the operator $\mathbb{T}$ is a self-mapping on the bounded, closed and convex set $B$ as defined in (44) and $\mathbb{T}$ is weakly* continuous in $X$ as defined in (33). Hence $\mathbb{T}f$ has a fixed point $f\in B$. If this fixed point $f$

is monotonically decreasing and

satisfies (31), then $f$ solves the inverse problem (2), (3), (4), .*

Proof.

It only remains to prove that under conditions (31) and (40), a fixed point $f$ solves (2), (3), (4). By Remark (3.1) we have $f(g)=y=f(Pu(\cdot,T;f)+\mathbb{L}(u(\cdot,T;f)-g)$, hence the difference $w=u(\cdot,T;f)-g$ satisfies the elliptic PDE

[TABLE]

with $\bar{y}=\int_{0}^{1}f^{\prime}(g+\theta Pw))P\,d\theta\leq 0$ and homogeneous boundary conditions, and therefore has to vanish. ∎

Below we will make a series of short remarks.

The question arises whether a similar result can achieved for $\alpha<\frac{4}{5}$. we have not investigated all possibilities to prove this, but also have no evidence that it cannot be done.

The assumption of strict monotonicity of $g$ is used to guarantee attainability (40), cf. Remark (3.1). Our iteration algorithm can be implemented without it and in practice our scheme worked perfectly well also for non-monotone $g$.

Estimate (43) shows that the influence of the initial data indeed decreases for larger $T$. This confirms the computational observation of faster convergence for larger $T$.

3.3 Contractivity for monotone nonlinearities in the parabolic case – uniqueness

To establish a contraction property of $\mathbb{T}$ in the parabolic case, we will make use of the fact that the solution $u$ to the forward problem (4), (2) and even its time derivative decays exponentially, provided $f$ is monotonically decreasing. In order to prove this, we restrict ourselves to $\mathbb{L}=\nabla\cdot(\mathbf{A}\nabla)$, where $\mathbf{A}\in\mathbb{R}^{n\times n}$ is a symmetric positive definite matrix, i.e., the coefficients of the elliptic differential operator $\mathbb{L}$ are supposed to be constant, and $r=0$.

Assume that $f^{\prime}\leq 0$, then from

[TABLE]

we obtain that $v=(u_{t})^{2}$ satisfies

[TABLE]

where $v_{0}=(\mathbb{L}u_{0}+f(u_{0}))^{2}$. Now if

[TABLE]

then the maximum principle shows that $\overline{v}\geq 0$ but more crucially, $v\leq\overline{v}$. For $\overline{v}$ we can use the representation formula (12), (13), (14), which yields

[TABLE]

and

[TABLE]

Now we will employ (52) in order to show that provided we have $f_{1},\;f_{2}\in W^{2,\infty}(I)$, then $\|\mathbb{T}(f_{1})-\mathbb{T}(f_{2})\|_{W^{1,\infty}(I)}\leq Ce^{-\lambda_{1}T}\|f_{1}-f_{2}\|_{W^{1,\infty}(I)}$ holds.

For $f_{1}$, $f_{2}\in X$, when taking the difference $\mathbb{T}(f_{1})-\mathbb{T}(f_{2})=f_{1}^{+}-f_{2}^{+}$, the terms containing $g$ explicitly and $r$ cancel, so that with $u^{(i)}(x,t)=u(x,t;f_{i})$ we have $f_{1}^{+}(g)-f_{2}^{+}(g)=D_{t}^{\alpha}(u^{(1)}-u^{(2)})=D_{t}^{\alpha}z$ where $z$ solves

[TABLE]

and we can estimate $\|(\mathbb{T}(f_{1})-\mathbb{T}(f_{2}))(g)\|_{Z}$ by estimating $D_{t}^{\alpha}z$.

In the parabolic case $\alpha=1$ this can be done by differentiating (53) with respect to time to obtain that $w=D_{t}z$ satisfies

[TABLE]

with $\tilde{y}=(f_{1}^{\prime}(Pu^{(1)})-f_{2}^{\prime}(Pu^{(2)}))P=({\textstyle\int_{0}^{1}f_{1}^{\prime\prime}(P(u^{(1)}+\theta z))\,d\theta}P\,z+(f_{1}-f_{2})^{\prime}(Pu^{(2)}))P$.

We assume that there exist potentials $\hat{q},\check{q}\in L^{\infty}(\Omega)$, $\hat{q},\check{q}\geq 0$ such that

[TABLE]

for some $\hat{c}_{1}$, $\check{c}_{1}$, such that

$4\hat{c}_{1}<\hat{\lambda}_{1}^{2}$, $4\check{c}_{1}<\check{\lambda}_{1}^{2}$,

where $\hat{\lambda}_{1}$ is the smallest eigenvalue of $\mathbb{L}_{\hat{q}}=\mathbb{L}+\check{q}$ (and analogously for $\check{q}$). We rewrite both initial boundary value problems in terms of $\mathbb{L}_{\hat{q}}$, and $\mathbb{L}_{\check{q}}$, respectively, and use corresponding eigensystems $\hat{\lambda}_{j}$, $\hat{\varphi}_{j}$, $\check{\lambda}_{j}$, $\check{\varphi}_{j}$. This yields

[TABLE]

hence abbreviating $df_{j}(s)=|(\tilde{y}(s)\,u^{(2)}_{t}(s),\hat{\varphi}_{j})|$

[TABLE]

where the first term on the right hand side can be estimated by

[TABLE]

This shows that we cannot take $r$ larger than one here, in order to avoid derivatives acting on the term $\tilde{y}(s)\,u^{(2)}_{t}(s)$ (note that the maximum principle argument for exponential decay of $u_{t}$ works for Lebesgue norms but not for higher order Sobolev norms). As a consequence, an embedding into $W^{1,\infty}(\Omega)$ will not be possible here and we have to work with a weaker norm on the differences $f_{1}-f_{2}$ which, however, should still be strong enough to capture the $f_{1}-f_{2}$ difference term appearing in $\tilde{y}$. We will use the norm $|||f_{1}-f_{2}|||:=\|(f_{1}-f_{2})(g)\|_{\dot{H}^{1}(\Omega)}$ and make the assumptions (36) and $u_{0}\in W^{1,\infty}(\Omega)$, which, together with the fact that by Theorem 3.1 $\|\nabla u^{(2)}\|_{L^{\infty}(0,T;L^{\infty}(\Omega))}\leq TC_{\dot{H}^{\sigma},W^{1,\infty}}^{\Omega}\rho$, implies the existence of constants $\bar{C}_{-1}(g)$, $\bar{C}_{0}(g)$, $\bar{C}_{1}(g)$ such that

[TABLE]

Altogether we have

[TABLE]

For $\eta(s)=\|w(s)\|_{\dot{H}^{r}(\Omega)}^{2}+\frac{4}{\hat{c}_{3}}\|\tilde{y}(s)\,u^{(2)}_{t}(s)\|_{\dot{H}^{r-1}(\Omega)}^{2}$ with $\hat{c}_{3}=\frac{4\hat{c}_{1}}{\hat{\lambda}_{1}}$, the above inequality becomes

[TABLE]

The convolution version of Gronwall’s inequality Lemma 2.5 yields

[TABLE]

hence

[TABLE]

Analogously to (55), we obtain for (53) the identity

[TABLE]

hence, using the replacements in (58)

[TABLE]

we obtain

[TABLE]

which for sufficiently large $\tilde{r}$ would even yield

[TABLE]

but actually we will only need $\tilde{r}=0$ below. This allows us to estimate the term $\tilde{y}u^{(2)}_{t}$ as follows

[TABLE]

Inserting this together with (52) into (58), setting $r=1$, and using our assumption (56), we obtain

[TABLE]

for C=C_{1}C_{2}+\Bigl{(}1+e^{-(\hat{\lambda}_{1}+\hat{c}_{3})T}\Bigr{)}\bar{C}_{0}(g)^{2} where C_{1}=4\Bigl{(}\|f_{1}^{\prime\prime}\|_{L^{\infty}(I)}\bar{C}_{-1}(g)+\bar{C}_{1}(g)\Bigr{)}^{2}, $C_{2}=(C_{\dot{H}^{\mu},L^{\infty}}^{\Omega})^{2}\|(\mathbb{L}u_{0}+f(u_{0}))^{2}\|_{{\dot{H}^{\mu}(\Omega)}}\tfrac{1}{\hat{\lambda}_{1}+\hat{c}_{3}}$.

Altogether we have proven:

Theorem 3.2**.**

*Let the assumptions of Theorem 3.1 and additionally $\alpha=1$, $r=0$, $f_{1},f_{2}\in W^{2,\infty}_{-}(I):=\{f\in W^{2,\infty}(I)\,:f^{\prime}\leq 0\mbox{ and }(\mathbb{L}u_{0}+f(u_{0}))^{2}\in{\dot{H}^{\mu}(\Omega)}\}$ for some $\mu\in(\frac{1}{2},2]$ hold.

Then there exists a constant $C>0$ depending only on $\|f_{1}^{\prime\prime}\|_{L^{\infty}(I)}$, $\|(\mathbb{L}u_{0}+f(u_{0}))^{2}\|_{{\dot{H}^{\mu}(\Omega)}}$ and the constants $\bar{C}_{-1}(g)$, $\bar{C}_{0}(g)$, $\bar{C}_{1}(g)$ in (56), such that*

[TABLE]

with $\hat{c}_{1}$ as in (54).

As a consequence of Theorem 3.2 we have the following conditional uniqueness result

Corollary 3.1**.**

*Let the assumptions of Theorem 3.2 be satisfied.

There is a $T^{\star}>0$ such that if $T>T^{\star}$, then there exists at most one fixed point of $\mathbb{T}$ in $W^{2,\infty}_{-}(I)$, thus a solution of the inverse problem (4), (2), (3) is unique in $W^{2,\infty}_{-}(I)$.*

3.4 Some thoughts and comments

It is instructive to compare the inverse problem for (1) with the initial-boundary conditions (2) studied in this paper with earlier work in the references [PilantRundell:1986, PilantRundell:1987, PilantRundell:1988] for similar equations and initial-boundary conditions. Both approaches used an iterative method obtained by projecting the equation onto a boundary where overposed data is available. In the current case this represents a spatial measurement of $u$ at a final time $t=T$, $u(x,T)$ for $x\in\Omega$, and in the former a time-measurement $u(x_{0},t)$ for $x_{0}\in\partial\Omega$. In each case the method relies on a mapping of $f$ to the differential operator projected onto the overposed boundary. In the current case this is $f\to u_{t}(x,T)$ and in the previous it is (in one space dimension) $f\to u_{xx}(x_{0},t)$. If one sets the problem in Hölder spaces $C^{k,\beta}$ then for $f\in C^{k_{0},\beta}$ the solution $u\in C^{k_{0}+2,\beta}(\Omega)\times C^{k_{0}+1,\beta/2}(0,T)$ Thus we have $u_{xx}(x_{0},t)\in C^{k_{0},\beta}$ and there is a relatively straightforward self mapping property. However, $u_{t}(x,T)\in C^{k_{0},\beta/2}$ causing difficulties due to the lack of regularity pickup as noted in the previous two subsections.

3.5 Data smoothing and noise propagation

Typically, we will be given noisy data $g^{\delta}$ and sometimes also some information on the statistics of the noise or, in the deterministic setting considered here, on the noise level with respect to the $L^{2}$ norm

[TABLE]

Application of the iterative scheme requires sufficiently smooth data, so instead of $g^{\delta}$ we employ a smoothed version $\tilde{g}^{\delta}$ satisfying the conditions on $g$ of Theorems 3.1, 3.2 as well as

[TABLE]

for some $m>0$, where $\tilde{\delta}\to 0$ as $\delta\to 0$, cf., e.g., [KaltenbacherRundell:2019a, Section 6.1].

In the situation of Theorem 3.1 with $m=\sigma$, we can make use of uniform boundedness by $\rho$ of the fixed point $f^{\tilde{g}^{\delta}}$ of the operator $\mathbb{T}_{\tilde{g}^{\delta}}$ corresponding to data $\tilde{g}^{\delta}$ as $\tilde{\delta\to 0}$, as well as weak* continuity of $\mathbb{T}$ to prove subsequential weak* convergence in $X$ of $f^{\tilde{g}^{\delta}}$ to $f_{act}$ as $\delta\to 0$, where $f^{\tilde{g}^{\delta}}$ is the fixed point corresponding to data $\tilde{g}^{\delta}$. Namely, uniform boundedness of $(\|f^{\tilde{g}^{\delta}}\|_{X})_{\delta\in(0,\bar{\delta}]}=(\|\mathbb{T}_{\tilde{g}^{\delta}}(f^{\tilde{g}^{\delta}})\|_{X})_{\delta\in(0,\bar{\delta}]}$ implies existence of weakly* in $W^{1,\infty}(I)$ and strongly in $L^{\infty}(I)$ convergent subsequences $\tilde{f}_{k}:=f^{\tilde{g}^{\delta}_{k}}$, $\mathbb{T}_{\tilde{g}^{\delta}_{k}}(f^{\tilde{g}^{\delta}_{k}})$ with limits $\bar{f}$ and $\hat{f}$, respectively. It remains to prove that $\bar{f}=\mathbb{T}_{g}(\bar{f})$, which means that $\bar{f}$ is a solution to the inverse problem with exact data $g$, provided it satisfies the range condition (31). To this end, we abbreviate $\tilde{f}_{k}:=f^{\tilde{g}^{\delta}_{k}}$ and $\tilde{g}_{k}:=\tilde{g}^{\delta}_{k}$ and consider the following decomposition

[TABLE]

where $y_{k}=D_{t}^{\alpha}u(\cdot,T;\tilde{f}_{k})-\mathbb{L}g-r(\cdot,T)$, $\tilde{y}_{k}=D_{t}^{\alpha}u(\cdot,T;\tilde{f}_{k})-\mathbb{L}\tilde{g}_{k}-r(\cdot,T)$,

$\mathbb{P}_{g}y\in\mbox{argmin}\{\|f^{\sharp}(g)-y\|_{Z}\,:\,f^{\sharp}\in X\mbox{ and }\|f^{\sharp}(u_{0})+r(\cdot,0)\|_{\dot{H}^{\sigma}(\Omega)}\leq\rho_{0}\}$, and we have used the fact that $\mathbb{T}_{\tilde{g}_{k}}(\tilde{f}_{k})=\mathbb{P}_{\tilde{g}_{k}}(\tilde{y}_{k})=\tilde{f}_{k}$. The first and last difference $\mathbb{T}_{g}(\bar{f})-\mathbb{T}_{g}(\tilde{f}_{k})$ and $\tilde{f}_{k}-\bar{f}$ tend to zero weakly* in $X$, due to the weak* continuity of $\mathbb{T}_{g}$. For the second and third differences, we abbreviate the projections defined by (32) by $f_{k}^{+}=\mathbb{P}_{g}(y_{k})$, $\tilde{f}_{k}^{+}=\mathbb{P}_{g}(\tilde{y}_{k})$, and $\tilde{\tilde{f}}_{k}^{+}=\mathbb{P}_{\tilde{g}_{k}}(\tilde{y}_{k})$, i.e., $f_{k}^{+}(g)=y_{k}$ and $\tilde{f}_{k}^{+}(g)=\tilde{y}_{k}=\tilde{\tilde{f}}_{k}^{+}(\tilde{g}_{k})$. Hence, the second difference in (62) satisfies

[TABLE]

and the third difference can be estimated by

[TABLE]

Altogether, the right hand side of (62) tends to zero as $k\to\infty$. A subsequence-subsequence argument therefore yields stability.

Proposition 3.1**.**

*Let $(g^{\delta})_{\delta\in(0,\bar{\delta}]}\subseteq L^{2}(\Omega)$ be a family of data satisfying (60) with corresponding smoothed versions $\tilde{g}^{\delta}$ satisfying (36) and (61) with $m=\sigma$ such that $\tilde{\delta}\to 0$ as $\delta\to 0$ and $\sigma>\frac{3}{2}$.

Then the sequence of fixed points $f^{\tilde{g}^{\delta}}$ of the perturbed operators $\mathbb{T}_{\tilde{g}^{\delta}}$ converges weakly* subsequentially in $X$ to the set of fixed points of the unperturbed operator $\mathbb{T}=\mathbb{T}_{g}$ as $\delta\to 0$.333i.e., every subsequence has an $X$ weakly* convergent subsequence and the limit of every $X$ weakly* convergent subsequence is a fixed point of $\mathbb{T}_{g}$. If the fixed point of $\mathbb{T}_{g}$ is unique, than the whole family $f^{\tilde{g}^{\delta}}$ converges weakly* in $X$ to this fixed point as $\delta\to 0$.*

Under the conditions of Theorem 3.2 with $T$ sufficiently large so that $Ce^{-\hat{\lambda}_{1}T/2}=:c<1$ and $m=3$, we can even provide a quantitative stability estimate as follows.

[TABLE]

hence

[TABLE]

4 Newton type methods

4.1 Operator formulation of nonlinearity identification in reaction-diffusion equation

Denote by $u_{act}(x,t)=u(x,t;f_{act})$ the exact solution of the inverse problem (4), (2), (3).

Again, we assume the range condition (31) to hold for $f=f_{act}$ and use the projection $P$ in the reformulation (35) of the direct problem. Moreover, as in the previous section, we use $Y=L^{2}(\Omega)$ as data space and a continuously embedded subspace $X$ of $W^{1,\infty}(I)$ as preimage space, e.g., just $X=W^{1,\infty}(I)$, or, in order to work with a Hilbert space setting, $X=H^{\sigma}(I)$ with $\sigma>\frac{3}{2}$) for $F:X\to Y$, $f\mapsto u(\cdot,T)$ where $u$ solves

[TABLE]

The operator $F$ is well-defined due to Theorem 2.3 (replacing $f$ there by $f\circ P$ here).

This allows to write the inverse problem as

[TABLE]

We denote the exact solution by $f_{act}$ and some initial guess by $f_{0}$.

The derivative of the forward operator is given by $F^{\prime}(f)h=v(\cdot,T)$ where $v$ solves

[TABLE]

Well-definedness of $v$ can be concluded from results in [SakamotoYamamoto:2011a], which has to be applied together with some fixed point argument, since the potential $f^{\prime}(Pu)$ depends on $x$ and $t$. Its adjoint (as appearing in the explicit version of Tikhonov regularized Newton or in source conditions) can be derived by means of the usual integration by parts procedure using the co-area formula. However, this leads to quite opaque formulas and hence we avoid it by stating Newton’s method and source conditions in a variational manner below.

Further regularity of $F$, such as Lipschitz continuity at $f_{act}$, can be verified by using the fact that $F^{\prime}(f_{act})-F^{\prime}(f)=w(\cdot,T)$ where $w$ solves

[TABLE]

where

[TABLE]

$u$, $v$ solve (63), (65), and $u_{act}$ solves (63), (65) with $f$ replaced by $f_{act}$, respectively, see Subsection 4.4 below.

4.2 Formulation of Newton type methods

For simplicity we again just use a deterministic noise bound

[TABLE]

while the noise itself might of course be random.

Consider two versions of regularized Newton iterations, cf., e.g., [Baku92, Hohage97, Jin:2010, KNS08, KaltenbacherPreviatti18, SchusterKaltenbacherHofmannKazimierski:2012, Werner15] and the references therein

Tikhonov regularized Newton:

[TABLE]

with a priori choice of the regularization parameter $\gamma_{k}=\gamma_{0}\kappa^{k}$ for some $\kappa\in(0,1)$.

Ivanov regularized Newton:

[TABLE]

with $\varrho$ chosen so that the exact solution is admissible, i.e., $\|f_{act}-f_{0}\|\leq\varrho$.

In order to guarantee existence of minimizers of (69) or (70) respectively, it is crucial to use a space $X$ that guarantees weak() compactness of sublevel sets as well as weak() lower semicontinuity of the cost functional (and constraint). Therefore a choice $X=H^{\sigma}(I)$ or $X=W^{1,\infty}(I)$ is admissible, while $X=C^{1}(I)$ would not satisfy these requirements.

Both versions of Newton’s method are stopped according to the discrepancy principle:

[TABLE]

for some $\tau>1$ (independent of $\delta$).

The advantage of (70) as compared to (69) is that it always guarantees a bounded sequence of iterates, no matter whether further conditions guaranteeing convergence are satisfied, see Section 4.3 below. Its drawback is that an upper bound of the distance between $f_{0}$ and $f_{act}$ should be known in order to have a reasonable choice of $\varrho$. If for all $\delta$ sufficiently small the sequence reaches the stopping index defined by the discrepancy principle, then weak compactness of $B_{\varrho}(f_{0})$ also guarantees weak subsequential convergence of $f_{k_{*}(\delta)}$ to a solution $f_{act}$, provided $F$ is weakly continuous.

4.3 Convergence of Newton type methods under a Lipschitz condition at $f_{act}$ and a variational source condition

Structural constraints on the forward operator, such as the tangential cone condition, are often not satisfied for parameter identification problems in PDEs unless complete observations of the state (i.e., in our case observations of $u$ on all of $\Omega\times(0,T)$) are available. So they are likely not satisfied for the inverse problem under consideration here.

In contrast to this, a Lipschitz condition on $F$

[TABLE]

can usually be verified in a much more standard way. Note that one of the two points where $F^{\prime}$ is evaluated is fixed to $f_{act}$, which enables to establish (72) without the need for imposing a stronger space $X$ (which would make the problem more ill-posed) but just assuming additional regularity on $f_{act}$, see Section 4.4.

A convergence proof of Newton type methods without structural condition (so with just the Lipschitz condition on $F^{\prime}$) requires $f_{act}$ to satisfy a source condition. The explicit version of the source condition would involve an adjoint, which we avoid by considering the variational version (cf., e.g., [HKPS07])

[TABLE]

for some $\beta<\frac{1}{L}$. Note that (73) implies that $f_{act}-f_{0}$ lies in the orthogonal complement of the kernel of $F^{\prime}(f_{act})$ so its verification is related to the uniqueness question for the linearized problem. We here restrict ourselves to the case of a Hilbert space $X$ with inner product $(\cdot,\cdot)_{X}$ for simplicity. The considerations below can be extended to a Banach space setting by using Bregman distances, (cf., e.g., [HKPS07]).

Convergence of (69) under conditions (72) and (73) so that the stopping index defined by the discrepancy principle (71) is reached after a finite number of steps, follows from existing results (e.g., [KNS08, Lemma 4.1], see also [Baku92] for the Hilbert space setting and [SchusterKaltenbacherHofmannKazimierski:2012, Werner15] in Banach spaces).

In case of (70), we can argue as follows, assuming that $\varrho$ has been chosen such that $\|f_{act}-f_{0}\|^{2}\leq\varrho\leq\|f_{act}-f_{0}\|^{2}+C\delta$ and omitting the norm subscripts to simplify the notation.

Minimality of $f_{k+1}$ and admissibility of $f_{act}$ implies

[TABLE]

where the left and right hand sides can further be estimated by means of the Lipschitz condition on $F^{\prime}$

[TABLE]

and

[TABLE]

so inserting $F_{p}(k)$ and $F_{pa}(k)$ into (74) we get

[TABLE]

Here the quadratic terms can be estimated by means of the variational source condition and the fact that $\|f_{\ell}-f_{0}\|^{2}\leq\varrho\leq\|f_{0}-f_{act}\|^{2}+C\delta$

[TABLE]

where we have again used the Lipschitz condition in the last estimate, hence

[TABLE]

Using this for $\ell=k$ and $\ell=k+1$ altogether yields

[TABLE]

i.e., an estimate of the form

[TABLE]

with $m=\frac{8\beta L}{1-4\beta\delta}$ smaller than one if $\beta$ is sufficiently small, so

[TABLE]

for

[TABLE]

which proves that $k_{*}(\delta)$ is finite provided $\tau$ has been chosen sufficiently large, $\tau>1+\tfrac{\bar{C}}{1-m}$.

Moreover, (71) together with (75) for $\ell=k_{*}(\delta)$ yields the convergence rate

[TABLE]

Proposition 4.1**.**

Let $Y$ be a Banach space, $X$ a Hilbert space, and $(g^{\delta})_{\delta\in(0,\bar{\delta}]}\subseteq Y$ a family of data satisfying (68) and let, for each $\delta\in(0,\bar{\delta}]$, the regularized approximation $f_{k_{*}(\delta)}$ be defined by an early stopped Newton type iteration according to either (69), (71) or (70), (71), with $\tau$ sufficiently large and for some differentiable forward operator $F:X\to Y$ satisfying (72), and some initial guess satisfying (73).

Then for each $\delta\in(0,\bar{\delta}]$, the stopping index is finite and the regularized approximations converge to the exact solution $f_{act}$ of $F(f)=g$ as appearing in (73) at the rate (76), i.e., $\|f_{k_{*}(\delta)}-f_{act}\|=O(\sqrt{\delta})$, $\|F(f_{k_{*}(\delta)})-g\|=O(\delta)$.

4.4 Back to the reaction-diffusion equation

We now verify the assumptions of Proposition 4.1 for the forward operator $F$ defined by (63).

4.4.1 Lipschitz continuity of $F^{\prime}$

Since the right hand side of (66), cf. (67), can be rewritten as

[TABLE]

we can remain with the assumption $X\subseteq W^{1,\infty}(I)$ on the preimage space, while assuming twice differentiability only on $f_{act}$. Assuming that $f^{\prime}(Pu_{act})P$ and $f^{\prime}(Pu)P$ are not too far away from an only space dependent function, so that there exists -q(x)\in\mathcal{D}(F):=\Bigl{\{}\tilde{q}\in L^{2}(\Omega)\,:\,\exists\bar{q}>0\,:\ \|\tilde{q}-\bar{q}\|_{L^{2}(\Omega)}<\mu\Bigr{\}}=\bigcup_{\bar{q}>0}\mathcal{B}^{L^{2}}_{\mu}(\bar{q}) (cf. [KaltenbacherRundell:2019b]) with

[TABLE]

for all $\vartheta\in[0,1]$ and some $c_{1}$ sufficiently small, we can rewrite (66), (65), and the equation for the difference $\hat{u}:=u_{act}-u$ as

[TABLE]

[TABLE]

and

[TABLE]

where these functions are equipped with homogeneous initial and impedance boundary conditions, where $\overline{f^{\prime}}=\int_{0}^{1}f_{act}^{\prime}(P(u_{act}+\vartheta(u-u_{act})))P\,d\vartheta$. Note that (78) holds for all $f\in B_{\varrho}(f_{0})$ provided

[TABLE]

and $\varrho$ is sufficiently small. Using eigenvalues and eigenfunctions of $-\mathbb{L}_{q}=-\mathbb{L}+q(x)$ we apply the representation formula (8). Together with (16) with $p=\infty$, $q=r=2$ and (41) for $\theta\in(0,2-1/\alpha)$, this yields, with $p^{*}=\frac{p}{p-1}$, $p^{**}=\frac{2p}{p-2}$,

[TABLE]

For $\|f_{act}^{\prime}(Pu_{act})P+q\|_{L^{2}(0,T;L^{p^{**}}(\Omega))}$ sufficiently small, this gives an estimate on the term $\|w\|_{C((0,T);L^{2}(\Omega))}$. Similarly, we get estimates on $\|\hat{u}\|_{L^{2}(0,T;L^{p^{*}}(\Omega))}$ $\|v\|_{L^{2}(0,T;L^{p^{*}}(\Omega))}$ appearing in the right hand side, using (16) with $p=r=2$, $q=1$ and $\int_{0}^{t}|\tau^{\alpha-1}E_{\alpha,\alpha}(-\lambda\tau^{\alpha})|\,dt\leq\frac{1}{\lambda}$,

[TABLE]

and likewise

[TABLE]

4.4.2 Variational source condition

Recall that $F^{\prime}(f_{act})h=v(\cdot,T)=:j$ where $v$ solves (65). Thus, using the operators $\mathbb{A}$, $\mathbb{B}$, $\mathbb{C}$, defined by $\mathbb{A}w=-\mathbb{L}w-f_{act}^{\prime}(g)w$ with with impedance boundary conditions and $g(x)=u_{act}(x,T)$, $\mathbb{B}h:=h(g)$, and $\mathbb{C}h=D_{t}^{\alpha}z(\cdot,T)$, where $z$ solves

[TABLE]

we can write

[TABLE]

Here we assume $q:=-f_{act}^{\prime}(g)$ to be nonnegative and contained in $L^{\infty}(\Omega)$ so that $\mathbb{A}$ is an isomorphism between $\dot{H}^{2}_{q}(\Omega)$ and $L^{2}(\Omega)$, where $\dot{H}^{2}_{q}(\Omega)$ is defined analogously to $\dot{H}^{s}(\Omega)$ cf. (6) with $\mathbb{L}_{q}:=\mathbb{L}-q$ in place of $\mathbb{L}$.

If $\mathbb{C}$ can be shown to be a contraction in the sense that

[TABLE]

then (73) holds with

[TABLE]

The spaces $V(\Omega)$ and $W(I)$ need to be appropriately chosen for this purpose. In particular we can choose the norms on $W(I)$ and $V(\Omega)$ weak in order to enable all these estimates to hold. (The weaker we choose $W(I)$, the stronger we have to pay for it by the smoothness condition of $\|f_{act}-f_{0}\|_{W(I)^{*}}$ being finite and small enough.)

Using a system of eigenvalues and eigenfunctions $\lambda_{j}^{q}$, $\varphi_{j}^{q}$ of $-\mathbb{L}+q$ we get

[TABLE]

hence

assuming $P(0)=0$, i.e., $g_{\min}\leq 0\leq g_{\max}$

[TABLE]

We therefore get, analogously to (42) and estimating

[TABLE]

that

[TABLE]

for $Q\leq\frac{2(2-\theta)}{2-2\theta+\sigma}$. Gronwall’s inequality (17) then yields

[TABLE]

where

[TABLE]

and

[TABLE]

Assuming

[TABLE]

for all $h\in X$ and $c_{2}$ or $c_{3}$ as well as $\|Pu_{act\,,t}\|_{L^{Q^{*}}(0,t;L^{2}(\Omega))}$ sufficiently small so that

[TABLE]

we obtain contractivity (80) with $V(\Omega)=\dot{H}^{\sigma}_{q}(\Omega)$.

From Proposition 4.1 we can conclude the following convergence result.

Corollary 4.1**.**

*Assume that $\Omega\subseteq\mathbb{R}^{1}$, $\alpha\in(\frac{4}{5},1]$, $f_{act}^{\prime}\leq 0$,

$g_{\min}\leq 0\leq g_{\max}$,

$\|f_{act}-f_{0}\|\leq\varrho$ sufficiently small and that (79), (81) holds with $c_{2}$ or $c_{3}$ and $\|Pu_{act\,,t}\|_{L^{Q^{*}}(0,t;L^{2}(\Omega))}$ sufficiently small so that (82) is satisfied. Let $F:X=\dot{H}^{\sigma}_{q}(\Omega)\to Y=L^{2}(\Omega)$ with $\sigma>\frac{3}{2}$ and $q$ as in (79) be defined by $F(f)=u(\cdot,T)$ where $u$ solves (63), and let $(g^{\delta})_{\delta\in(0,\bar{\delta}]}\subseteq Y$ be a family of data satisfying (68) and let, for each $\delta\in(0,\bar{\delta}]$, the regularized approximation $f_{k_{*}(\delta)}$ be defined by an early stopped Newton type iteration according to either (69), (71) or (70), (71), with $\tau$ sufficiently large.*

Then for each $\delta\in(0,\bar{\delta}]$, the stopping index is finite and the regularized approximations converge to the exact solution $f_{act}$ of $F(f)=g$ as appearing in (73) at the rate (76), i.e., $\|f_{k_{*}(\delta)}-f_{act}\|=O(\sqrt{\delta})$, $\|F(f_{k_{*}(\delta)})-g\|=O(\delta)$.

5 Reconstructions