Barbara Kaltenbacher

Department of Mathematics,
Alpen-Adria-Universität Klagenfurt

9020 Klagenfurt, Austria

William Rundell$^*$

Department of Mathematics,
Texas A&M University

College Station, Texas 77843, USA

(Communicated by Fioralba Cakoni)

This paper considers the inverse problem of recovering both the
unknown, spatially-dependent conductivity $a(x)$ and the nonlinear reaction
term $f(u)$ in a reaction-diffusion  equation
from overposed data.
These measurements can  consist of:
the value of two different solution measurements taken at a later time $T$;
time-trace profiles from two solutions;
or both final time and time-trace measurements from a single forwards
solve  data run.
We prove both uniqueness results and the convergence of iteration
schemes designed to recover these coefficients.
The last section of the paper shows numerical reconstructions based on these

§ INTRODUCTION

In their most basic form reaction-diffusion equations are nonlinear
parabolic equation of the type
$ u_t(x,t) + \mathbb{L}u = f(u)$ where $\,\mathbb{L}$ is an elliptic operator
defined on a domain $\Omega\subset \mathbb{R}^d$, $d\geq 1$
and $f(u)$ is a function only of the state or dependent variable $u$
which is defined on $\Omega\times [0,t]$ for some $T>0$.

Such equations arise in a wide variety of applications
and two examples of these models are in ecology where $u$
represents the population
density at a fixed point $x$ and time $t$ and $f(u)$ is frequently taken
to be quadratic in $u$ as in the Fisher model;
or in chemical reactions where $f$ is cubic in
the case of the combustion theory of Zeldovich and Frank-Kamenetskii,
[12, 23, 19].
In these models $\mathbb{L}$ is of second order but higher order operators
also feature in significant applications, for example in the Cahn-Hilliard
equation [5, 8].

The physical interpretation of these equations is the time rate
of change of  $u$ is a sum of two components:
the diffusion term given by $\mathbb{L}u$ and the driving term $f(u)$.
Often these are in competition, the former tending to disperse the value
of $u$ and the latter, if $f>0$, to increase it.
If $f$ is negative or positive but sublinear, then on a bounded domain
the solution will usually decay exponentially (but depending on the boundary
conditions/values and the eigenvalues of the combined operator $\mathbb{L}$
being positive).
At the other extreme if the growth rate of $f$ is sufficiently large
then global existence in time of the solution can fail.
This is a complex but well-documented situation, see for example,
[1, 20] and the references within.
More general references for semilinear parabolic equations
are [9, 10].

In most applications one assumes these terms in the equation are known,
that is, the coefficients in the elliptic operator and the exact form
of $f(u)$.  Indeed, frequently, $f(u)$ is taken to be a polynomial
or possibly a rational function of low degree so that only a few parameters
have to be determined by a least squares fit to the data.
Even if the diffusion constant has also to be determined the method remains
fundamentally the same.
In this paper we make no such assumptions: both our diffusion coefficient and
the nonlinear term can be arbitrary functions within a given smoothness class.
We shall take $\mathbb{L}u = -\nabla\cdot\bigl(a(x)\nabla u\bigr)$ where
the diffusion coefficient $a(x)$ is unknown.

With this background setting our basic model equation is thus

\begin{equation}\label{eqn:basic_pde_parabolic}
u_t(x,t) -\nabla\cdot(a(x)\nabla u(x,t)) = f(u) + r(x,t,u)
\end{equation}

where $r(x,t,u)$ is a known forcing function.
Boundary conditions for (<ref>)
will be of the impedance form

\begin{equation}\label{eqn:bdry_cond}
a\frac{\partial u}{\partial\vec{n}} + \gamma u = {b},\qquad
x\in\partial\Omega,\quad t>0,
\end{equation}

with possibly space dependent but known coefficient $\gamma$ (where $\gamma=\infty$ indicates the case of Dirichlet boundary conditions) as well as given space and time dependent right hand side $b$.
We impose the initial condition

\begin{equation}\label{eqn:init_cond}
u(x,0) = u_0(x),\qquad x\in\Omega.
\end{equation}

Now of course the recovery of the coefficient $a$ and the term $f(u)$ requires
over-posed data and we shall consider several different forms and combinations.

The two basic forms are final time and boundary measurements of $u$.
The former consists of the value $g(x) = u(x,T)$ at some later time $T$
and corresponds to such conditions as a thermal or concentration  map of $u$
within the region $\Omega$ or census data taken at a fixed time.
The latter consists of a time trace $h(t) = u(\tilde x,t)$ for some fixed
point $\tilde x$ on $\partial\Omega$.

Combinations include measurements of one of these for two different
sets of initial / boundary conditions corresponding to solutions $u(x,t$ and
$v(x,t)$ or a measurement of both of these for a single data run.
As we will see, some of these combinations are superior to others and the
difference can be substantial.

The main technique we will use follows that of previous work by the
authors and others:  we project the differential equation onto the
boundaries containing the overposed data and obtain nonlinear equations
for our unknowns $a(x)$ and $f(u)$ which we then solve by iterative methods.
This idea for recovering the nonlinearity $f(u)$ in equations such as
was first used in [26, 25]
under time-trace data.
References for recovering coefficients in parabolic equations in general,
but specifically including final-time data,
are [13, 14].
A recent such problem for the subdiffusion case based on fractional derivatives
is found in [31].

More recently, the authors have considered several inverse problems
based on the above equation.
In [15] the equation taken was
$u_t - \triangle u = q(x)g(u)$ where $g(u)$ was assumed known and $q(x)$
which determines the unknown spatial strength of $g$, had to be recovered
from final time data.
In [16] the equation taken was
$u_t - \triangle u = f(u)$ and recovery of the unknown $f(u)$ was again from
final time data.
The case of $\mathbb{L}u = -\nabla\cdot\bigl(a(x)\nabla u\bigr) + q(x)u$
where both  $a(x)$ and $q(x)$ are unknown and have to be determined
from a pair of runs each providing final-time data was considered in
In this work the equation taken was linear, so that $f(u)$ did not appear.

Our main tool is based on the technique of projecting the differential
operator onto that segment of the boundary where overposed data lies
and thereby obtain an operator equation suitable for fixed point iteration
and the attendant methods for showing convergence of such schemes.
Other possibilities exist and have been widely used such as
Newton-type schemes.
For the problems at hand these techniques are certainly viable candidates
but require, often considerably, more complex analysis,
in particular, convergence conditions that are likely
not satisfied by highly nonlinear problems.
In a departure from most of the recent work in this area we will
use a Schauder space setting for the parabolic equation.
This offers several advantages including the fact that regularity estimates
for parabolic equations set in Schauder space are independent of spatial

In each of the above papers the diffusion process was extended to include
the subdiffusion case where the time-derivative was of fractional type.
The very rough picture that evolves from this modification is no change
to the uniqueness questions of recovery of terms but possibly
substantial changes in the degree of ill-conditioning, especially when
either only very short or when very long time values are present.
Certainly, the analysis is more complex when fractional derivatives are
In the present work we will consider only the classical parabolic situation
as we don't see sufficient novelty of outcome in the central ideas
to merit the analytical complexities that a comprehensive inclusion would
entail although we will sometimes comment on particular instances.

We make a simple, but essential, overview comment on inverse problems
seeking to recover a term that involves the dependent variable $u$.
One can only recover such a term over a range of values that are contained
in the data measurement range and, further, during the evolution process
the range of the solution $u$ cannot have taken values outside of this range.
This makes such problems quite distinct from determining coefficients
where the required range is known from the geometry of the domain $\Omega$.
From a practical perspective this means that the experimental set up
must take this point into consideration and an arbitrary mix of initial
and boundary conditions and values is unlikely to satisfy this range

There is some way around this issue by representing the unknown function,
here $f(u)$, in a chosen basis set.
If $f$ is known to have a simple form then a few well-chosen elements
may suffice, but otherwise the required extrapolation will lead to
a severely ill-posed problem.
We are not considering this case as we seek a methodology that includes
quite complex situations.
In addition, the iteration schemes we propose are best formulated
in a pointwise form and this would be precluded without a range condition.

In the current work we seek to determine a pair of unknown coefficients.
One of these, $a(x)$ is only spatially dependent but strongly influences the
resulting solution profile.
The reaction term $f(u)$ also might have a strong influence on the
solution profile but
it can also be such that it is small in magnitude and its influence on the
solution is dominated by the diffusion process governed by $a(x)$.
This case is actually more difficult than when the two terms have roughly
equal effect on the solution.
In either situation the range condition must still be in place.
In section sect:recons we will see a situation where violation of
this can occur in a rather subtle way leading to incorrect reconstructions.

The plan of the paper is to introduce the iterative schemes that will
be used to recover both $a$ and $f$ in section sec:conv
and present conditions under which we are able to show a convergence analysis.
This section is broken down into three cases depending on the data type being
measured: two runs with different initial/boundary data and each
measuring a final time $g(x) = u(x,T)$ for some fixed $t=T$;
two runs each measuring a time trace of the solution $u$ at a fixed
point $\tilde x\in\partial\Omega$; and a single run with measurements of
both a time trace and a final time.
We will also provide a convergence analysis of each of these schemes.
The following section sect:recons  will show actual reconstructions
using these iterative schemes and demonstrate both their capabilities and their

Before proceeding to the above agenda we should make some comments
on physical scaling as this will be very relevant to our reconstructions in
section sect:recons.
${\mathbb L}$ has been written with a single coefficient $a(x)$ but
more general possibilities would include:
$\;{\mathbb L}u := -\nabla\cdot(a\nabla u) + d(x)\nabla_i u + q(x) u$.
In this setting $d$ would correspond to a drift term and $q$ a potential
or as a linear forcing term with a space-dependent modifier of its
Equation (<ref>) then appears with several
terms all of which have specific physical attributes and in consequence
scale factors depending on the situation.
By this we mean that depending on the context the range of value
taken by the various terms, including our assumed unknowns $a(x)$ and $f(u)$
may have enormous differences.
Examples for the diffusion coefficient range from around
$10^{-6}$ to $10^{-5}$ metres$^2$/sec for molecules in a gas to about
$10^{-10}$ metres$^2$/sec for molecules in a liquid such as water,
and in this case the effect of drift represented by $d(x)$
can overshadow the diffusion effect
(of course the domain size of $\Omega$ also plays a role).
Because of this we must be careful in interpreting the results of
numerical experiments which have taken generic order unity values
as the default.

§ CONVERGENCE ANALYSIS

We consider the inverse problem of recovering the spatially varying diffusion $a(x)$ and the nonlinearity $f(u)$ in the context of three different observation settings:

* observation of two states (driven by different initial and boundary conditions as well as forcing terms) at the same final time $T$

* observation of one or two states at final time and at a (boundary or interior) point $\tilde{x}$ over time

* observation of two different states at $\tilde{x}$ over time

We discuss each of these cases in a separate subsection below.

In the first two cases we can follow the approach of projecting the PDE onto the
observation manifold; in the third case this is only partly possible, since
time trace measurements do not map directly into space dependence (of $a$).
While the methods themselves clearly extend to higher space dimensions in the
first two cases (in the third case they cannot, by an elementary dimension count)
we provide an analysis in one space dimension $\Omega=(0,L)$ only.
Possibilities and limitations with respect to to higher space dimensions will
be indicated in some remarks.

Since well-definedness of the methods to be discussed relies on existence, uniqueness and regularity of solutions to the semilinear parabolic initial boundary value problem (<ref>), (<ref>), (<ref>), we will provide a statement from the literature [9] on well-posedness of this forward problem.
For our purposes it suffices to consider the case of homogeneous boundary conditions ($b=0$ in (<ref>)). As a matter of fact, it would be enough to look at the spatially 1-d case for which we do the analysis; however this would not change so much in the assumptions, since the Schauder space results in [9] are independent of the space dimensions.
In case $\gamma=\infty$ (Dirichlet boundary value problem) we can directly make use of <cit.> for existence and of <cit.> for uniqueness.
In case of Neumann or impedance conditions $\gamma<\infty$, these can be extended by means of the results from <cit.>, noting that these corresond to the second boundary value problem in [9] with $\beta=\frac{\gamma}{a}$, $g=\frac{b}{a}$.

Let $\Omega$ be a $C^{2+s}$ domain for some $s>0$, let $f$ and $r$ be Hölder continuous on bounded sets and $f$ satisfy the growth conditions
$uf(u)\leq A_1 u^2+A_2$
$|f(u)|\leq A(u)$ for some positive monitonically increasing function $A$
and let $a$, $\nabla a$ be Hölder continuous (note that the differential operator in [9] is in nondivergence form) with $a$ being positive and bounded away from zero on $\overline{\Omega}$.
Moreover, let $u_0\in C^{2+s}(\overline{\Omega})$ satisfy the compatibility conditions
$u_0=0$ and $\nabla\cdot(a(x)\nabla u_0(x)) = f(u_0(x)) + r(x,0,u_0)$ for $x\in\partial\Omega$ if $\gamma=\infty$,
or vanish in a neighborhood of $\partial \Omega$ in case $\gamma<\infty$.

Then there exists a classical solution $u$ with $u_t,\,u_{x_i},\,u_{x_ix_j}$ Hölder continuous of (<ref>), (<ref>), (<ref>) with $b=0$.

If additionally $f$ is Lipschitz continuous, then this solution is unique.

Note that the results taken from [9] deal with solutions on bounded time domains, which allows to impose less restrictive growth conditions than those needed on all of $[0,\infty)$ as a time interval.
Since we will consider functions $f$ that are nonconstant on a bounded domain only, Lipschitz continuity of $f$ is sufficient for the smoothness and growth conditions in Theorem th:wellposedforward.

§.§ Final time observations of two states

Identify $a(x)$, $f(u)$ in

\begin{eqnarray}
&&u_t-(a u_x)_x=f(u)+r^u \quad t\in(0,T)\,, \quad u(0)=u_0\label{eqn:u-f}\\
&&v_t-(a v_x)_x=f(v)+r^v \quad t\in(0,T)\,,  \quad v(0)=v_0\label{eqn:v-f}
\end{eqnarray}

with homogeneous impedance boundary conditions

\begin{equation}\label{eqn:bc}
a\partial_{\vec{n}} u+\gamma^u u =0\,, \quad a\partial_{\vec{n}} v+\gamma^v v =0 \quad\mbox{ on }\partial\Omega
\end{equation}

from final time data

\begin{equation}\label{eqn:finaltimedata}
g^u(x)=u(x,T)\,, \quad g^v=v(x,T)\,, \quad x\in\Omega,
\end{equation}

where $\Omega=(0,L)$.

In order for this data to contain enough information to enable identification
of both $a$ and $f$, it is essential that the two solutions $u$ and $v$ are
sufficiently different, which is mainly achieved by the choice of the two
different driving terms $r^u$ and $r^v$, while the boundary conditions and
later on also the initial data will just be homogeneous.

Projecting the PDEs (<ref>), (<ref>) onto the measurement
manifold $\Omega\times\{T\}$ and denoting by $u(x,t;a,f)$, $v(x,t;a,f)$ their
solutions with coefficients $a$ and $f$, we obtain equivalence of the inverse
problem (<ref>)–(<ref>) to a fixed point equation
with the fixed point operator
$\mathbb{T}$ defined by $(a^+,f^+)=\mathbb{T}(a,f)$ with

\begin{equation}\label{eqn:fp_fiti-fiti}
\begin{aligned}
&(a^+ g^u_x)_x(x)+f^+(g^u(x))=D_tu(x,T;a,f)-r^u(x,T) \quad x\in\Omega\\
&(a^+ g^v_x)_x(x)+f^+(g^v(x))=D_tv(x,T;a,f)-r^v(x,T) \quad x\in\Omega\,.
\end{aligned}
\end{equation}

Thus, the iterates are implicitly defined by

\[
\mathbb{M}(a^+,f^+)=\mathbb{F}(a,f) +\mathrm{b}\,,
\]

with the linear operator $\mathbb{M}$, the nonlinear operator $\mathbb{F}$, and the inhomogeneity $\mathrm{b}$ defined by

\begin{equation}\label{eqn:MFb}
\begin{aligned}
\mathbb{M}(a^+\!,f^+)=
\left(\begin{array}{c}\!(a^+ g^u_x)_x+f^+(g^u)\\(a^+ g^v_x)_x+f^+(g^v)\end{array}\right)\,, \quad
\\
\mathbb{F}(a,f)=
\left(\begin{array}{c}u_t(\cdot,T;a,f)\\ v_t(\cdot,T;a,f)\end{array}\right)\,, \quad
\mathrm{b}=
-\left(\begin{array}{c}\!r^u(\cdot,T)\!\\ \!r^v(\cdot,T)\end{array}\right)\,.
\end{aligned}\end{equation}

The question of self mapping and contractivity of $\mathbb{T}$ in a neighbourhood of an exact solution $(a_{ex},f_{ex})$ leads us to consider the differences $\hat{a}^+=a^+-a_{ex}$, $\hat{f}^+=f^+-f_{ex}$, $\hat{a}=a-a_{ex}$, $\hat{u}=u(\cdot,\cdot;a,f)-u_{ex}$, $\hat{v}=v(\cdot,\cdot;a,f)-v_{ex}$,  along with the identity

\[
\mathbb{M}(\hat{a}^+,\hat{f}^+)=\mathbb{F}(a,f)-\mathbb{F}(a_{ex},f_{ex})
\]

(note that the inhomogeneities cancel out)
so that

\begin{eqnarray*}
&&(\hat{a}^+ g^u_x)_x(x)+\hat{f}^+(g^u(x))=D_t\hat{u}(x,T) \quad x\in\Omega \\
&&(\hat{a}^+ g^v_x)_x(x)+\hat{f}^+(g^v(x))=D_t\hat{v}(x,T) \quad x\in\Omega \,.
\end{eqnarray*}

As in [17], the convergence estimates consist of two steps:

(a) prove bounded invertibility of $\mathbb{M}$, i.e., estimate appropriate norms of $a^+$, $f^+$

\begin{eqnarray}
&&(\hat{a}^+ g^u_x)_x(x)+\hat{f}^+(g^u(x))=\check{u}(x,T) \quad x\in\Omega\label{eqn:diffu}\\
&&(\hat{a}^+ g^v_x)_x(x)+\hat{f}^+(g^v(x))=\check{v}(x,T) \quad x\in\Omega\label{eqn:diffv}
\end{eqnarray}

by appropriate norms of $\check{u}$, $\check{v}$.

(b) prove smallness of $\mathbb{F}(a,f)-\mathbb{F}(a_{ex},f_{ex})$, i.e., estimate $\check{u}$, $\check{v}$ in terms of the above chosen norms of $\hat{a},\hat{f}$, using the fact that $\check{u}$, $\check{v}$ satisfy certain

Step (a), bounding $\mathbb{M}^{-1}$:
Unlike that of [17],
the elimination strategies are limited.
In this reference  $\hat{q}^+$ was eliminated in order to estimate
$\hat{a}^+$ and then $\hat{a}^+$ was eliminated for estimating $\hat{q}^+$.
In the current situation we cannot eliminate $\hat{f}^+$ any more;
but we can in fact eliminate $\hat{a}^+$ and express $\hat{f}^+$ in terms of $\check{u}$ and $\check{v}$.
This step is considerably more complicated than the $\hat{a}^+$
elimination in [17]
and is carried out by the following steps.

We eliminate $\hat{a}^+$ by integrating (<ref>), (<ref>) with respect to $x$ and multiplying with $g^v_x(x)$ and $-g^u_x(x)$, respectively,

\[
\begin{aligned}
&g^v_x(x) \int_x^L \hat{f}^+(g^u(\xi))\, d\xi - g^u_x(x) \int_x^L \hat{f}^+(g^v(\xi))\, d\xi \\
&= g^v_x(x) \int_x^L \check{u}(\xi,T)\, d\xi - g^u_x(x) \int_x^L \check{v}(\xi,T)\, d\xi
\end{aligned}
\]

Here we have avoided boundary terms involving $\hat{a}$ by assuming that either $a(L)$ is known (so that $\hat{a}(L)=0$) or homogeneous Neumann boundary conditions are imposed on the right hand boundary (so that $g^u_x(L)=g^v_x(L)=0$).
In order to arrive at a first kind Volterra integral equation for $\hat{f}^+$, we assume that $g^u$ is strictly monotone (w.l.o.g. increasing) with

\begin{equation}\label{eqn:mu}
g^u_x(x)\geq\mu>0  \quad x\in\Omega
\end{equation}

and divide by $g^u_x(x)$ to get

\[
\begin{aligned}
&\tfrac{g^v_x}{g^u_x}(x) \int_x^L \hat{f}^+(g^u(\xi))\, d\xi - \int_x^L \hat{f}^+(g^v(\xi))\, d\xi
= \tfrac{g^v_x}{g^u_x}(x) \int_x^L \check{u}(\xi,T)\, d\xi - \int_x^L \check{v}(\xi,T)\, d\xi
\end{aligned}
\]

and after differentiation

\begin{equation}\label{eqn:fhatplusVolterra}
\begin{aligned}
& \hat{f}^+(g^v(x))
-\tfrac{g^v_x}{g^u_x}(x) \hat{f}^+(g^u(x))
+(\tfrac{g^v_x}{g^u_x})_x(x)\int_x^L \hat{f}^+(g^u(\xi))\, d\xi
\\
&= \check{v}(x,T)
- \tfrac{g^v_x}{g^u_x}(x) \check{u}(x,T)
+(\tfrac{g^v_x}{g^u_x})_x(x)\int_x^L \check{u}(\xi,T)\, d\xi \,.
\end{aligned}
\end{equation}

Further differentiations of (<ref>) yield the identities

\begin{equation}\label{eqn:fhatplusprime}
\begin{aligned}
& (\fhatplusp (g^v(x)) - \fhatplusp (g^u(x))) g^v_x(x)
-2(\tfrac{g^v_x}{g^u_x})_x(x) \hat{f}^+(g^u(x))
+(\tfrac{g^v_x}{g^u_x})_{xx}(x)\int_x^L \hat{f}^+(g^u(\xi))\, d\xi
\\
&= \check{v}_x(x,T)
- \tfrac{g^v_x}{g^u_x}(x) \check{u}_x(x,T)
- 2(\tfrac{g^v_x}{g^u_x})_x(x) \check{u}(x,T)
+(\tfrac{g^v_x}{g^u_x})_{xx}(x)\int_x^L \check{u}(\xi,T)\, d\xi \,.
\end{aligned}
\end{equation}

\begin{equation}\label{eqn:fhatplusprimeprime}
\begin{aligned}
(&\fhatpluspp(g^v(x))g^v_x(x) - \fhatpluspp(g^u(x))g^u_x(x)) g^v_x(x)\\
&\quad=-(\fhatplusp (g^v(x)) - \fhatplusp (g^u(x))) g^v_{xx}(x)
+2(\tfrac{g^v_x}{g^u_x})_x(x)g^u_x(x)\, \fhatplusp (g^u(x))\\
&\qquad+3(\tfrac{g^v_x}{g^u_x})_{xx}(x) \hat{f}^+(g^u(x))
-(\tfrac{g^v_x}{g^u_x})_{xxx}(x)\int_x^L \hat{f}^+(g^u(\xi))\, d\xi
\\
&\qquad+ \check{v}_{xx}(x,T)
- \tfrac{g^v_x}{g^u_x}(x) \check{u}_{xx}(x,T)
- 3(\tfrac{g^v_x}{g^u_x})_x(x) \check{u}_x(x,T)\\
&\qquad- 3(\tfrac{g^v_x}{g^u_x})_{xx}(x) \check{u}(x,T)
+(\tfrac{g^v_x}{g^u_x})_{xxx}(x)\int_x^L \check{u}(\xi,T)\, d\xi \\
&\quad =:\ \Phi(\fhatplusp ,{\hat{f}^+})(x)+ b(x)\,.
\end{aligned}
\end{equation}

To analyse this integral equation we assume that the slopes of $g^u$ and $g^v$ are sufficiently different in the sense that

\begin{equation}\label{eqn:kappa}
\left|\frac{g^v_x(x)}{g^u_x(x)}\right|\leq \kappa<1 \quad x\in\Omega
\end{equation}

and at the same time still ensuring that

\begin{equation}\label{eqn:delta}
g^v_x(x)\geq\delta  \quad x\in\Omega.
\end{equation}

Moreover, since the strategy is to control the full $C^2$ norm by bounding $\fhatpluspp$, we assume $f(g^u(0))$, $f(g^v(0))$, $f'(g^u(0))$, $f'(g^v(0))$ to be known so that

\begin{equation}\label{eqn:fhatplus0}
{\hat{f}^+}(g^u(0))={\hat{f}^+}(g^v(0))=\fhatplusp (g^u(0))=\fhatplusp (g^v(0))=0.
\end{equation}

As was found in
[16], [18],
we have to impose range conditions

\begin{equation}\label{eqn:rangecondition_finaltime}
J=g^u(\Omega)\supseteq g^v(\Omega)\,, \quad
J\supseteq u_{ex}(\Omega\times(0,T))\,, \quad J\supseteq v_{ex}(\Omega\times(0,T))\,.
\end{equation}

Of course the roles of $g^u$ and $g^v$ could be reversed here. Note however, that the assumption of $g^u$ being the data with the larger range conforms with the condition (<ref>) of $g^u$ being the function with the steeper slope.

Let $\beta\in[0,1]$, $g^u,g^v\in C^4(\Omega)$, $\check{u}(T),\check{v}(T)\in C^{2,\beta}(\Omega)$, and assume that (<ref>), (<ref>), (<ref>), (<ref>) hold.
Then (<ref>) is uniquely solvable in
$C^{2,\beta}(J)$ and

\[
\|\fhatpluspp\|_{C^{0,\beta}(J)}
\leq C(g^u,g^v) \Bigl(\|\check{u}(T)\|_{C^{2,\beta}(\Omega)}+\|\check{v}(T)\|_{C^{2,\beta}(\Omega)}\Bigr)
\]

and therefore, by (<ref>),

\begin{equation}\label{eqn:estf}
\begin{aligned}
\|\hat{f}^+\|_{C^{2,\beta}(J)}\leq \tilde{C} \|\fhatpluspp\|_{C^{0,\beta}(J)}
&\leq \bar{C}(g^u,g^v) \Bigl(\|\check{u}(T)\|_{C^{2,\beta}(\Omega)}+\|\check{v}(T)\|_{C^{2,\beta}(\Omega)}\Bigr)\,,
\end{aligned}\end{equation}

where $\bar{C}(g^u,g^v)$ only depends on smoothness of $g^u,g^v$, as well as on $\kappa, \delta, \mu$.

By (<ref>), (<ref>) we get, after division by $-g^u_x(x) g^v_x(x)$,

\[
\hat{f}^{+\,\prime\prime}(g^u(x))-\hat{f}^{+\,\prime\prime}(g^v(x))\tfrac{g^v_x}{g^u_x}(x)
= -\frac{1}{g^v_x(x)g^u_x(x)} \Bigl(\Phi(\hat{f}^{+\,\prime},{\hat{f}^+})(x)+b(x)\Bigr)\,.
\]

Hence, taking the supremum over $x\in\Omega$ and using (<ref>)
we obtain, using (<ref>), (<ref>), (<ref>),

\[
\|\hat{f}^{+\,\prime\prime}\|_{C(J)}\leq \kappa \|\hat{f}^{+\,\prime\prime}\|_{C(J)} + \frac{1}{\mu\delta}
\|\Phi(\hat{f}^{+\,\prime},{\hat{f}^+})+b\|_{C(\Omega)}
\]

that is,

\[
\|\hat{f}^{+\,\prime\prime}\|_{C(J)}\leq \frac{1}{\mu\delta(1-\kappa)}
\|\Phi(\hat{f}^{+\,\prime},{\hat{f}^+})\|_{C(\Omega)}\,.
\]

By (<ref>) we can write

\[
\begin{aligned}
\hat{f}^{+\,\prime}(g^u(x))&=\int_0^x \hat{f}^{+\,\prime\prime}(g^u(\xi))g^u_x(\xi)\, d\xi\,, \\
%=\int_0^x \hat{f}^{+\,\prime}(g^u(\xi))g^u_x(\xi)\, d\xi
&=\int_0^x \int_0^\xi \hat{f}^{+\,\prime\prime}(g^u(\sigma))g^u_x(\sigma)\, d\sigma g^u_x(\xi)\, d\xi \,,\\
\int_x^L {\hat{f}^+}(g^u(\xi)) \, d\xi
&= \int_x^L \int_0^\xi \int_0^\sigma \hat{f}^{+\,\prime\prime}(g^u(\tau))g^u_x(\tau)\, d\tau g^u_x(\sigma)\, d\sigma \, d\xi\,,
\end{aligned}
\]

\[
\Phi(\hat{f}^{+\,\prime},{\hat{f}^+}) = K \hat{f}^{+\,\prime\prime}
\]

with a compact operator $K:C(J)\to C(\Omega)$.
On the other hand, with the isomorphism  $B:C(J)\to C(\Omega)$ defined by

\[
(Bj)(x):= (j(g^v(x))g^v_x(x)-j(g^u(x))g^u_x(x))g^v_x(x)
= -\Bigl(j(g^u(x))-j(g^v(x))\tfrac{g^v_x(x)}{g^u_x(x)}\Bigl)g^u_x(x) g^v_x(x)
\]

(and noting that $\,\|B^{-1}\|_{C(\Omega)\to C(J)}\leq \frac{1}{\mu\delta(1-\kappa)}$, see above),
we can write the problem of recovering $\hat{f}^{+\,\prime\prime}$ as a second kind Fredholm equation

\[
j - B^{-1}K j = B^{-1}b
\]

for $j=\hat{f}^{+\,\prime\prime}$ and apply the Fredholm alternative in $C(J)$. To this end, we have to prove that the kernel of $I-B^{-1}K$ is trivial. For this purpose we have to be able to conclude $j\equiv0$ from
$J\equiv0$ where

\[
\begin{aligned}
&J(x)=: j\bigl((g^v(x))g^v_x(x) -(g^u(x))g^u_x(x)\bigr) g^v_x(x)\\
&\quad+\int_0^x\! j(g^v(\xi))g^v_x(\xi)\, d\xi -
g^v_{xx}(x)\int_0^x\! j(g^u(\xi))g^u_x(\xi)\, d\xi
\\
&\quad-2(\tfrac{g^v_x}{g^u_x})_x(x) g^u_x(x) \int_0^x\!j(g^u(\xi))g^u_x(\xi)\, d\xi -3(\tfrac{g^v_x}{g^u_x})_{xx}(x)
\int_0^x \!\int_0^\xi\!j(g^u(\sigma))g^u_x(\sigma)\, d\sigma g^u_x(\xi)\, d\xi
\\
\int_x^L\!\int_0^\xi\!\int_0^\sigma\!j(g^u(\tau))g^u_x(\tau)\, d\tau g^u_x(\sigma)\, d\sigma \, d\xi.
\end{aligned}
\]

We then have

\[
\begin{aligned}
J(x) =& j\bigl(g^v(x))g^v_x(x) - (g^u(x))g^u_x(x)\bigr) g^v_x(x)
+\int_{g^v(0)}^{g^v(x)}\!j(z)\, dz\\
&\quad- (g^v_{xx}(x)+2(\tfrac{g^v_x}{g^u_x})_x(x))
\int_{g^u(0)}^{g^u(x)}\!j(z)\, dz
\int_{g^u(0)}^{g^u(x)}\!\int_{g^u(0)}^{z} j(y)\, dy \, dz
\\
\int_x^L
\int_{g^u(0)}^{g^u(\xi)}\!\int_{g^u(0)}^{z} j(y)\, dy \, dz \, d\xi \,.
\end{aligned}
\]

After division by $-(g^u_x\cdot g^v_x)(x)$ and the change of variables $\zeta=g^u(x)$, this is equivalent to

\[
\begin{aligned}
0=& j(\zeta)-\mathrm{a}(\zeta)\,j(\mathrm{b}(\zeta))
+\int_{\underline{g}}^\zeta \mathrm{k} (\zeta,z) j(z)\, dz\\
&+\mathrm{c}(\zeta)\int_{g^u(0)}^\zeta \int_{g^u(0)}^{z} j(y)\, dy \, dz
+\mathrm{d}(\zeta)\int_\zeta^{g^u(L)}\!\frac{1}{g^u_x((g^u)^{-1}(r)}\int_{g^u(0)}^\zeta \int_{g^u(0)}^{z} j(y)\, dy \, dz\, dr
\end{aligned}
\]

\[
\begin{aligned}
&\mathrm{a}(g^u(x))=(\tfrac{g^v_x}{g^u_x})(x), \qquad\mathrm{b}(g^u(x))=g^v(x), \\
\frac{g^v_{xx}+2(\tfrac{g^v_x}{g^u_x})_x}{(g^u_x\cdot g^v_x)}((g^u)^{-1}(\zeta))
-\frac{1}{(g^u_x\cdot g^v_x)}((g^u)^{-1}(\zeta)) 1\!\mathrm{I}_{[g^v(0),\mathrm{b}(\zeta)]}(z),
\end{aligned}
\]

that is, to an integral equation of the form

\[
0= j(\zeta)-\mathrm{a}(\zeta)\,j(\mathrm{b}(\zeta))
+\int_{\underline{g}}^\zeta \widetilde{\mathrm{k}} (\zeta,z) j(z)\, dz
\]

for all $\zeta\in [g^u(0),g^u(L)]$.
Since this is a homogeneous second kind Volterra integral equation, using Gronwall's inequality and the fact that $0\leq\mathrm{a}\leq\kappa<1$, it follows that $j\equiv0$.

Thus we have

\[
\|\hat{f}^{+\,\prime\prime}\|_{C(J)}\leq C \|b\|_{C(\Omega)}\,.
\]

Moreover we can replace $C(J)$ by $C^\beta(J)$ since the operator $B$ defined above is an isomorphism also between $C^\beta(J)$ and $C^\beta(\Omega)$ as $Bj=b$ implies

\[
\begin{aligned}
&\frac{|j(g^u(x))-j(g^u(y))|}{|g^u(x)-g^u(y)|^\beta} \\
&= g^v_x(x))g^u_x(x)^{-1}\frac{|j(g^v(x))-j(g^v(y))|}{|g^u(x)-g^u(y)|^\beta}
+\frac{g^v_x(x))g^u_x(x)^{-1}-g^v_x(y))g^u_x(y)^{-1}}{|g^u(x)-g^u(y)|^\beta} j(g^v(y))\\
\end{aligned}
\]

hence, by taking the supremum over $x\in\Omega$ and using the range condition as well as the fact that by (<ref>), $|g^u(x)-g^u(y)|^\beta\geq |g^v(x)-g^v(y)|^\beta$,

\[
\begin{aligned}
&\leq \frac{1}{1-\kappa} \Bigl( C_1 \|j\|_{C(J)} + \|b\|_{C^{0,\beta}(\Omega)}\Bigr)
\leq \frac{1}{1-\kappa} \Bigl( \frac{C_1}{\mu\delta(1-\kappa)} \|b\|_{C(\Omega)} + \|b\|_{C^{0,\beta}(\Omega)}\Bigr)
\end{aligned}
\]

under the smoothness assumptions made on $g^u$, $g^v$.

Here in case $\beta=0$, $C^{0,0}(\Omega)$ could as well be replaced by $L^\infty(\Omega)$, so we have the choice of spaces

\[
X_f= C^{2,\beta}(J) \mbox{ for some }\beta\in[0,1]\mbox{ or }X_f=W^{2,\infty}(J)
\]

Now we bound $\hat{a}^+$ in

\begin{equation}\label{eqn:diffu1}
(\hat{a}^+ g^u_x)_x(x)+\hat{f}^+(g^u(x))=\check{u}(x,T) \quad x\in\Omega
\end{equation}

(cf. (<ref>)), in terms of $\check{u}$, $\check{v}$ and $\hat{f}^+$ where the latter is then bounded by means of (<ref>).
Here we first of all try to choose the space for $\hat{a}^+$ analogously to [17] as

\[
X_a = \{d\in H^1(\Omega)\cap L^\infty(\Omega)\,:\, d(L)=0\} \,,
\]

where we have assumed to know $a$ at one of the boundary points (we take the right hand one again) for otherwise (<ref>) would not uniquely determine $\hat{a}^+$.
The estimate can simply be carried out as follows.
After integration with respect to space and dividing by $g^u_x$ we get

\begin{equation}\label{eqn:idahatplus}
\begin{aligned}
\hat{a}^+(x)&=\frac{1}{g^u_x(x)}\int_x^L (\check{u}(\xi,T)-\hat{f}^+(g^u(\xi)))\, d\xi\\
\hat{a}^+_x(x)&=
-\frac{g^u_{xx}(x)}{(g^u_x(x))^2}\int_x^L (\check{u}(\xi,T)-\hat{f}^+(g^u(\xi)))\, d\xi%\\
\end{aligned}
\end{equation}

and therefore

\[
\begin{aligned}
\|\hat{a}^+\|_{L^\infty(\Omega)}&\leq\frac{1}{\mu}\bigl(\|\check{u}(\cdot,T)\|_{L^1(\Omega)}
\|\hat{a}^+_x\|_{L^2(\Omega)}&\leq
\frac{\sqrt{L}\|g^u_{xx}\|_{L^2(\Omega)}}{\mu^2}\bigl(\|\check{u}(\cdot,T)\|_{L^1(\Omega)}
&\quad +\frac{1}{\mu}\bigl(\|\check{u}(\cdot,T)|_{L^2(\Omega)}+\sqrt{L}\|\hat{f}^+\|_{L^\infty(J)}\bigr)
\end{aligned}
\]

This gives an obvious mismatch with the much higher norm of $\hat{f}$ (and consequently of $\check{u}$, $\check{v}$) in the estimate (<ref>) which from [16] and [18] we know to be needed, though.
Also, Lemma <ref> would not work when applied to the first in place of the second derivative of $\hat{f}^+$, as both $\fhatplusp $ terms in the left hand side of (<ref>) have the same factor.

Thus we next try a space for $a^+$ that is more aligned to estimate (<ref>).
From (<ref>) we get, using the fact that $C^{2,\beta}(\Omega)$ is a Banach algebra and $\|j(g)\|_{C^{2,\beta}(\Omega))}\leq \|j\|_{C^{2,\beta}(J))}\|g\|_{C^{2,1}(\Omega))}$,

\begin{equation}\label{eqn:esta_withf}
\begin{aligned}
\Bigl(L\|\bigl(\tfrac{1}{g^u_x}\bigr)_x\|_{C^{2,\beta}(\Omega)}+ \|\tfrac{1}{g^u_x}\|_{C^{2,\beta}(\Omega)}\Bigr) \bigl(\|\check{u}(T)\|_{C^{2,\beta}(\Omega)}+\|g^u\|_{C^{2,1}(\Omega))} \|\hat{f}^+\|_{C^{2,\beta}(J))}\bigr)
\end{aligned}
\end{equation}

and thus, together with (<ref>),

\begin{equation}\label{eqn:esta}
\begin{aligned}
\|\hat{a}^+\|_{C^{3,\beta}(\Omega)}\leq \bar{C}(g^u,g^v) \Bigl(\|\check{u}(T)\|_{C^{2,\beta}(\Omega)}+\|\check{v}(T)\|_{C^{2,\beta}(\Omega)}\Bigr)\,,
\end{aligned}
\end{equation}

where again $\bar{C}(g^u,g^v)$ only depends on the smoothness of $g^u,g^v$, as well as on and $\kappa, \delta, \mu$.
This completes step (a).

For $g^u\in C^{4,\beta}(\Omega)$, $g^v\in C^4(\Omega)$ satisfying (<ref>), (<ref>), (<ref>), (<ref>), there exists a constant $\bar{C}(g^u, g^v)$ depending only on
$\|g^u\|_{C^{4,\beta}(\Omega)}$, $\|g^v\|_{C^4(\Omega)}$, $\mu$, $\kappa$, $\delta$ such that the operator $\mathbb{M}:X_a\times X_f\to C^{2,\beta}(\Omega)$ as defined in (<ref>) with

\begin{equation}\label{eqn:XaXf}
\begin{aligned}
&X_a=\{d\in C^{3,\beta}(\Omega)\, : \, d(L)=0\}\,,\\
&X_f=\{j\in C^{2,\beta}(J)\, : \, j(g^u(0))=j(g^v(0))=j'(g^u(0))=j'(g^v(0))=0\}
\end{aligned}
\end{equation}

is bounded and invertible with

\[
\|\mathbb{M}^{-1}\|_{C^{2,\beta}(\Omega)\to X_a\times X_f}\leq \bar{C}(g^u, g^v)\,.
\]

This part of the convergence proof would directly carry over to the case of a fractional
time derivative $D_t^\alpha$ in place of $D_t$ in
(<ref>), (<ref>).
However, further bounding $\bar{u}=D^\alpha u$ using the fact that it satisfies
the pde
$D_t^\alpha \bar{u}-(a \bar{u}_x)_x=D_t^\alpha f(u)+ D_t^\alpha r^u =D_t^\alpha (\frac{f(u)}{u} \cdot u )+ D_t^\alpha r^u$ becomes problematic due to the lack of a chain or product rule for the fractional time derivative.

Extension to higher space dimensions does not appear to be possible, since the strategy of eliminating $\hat{a}^+$ relies on one dimensional integration.

Besides some assumptions directly on the the searched for coefficients $a$ and $f$ (see (<ref>)), the proof of Proposition (<ref>) shows that we also need to impose some conditions on the data $g^u$, $g^v$. The monotonicity assumption (<ref>) allows for inversion of $g^u$ in order to recover values of $f$ from values of $f(g^u)$, avoiding potentially contradictory double or multiple assignments; likewise for (<ref>) and $g^v$.
Also, it is clear that the range of the data must cover all values that will actually appear as arguments in $f$; note that the range condition (<ref>) is only imposed on the exact states. The range of the data is therefore the natural maximal domain for any reasonable reconstruction of $f$ and evaluation of $f$ outside this domain must be avoided since it could lead to false results in a forward simulation.
Finally, (<ref>) is a condition on sufficient deviation of slopes between $g^u$ and $g^v$, in order to allow for unique and stable recovery of $f^+$ according to Lemma (<ref>). It is probably the most technical and most difficult to realize of these condition. Still note that it is compatible with (<ref>) in the sense that the steeper function $g^u$ is the one with larger range. Note that no such condition as (<ref>) will be needed in the analysis of the final time - time trace observation setting of the next section.

All these conditions can in principle be achieved by appropriate excitations $r^u$, $r^v$ and their validity can be checked during the reconstruction process. Maximum and comparison principles can provide a good intuition for this choice. However, the concrete design of excitations to enable / optimize reconstruction is certainly a topic on its own.

Step (b), smallness of

The differences $\hat{u}$, $\hat{v}$ along with $\bar{u}=u_t$, $\bar{v}=v_t$ solve

\begin{eqnarray}
&&D_t\hat{u}-(a_{ex} \hat{u}_x)_x+q^u\,\hat{u} = (\hat{a} u_x)_x +\hat{f}(u) \quad t\in(0,T)\,, \quad \hat{u}(0)=0 \label{eqn:uhat}\\
&&D_t\hat{v}-(a_{ex} \hat{v}_x)_x+q^v\,\hat{v} = (\hat{a} v_x)_x +\hat{f}(v) \quad t\in(0,T)\,, \quad \hat{v}(0)=0 \label{eqn:vhat}
\end{eqnarray}

\[
\begin{aligned}
q^u&=-\frac{f_{ex}(u)-f_{ex}(u_{ex})}{u-u_{ex}}=-\int_0^1f_{ex}'(u_{ex}+\theta\hat{u})\, d\theta
\,, \quad
q^v&=-\frac{f_{ex}(v)-f_{ex}(v_{ex})}{v-v_{ex}}%=-\int_0^1f_{ex}'(v_{ex}+\theta\hat{v})\, d\theta
\end{aligned}
\]

\begin{eqnarray*}
&&\bar{u}_t-(a \bar{u}_x)_x-f'(u)\bar{u}=D_tr^u \quad t\in(0,T)\,, \quad
\bar{u}(0)=(a u_{0x})_x+f(u_0)+r^u(0)\\
&&\bar{v}_t-(a \bar{v}_x)_x-f'(v)\bar{v}=D_tr^v \quad t\in(0,T)\,,  \quad
\bar{v}(0)=(a v_{0x})_x+f(v_0)+r^v(0)\,.
\end{eqnarray*}

Moreover we can write $\check{u}=D_t\hat{u}$, $\check{v}=D_t\hat{v}$ as solutions to the pdes

\begin{eqnarray}
&&\check{u}_t-(a_{ex} \check{u}_x)_x-f_{ex}'(u_{ex})\,\check{u} = (\hat{a} \bar{u}_x)_x
+\Bigl(\hat{f}'(u) + \int_0^1f_{ex}''(u_{ex}+\theta\hat{u})\, d\theta\hat{u}\Bigr) \bar{u}
\,, \label{equ:ucheck}\\
&&\check{v}_t-(a_{ex} \check{v}_x)_x-f_{ex}'(v_{ex})\,\check{v} = (\hat{a} \bar{v}_x)_x
+\Bigl(\hat{f}'(v) + \int_0^1f_{ex}''(v_{ex}+\theta\hat{v})\, d\theta\hat{v}\Bigr) \bar{v}
\label{equ:vcheck}
\end{eqnarray}

with initial conditions

\begin{equation}\label{eqn:initucheckvcheck}
\check{u}(0)= (\hat{a} u_{0x})_x +\hat{f}(u_0)\,, \qquad
\check{v}(0)= (\hat{a} v_{0x})_x +\hat{f}(v_0)\,,
\end{equation}

where we have used the identity $f'(u)D_tu-f_{ex}'(u_{ex})D_tu_{ex}$
+ \Bigl((f-f_{ex})'(u)+(f_{ex}'(u)-f_{ex}'(u_{ex}))\Bigr)D_tu$.

We can estimate $\check{u}$, $\check{v}$ from (<ref>), (<ref>), (<ref>) in the same way as we estimate $z$ in the estimate preceding Theorem 3.4 of [18],
taking into account the additional term $(\hat{a} \bar{u}_x)_x$ in the right hand side and $(\hat{a} u_{0x})_x$ in the initial condition, which yields the following.

We split the solution of (<ref>) $\check{u}=\check{u}^r+\check{u}^0$ into a part $\check{u}^r$ satisfying the inhomogeneous pde

\[
\check{u}^r_t-(a_{ex} \check{u}^r_x)_x - {f_{ex}}'(u_{ex})\, \check{u}^r = (\hat{a} \bar{u}_x)_x
+\Bigl(\hat{f}'(u) + \int_0^1f_{ex}''(u_{ex}+\theta\hat{u})\, d\theta\hat{u}\Bigr) \bar{u}
\]

with homogeneous initial conditions $\check{u}^r(x,0)=0$
and a part $\check{u}^0$ satisfying the homogeneous pde
$\check{u}^0_t-(a_{ex} \check{u}^0_x)_x - {f_{ex}}'(u_{ex}) \check{u}^0 = 0$ with inhomogeneous initial conditions

\begin{equation}\label{eqn:checku0_init}
\check{u}^0 (x,0)=(\hat{a} u_{0x})_x(x) +\hat{f}(u_0(x)) \quad  x\in\Omega\,.
\end{equation}

Using a classical estimate from the book by Friedman <cit.>, and abbreviating $\Theta_t=(0,t)\times\Omega$, $\Theta=(0,T)\times\Omega$,
we estimate

\begin{equation}\label{eqn:estzr0}
\begin{aligned}
\leq
\sum_{|m|\leq2}\|D^m_x \check{u}^r\|_{C^{0,\beta}(\Theta_t)}\\
&\leq K\Bigl(\|(\hat{a} \bar{u}_x)_x\|_{C^{0,\beta}(\Theta_t)}+
\|\int_0^1 {f_{ex}}''(u_{ex}+\theta\hat{u})\, d\theta \, \hat{u} + \hat{f}'(u)\|_{C^{0,\beta}(\Theta)} \|\bar{u}\|_{C^{0,\beta}(\Theta_t)}\Bigr) \\
&\leq K\Bigl(\|\hat{a}\|_{C^{2,\beta}(\Theta_t)}+
\Bigr)
\|\bar{u}\|_{C^{2,\beta}(\Theta_t)}\,,
\end{aligned}
\end{equation}

\[
\begin{aligned}
\end{aligned}
\]

\[
\|\check{u}\|_{C^{2,\beta}(\Theta_t)}
\leq \|\check{u}^r\|_{C^{2,\beta}(\Theta_t)}+\|\check{u}^0\|_{C^{2,\beta}(\Theta_t)}
\leq \sum_{|m|\leq2}\|D^m_x \check{u}^r\|_{C^{0,\beta}(\Theta_t)}
\]

Note that we actually estimate $\|\check{u}(T)\|_{C^{2,\beta}(\Omega)}$ by $\|\check{u}\|_{C^{2,\beta}(\Omega\times(0,T))}:=\sum_{m\leq2}
\|D_x^k\check{u}\|_{C^{\beta}(\Omega\times(0,T))}$
with the definition of $\|\cdot\|_{C^{\beta}(\Omega\times(0,T))}$ from [9] (note that $\|\cdot\|_{C^{2,\beta}(\Omega\times(0,T))}$ has a different meaning in [9]).

Here, $\hat{u}$ solves (<ref>), so applying again <cit.> we obtain

\[
\begin{aligned}
\|\hat{u}\|_{C^{2,\beta}(\Theta_t)}&\leq
\sum_{|m|\leq2}\|D^m_x \hat{u}\|_{C^{0,\beta}(\Theta_t)}
\leq K \|(\hat{a}u_x)_x+\hat{f}(u)\|_{C^{0,\beta}(\Theta_t)} \\
K \Bigl(\|\hat{a}\|_{C^{2,\beta}(\Theta_t)}\|\hat{u}\|_{C^{2,\beta}(\Theta_t)}+\|\hat{f}\|_{C^{0,\beta}(J)} (1+\|u_{ex}\|_{C^{0,1}(\Theta_t)}+\|\hat{u}\|_{C^{0,1}(\Theta_t)})\Bigr)\\
K \Bigl(\|\hat{a}\|_{C^{2,\beta}(\Theta_t)}+\|\hat{f}\|_{C^{0,\beta}(J)}\Bigr) (1+\|u_{ex}\|_{C^{0,1}(\Theta_t)}+\|\hat{u}\|_{C^{2,\beta}(\Theta_t)})\Bigr),
\end{aligned}
\]

thus, for $\|\hat{a}\|_{C^{2,\beta}(\Theta_t)}+\|\hat{f}\|_{C^{0,\beta_0}(J)}\leq\rho_0<\frac{1}{K}$,

\[
\|\hat{u}\|_{C^{0,1}(\Theta_t)}
\leq \|\hat{u}\|_{C^{2,\beta}(\Theta_t)}
\leq \frac{1}{1-\rho_0 K}
\Bigl(\|\hat{a}\|_{C^{2,\beta}(\Theta_t)}+\|\hat{f}\|_{C^{0,\beta}(J)}\Bigr)
\]

Altogether, for $\|\hat{a}\|_{C^{2,\beta}(\Theta_t)}+\|\hat{f}'\|_{C^{0,\beta}(J)}\leq\rho_1$ with $\rho_1<\frac{1}{K}$, we end up with an estimate for $\check{u}^r$ of the form

\begin{equation}\label{eqn:estucheckr}
\begin{aligned}
\|\check{u}^r\|_{C([0,t];C^{2,\beta}(\Omega))}
\leq &C(K,\rho_0,\rho_1,\|{f_{ex}}''\|_{C^{0,\beta}(J)}) \Bigl(
\|\hat{a}\|_{C^{2,\beta}(\Omega)}
\end{aligned}
\end{equation}

To estimate $\check{u}^0$, we do not use the same result from <cit.>, since this would not give a contraction estimate with respect to $\hat{a}$, $\hat{f}$. Rather we attempt to employ dissipativity of the equation, however, this fails like in [18], as we will now illustrate.
Let us first of all point out that the difficulty here lies in the fact that the coefficient $f_{ex}'(u_{ex})$ is time-dependent and thus the abstract ODE corresponding to (<ref>) is non-autonomous.
Thus semigroup decay estimates would require the zero order coefficient to be
constant in time or at least time periodic,
see, e.g., <cit.>. To illustrate why indeed the general
form of $f_{ex}'(u_{ex})$ most probably prohibits decay of solutions, we follow
a perturbation approach, with a time constant potential
$q\approx -f_{ex}'(u_{ex})$, which we assume to be positive.
Using the series expansion of $\check{u}^0$ in terms of the eigenvalues
and -functions $(\lambda_n,\phi_n)_{n\in\mathbb{N}}$ of the elliptic operator
$\mathbb{L}$ defined by
$\mathbb{L}v=-(a_{ex}v_x)_x+q\,v$ as well as the induced Hilbert spaces

\[
\dot{H}^\sigma(\Omega)=\{v\in L^2(\Omega)\ : \ \sum_{n=1}^\infty \lambda_n^{\sigma/2} \langle v,\phi_n\rangle \phi_n \, \in L^2(\Omega)\}
\]

where $\langle\cdot,\cdot\rangle$ denotes the $L^2$ inner product on $\Omega$, with the norm

\|v\|_{\dot{H}^\sigma(\Omega)}
= \Bigl(\sum_{n=1}^\infty \lambda_n^\sigma \langle v,\phi_n\rangle^2 \Bigr)^{\!\frac{1}{2}}

that is equivalent to the $H^\sigma(\Omega)$ Sobolev norm, we can write

\begin{equation}\label{eqn:checku0}
\begin{aligned}
\check{u}^0(x,t)=&\sum_{n=1}^\infty \bigl(e^{-\lambda_n t} \langle (\hat{a} u_{0x})_x \!+\!\hat{f}(u_0),\phi_n\rangle
\!+\! \int_0^t e^{-\lambda_n (t-s)} \langle (f_{ex}'(u_{ex}(s))+q)\check{u}^0(s),\phi_n\rangle\, ds\bigr) \phi_n(x)\\
\end{aligned}
\end{equation}

From Sobolev's embedding theorem with $\sigma > d/2+2+\beta = 5/2 + \beta$,

\begin{equation}\label{eqn:estucheck01C2beta}
\begin{aligned}
\|\check{u}^{0,1}(t)\|_{C^{2,\beta}(\Omega)}&\leq C_{\dot{H}^\sigma,C^{2,\beta}}^\Omega
\Bigl(\sum_{n=1}^\infty \lambda_n^\sigma e^{-2\lambda_n t} \langle (\hat{a} u_{0x})_x+\hat{f}(u_0),\phi_n\rangle^2\Bigr)^{1/2}\\
&\leq  C_{\dot{H}^\sigma,C^{2,\beta}}^\Omega
\sup_{\lambda\geq\lambda_1}\lambda^{\sigma/2-1} e^{-\lambda t}
\|(\hat{a} u_{0x})_x+\hat{f}(u_0)\|_{\dot{H}^2(\Omega)}\\
&\leq  C_{\dot{H}^\sigma,C^{2,\beta}}^\Omega
\Psi(t;\sigma,\lambda_1)
\Bigl(\|\hat{a}\|_{C^2(\Omega)} \|u_{0\,xx}\|_{\dot{H}^2(\Omega)} \\
&\qquad +\|\hat{a}\|_{C^3(\Omega)} \|u_{0\,x}\|_{\dot{H}^2(\Omega)}
+ C_\mathbb{L}\|\hat{f}\|_{C^2(J)} \|u_0\|_{\dot{H}^2(\Omega)}\Bigr)
\end{aligned}
\end{equation}

\begin{equation}\label{eqn:estucheck02C2beta}
\begin{aligned}
\|&\check{u}^{0,2}(t)\|_{C^{2,\beta}(\Omega)}\leq C_{\dot{H}^\sigma,C^{2,\beta}}^\Omega
\Bigl(\int_0^\tau\sum_{n=1}^\infty \lambda_n^\sigma e^{-2\lambda_n (t-s)}
\langle ({f_{ex}}'(u_{ex}(s))+q)\check{u}^0(s),\phi_n\rangle^2\, ds\Bigr)^{1/2}\\
&\leq  C_{\dot{H}^\sigma,C^{2,\beta}}^\Omega C(\Omega)
\left(\int_0^t\Psi(\tau-s;\sigma,\lambda_1)\|{f_{ex}}'(u_{ex}(s))+q\|_{H^2(\Omega)}^2\, ds \right)^{1/2}
\|\check{u}^0\|_{C([0,t];C^2(\Omega))}\,,
\end{aligned}
\end{equation}

with $C_{\mathbb{L}}$ such that
$\|\mathbb{L} j(v)\|_{L^2(\Omega)}\leq C_{\mathbb{L}} \|j\|_{C^2(\mathbb{R})} \|v\|_{\dot{H}^2(\Omega)}$ for all $j\in C^2(\mathbb{R})$, $v\in \dot{H}^2(\Omega)$ and

\begin{equation}\label{eqn:Psi}
\Psi(t;\sigma,\lambda_1) = \begin{cases}
(\sigma/2-1)^{\sigma/2-1} e^{1-\sigma/2}\, t^{1-\sigma/2}&\mbox{ for }t\leq \frac{\sigma-2}{2\lambda_1}\\
\lambda_1^{\sigma/2-1} e^{-\lambda_1 t}&\mbox{ for }t\geq \frac{\sigma-2}{2\lambda_1}
\,.\end{cases}
\end{equation}

The appearance of $\|\check{u}^0\|_{C([0,t];C^2(\Omega))}$ on the right hand side of (<ref>) forces us to look at the supremum of $\|\check{u}^0(t)\|_{C^2(\Omega)}$ over $t\in[0,T]$, which, however, due to the singularity of $\Psi$ at $t=0$ cannot be estimated appropriately by (<ref>).
Thus the convolution term in (<ref>) inhibits exponential decay as would be expected from the estimate of the first term in (<ref>).

Thus we have no means of establishing contractivity in the presence of the
initial term (<ref>) and therefore will assume that
$u_0$ and $f(0)$ vanish.
Thus, putting this altogether

\begin{equation}\label{eqn:estucheck}
\begin{aligned}
\|\check{u}^r\|_{C([0,t];C^{2,\beta}(\Omega))}\\
&\leq C
\|u_{ex,t}\|_{C^{2,\beta}(\Theta)}
\,.
\end{aligned}
\end{equation}

provided $\|\hat{a}\|_{C^{2,\beta}(\Theta_t)}+\|\hat{f}\|_{C^{1,\beta}(J)}\leq\rho$ small enough.

The same estimate can be used for bounding $\|\check{v}\|_{C^{2,\beta}(\Omega\times(0,T))}$

This together with (<ref>), (<ref>) yields contractivity of $\mathbb{T}$ on a ball of sufficiently small radius $\rho$.

Under the assumptions of Proposition <ref> and if additionally
$\|u_{ex,t}\|_{C^{0,\beta}(\Theta)}$ is sufficiently small and $u_0=0$, $f_{ex}(0)=0$, $g^u(0)=0$, there exists $\rho>0$ such that $\mathbb{T}$ is a self-mapping on $B_\rho^{X_a\times X_f}(a_{ex},f_{ex})$ and the convergence estimate

\[
\|\mathbb{T}(a,f)-\mathbb{T}(a_{ex},f_{ex})\|_{X_a\times X_f}\leq
q \|(a,f)-(a_{ex},f_{ex})\|_{X_a\times X_f}
\]

holds for some $q\in(0,1)$ and $X_a$, $X_f$ as in (<ref>).

In <cit.> we have
alternatively proven contractivity for monotone $f$.
However, the approach taken there does not seem to go through here
for the following reasons.
Additionally to the exponential decay of $\bar{u}$,
we would need to show exponential decay of $(\hat{a}\, \bar{u}_x)_x$.
Another obstacle is that the strategy of showing contractivity of $\mathbb{T}$ by exploiting the regularity gain

\[
\begin{aligned}
\|\mathbb{T}^2(a,f)-\mathbb{T}^2(a_{ex},f_{ex})\|_{C^{3,\beta}(\Omega)\times C^{2,\beta}(J)}
&\leq C(g^u,g^v)\|(\check{u},\check{v})\|_{C^{2,\beta}(\Omega)^2}\\
&\leq \tilde{C}\,C(g^u,g^v)\|\mathbb{T}(a,f)-\mathbb{T}(a_{ex},f_{ex})\|_{C^{1,\beta}(\Omega)\times C^{1,\beta}(J)}
\end{aligned}
\]

does not work here.
This is because we have no possibility to estimate the $C^{1,\beta}(J)$
norm of the $f$ part of
$\|\mathbb{T}(a,f)-\mathbb{T}(a_{ex},f_{ex})\|_{C^{1,\beta}(\Omega)\times C^{1,\beta}(J)}$ in terms of the $C^{1,\beta}(\Omega)$ norm of $\check{u},\check{v}$. (So far the Volterra integral equation approach from step (a) only works for differentiability order $\geq2$.)
However, the maximal parabolic regularity approach in [18] relies on the embedding $W^{2,p}(\Omega)\to C^{1,\beta}(\Omega)$ applied to $\check{u},\check{v}$ and therefore only works with the $C^{1,\beta}$ norm of $\hat{f}$.

§.§ Final time / time trace observations of one or two states

Identify $a(x)$, $f(u)$ in (<ref>), (<ref>)
with homogeneous impedance boundary conditions (<ref>),
from observations

\begin{equation}\label{eqn:mixeddata}
g(x)=u(x,T), \quad x\in\Omega, \qquad h(t)=v(\tilde{x},t), \quad t\in\Theta\subseteq[0,T]\,,
\end{equation}

for some $\tilde{x}\in\overline{\Omega}$, typically on the boundary. We here assume that $\tilde{x}$ is the right hand boundary point $L$; other cases can be treated analogously.
However, if $\tilde{x}$ is an interior point, then some term that we can eliminate by using the boundary condition would have to be additionally be taken into account in the analysis.
More precisely, we make use of the right hand boundary condition that implies
$a^+_x(L) v_x(L,t;a,f) =- \frac{a^+_x(L)}{a^+(L)} \gamma^v v(L,t;a,f)$ and replace $v(L,t;a,f)$ by the data $h(t)$ here.
Thus, again projecting onto the measurement manifolds $\Omega\times\{T\}$ and $\{L\}\times(0,T)$,
we define the fixed point operator $\mathbb{T}$ by $(a^+,f^+)=\mathbb{T}(a,f)$ with

\begin{equation}\label{eqn:fp_fiti-titr}
\begin{aligned}
&(a^+ g_x)_x(x)+f^+(g(x))=D_tu(x,T;a,f)-r^u(x,T) \quad x\in\Omega\\
&a^+(L) v_{xx}(L,t;a,f) - \frac{a^+_x(L)}{a^+(L)} \gamma^v h(t)+f^+(h(t))=h_t(t)-r^v(L,t) \quad t\in\Theta\,,
\end{aligned}
\end{equation}

which, when assuming that $a^L=a_{ex}(L)$ is known and $\gamma^v(L)=0$, simplifies to

\begin{eqnarray}
(a^+ g_x)_x(x)+f^+(g(x))&=&D_tu(x,T;a,f)-r^u(x,T) \quad x\in\Omega \label{eqn:faplus_mixed_u}\\
f^+(h(t))&=&h_t(t)-r^v(L,t)-a^L v_{xx}(L,t;a,f) \quad t\in\Theta\,. \label{eqn:faplus_mixed_v}
\end{eqnarray}

The recursion for the the differences $\hat{a}^+=a^+-a_{ex}$, $\hat{f}^+=f^+-f_{ex}$, $\hat{a}=a-a_{ex}$,
reads as

\begin{eqnarray}
(\hat{a}^+ g_x)_x(x)+\hat{f}^+(g(x))&=&D_t\hat{u}(x,T) \quad x\in\Omega\label{eqn:ahat_mixed}\\
\hat{f}^+(h(t))&=&-a^L \hat{v}_{xx}(L,t;a,f) \quad t\in\Theta\,.\label{eqn:fhat_mixed}
\end{eqnarray}

where $\hat{u}$, $\hat{v}$ solve (<ref>), (<ref>), that is

\[
\mathbb{M}(\hat{a}^+,\hat{f}^+)=\mathbb{F}(a,f)-\mathbb{F}(a_{ex},f_{ex})\,,
\]

with the operators $\mathbb{M}$ and $\mathbb{F}$ defined by

\[
\mathbb{M}(a^+,f^+)=
\left(\begin{array}{c}(a^+ g_x)_x+f^+(g^u)\\f^+(h)\end{array}\right)\,, \quad
\mathbb{F}(a,f)=
\left(\begin{array}{c}u_t(\cdot,T;a,f)\\ -a^Lv_{xx}(L,\cdot;a,f)\end{array}\right)\,.
\]

Inverting $\mathbb{M}$ is therefore much easier this time. We can first invert (<ref>) for $f^+$ relying on the range condition

\begin{equation}\label{eqn:rangecondition_mixed}
J=h(\Theta)\,, \quad
J\supseteq u_{ex}(\Omega\times(0,T))\,, \quad J\supseteq v_{ex}(\Omega\times(0,T))
\end{equation}

(in place of (<ref>)) and then compute $a^+$ from (<ref>) in exactly the same way as we have done before, cf. (<ref>). It will turn out that a lower regularity estimate of $a^+$ and $f^+$ (or more precisely, of $\hat{a}^+$ and $\hat{f}^+$, in order to prove contractivity) suffices, that is, we will use the function spaces
$X_a=C^{2,\beta}(\Omega)$, $X_f=C^{1,\beta}(J)$ (with further boundary conditions).

Step (a), bounding $\mathbb{M}^{-1}$:
From (<ref>), we get, after differentiation,

\[
\|\fhatplusp (h)\, D_th\|_{C^{0,\beta}(\Theta)}\leq |a^L| \|D_t\hat{v}_{xx}\|_{C^{0,\beta}(\Omega\times(0,T))}
\]

thus, assuming strict monotonicity of $h$, that is,

\begin{equation}\label{eqn:mu_h}
\end{equation}

we get, using $\hat{f}^+(p)=0$ for some $p\in J$,

\begin{equation}\label{eqn:estf_mixed}
\begin{aligned}
\|\hat{f}^+\|_{C^{1,\beta}(J)}\leq C(h) |a^L| \|D_t\hat{v}\|_{C^{2,\beta}(\Omega\times(0,T))}
\end{aligned}
\end{equation}

where $C(h)$ only depends in $\mu$ and $\|h\|_{C^{1,1}(\Theta)}$.
Likewise, we could attempt to estimate a higher order norm of $f^+$ by differentiating a second time

\[
\|\fhatpluspp(h)\, (D_th)^2\|_{C^{0,\beta}(\Theta)}\leq
\|\fhatplusp (h)\, D_t^2h\|_{C^{0,\beta}(\Theta)}
+|a^L| \|D_t^2\hat{v}_{xx}\|_{C^{0,\beta}(\Omega\times(0,T))}\,.
\]

However, this would involve too high derivatives of $\hat{v}$.

From (<ref>), that is (<ref>), we get
when taking only the $C^{1,\beta}$ norm,
(i.e., the $C^{2,\beta}$ norm of $\hat{a}$)

\begin{equation}\label{eqn:esta_withf_mixed}
\begin{aligned}
\bigl(L\|\bigl(\tfrac{1}{g_x}\bigr)_x\|_{C^{1,\beta}(\Omega)}+ \|\tfrac{1}{g_x}\|_{C^{1,\beta}(\Omega)}\bigr) \bigl(\|D_t\hat{u}(T)\|_{C^{1,\beta}(\Omega)}+\|g\|_{C^{1,1}(\Omega))} \|\hat{f}^+\|_{C^{1,\beta}(J)}\bigr)
\end{aligned}
\end{equation}

in place of (<ref>), thus, together with (<ref>), using again the abbreviations $\check{u}=D_t\hat{u}$, $\check{v}=D_t\hat{v}$,

\begin{equation}\label{eqn:esta_mixed}
\begin{aligned}
\|\hat{a}^+\|_{C^{2,\beta}(\Omega)}\leq
\bar{C}(g,h)\Bigl(\|\check{u}(T)\|_{C^{1,\beta}(\Omega)}+
\|\check{v}\|_{C^{2,\beta}(\Omega\times(0,T))}\Bigr)
\end{aligned}
\end{equation}

that is an estimate in norms that are by one order lower than those in (<ref>).

For $g\in C^{3,\beta}(\Omega)$, $h\in C^{1,1}(\Theta)$ satisfying (<ref>), (<ref>), (<ref>), $p\in J$,
there exists a constant $\bar{C}(g^u, g^v)$ depending only on
$\|g\|_{C^{3,\beta}(\Omega)}$, $\|h\|_{C^{1,1}(\Theta)}$, $\mu$, $\delta$, such that the operator $\mathbb{M}:X_a\times X_f\to C^{1,\beta}(\Omega)$ as defined in (<ref>) with

\begin{equation*}
\begin{aligned}
&X_a=\{d\in C^{2,\beta}(\Omega)\, : \, d(L)=0\}\,,\\
&X_f=\{j\in C^{1,\beta}(J)\, : \, j(p)=0\}
\end{aligned}
\end{equation*}

is boundedly invertible with

\[
\|\mathbb{M}^{-1}\|_{C^{1,\beta}(\Omega)\to X_a\times X_f}\leq \bar{C}(g^u, g^v)\,.
\]

The obvious spatially higher dimensional version of (<ref>), (<ref>) is

\begin{eqnarray}
\nabla\cdot(\hat{a}^+ \nabla g)(x)+\hat{f}^+(g(x))&=&D_t\hat{u}(x,T;a,f) \quad x\in\Omega \label{eqn:ahat_mixed1}\\
\hat{f}^+(h(t))&=&-a^L \triangle\hat{v}(\tilde{x},t;a,f) \quad t\in\Theta\,.\label{eqn:fhat_mixed1}
\end{eqnarray}

However, (<ref>), which is a transport equation for $\hat{a}^+$, in general only yields an estimate of $\hat{a}^+$ in a weaker norm than the one needed here.

step (b) smallness of $\mathbb{F}(a,f)-\mathbb{F}(a_{ex},f_{ex})$:

For $D_t\hat{u}=\check{u}$ we can just use the estimate (<ref>),

and again assume that $u_0$ and $f(0)$ vanish.

For $g\in C^{3,\beta}(\Omega)$, $h\in C^{2,\beta}(\theta)$ satisfying (<ref>), (<ref>), (<ref>),  and if additionally $f_{ex}(0)=0$, $u_0=0$,
$\|u_{ex,t}\|_{C^{2,\beta}(\Theta)}$ is sufficiently small, there exists $\rho>0$ such that $\mathbb{T}$ is a self-mapping on $B_\rho^{X_a\times X_f}(a_{ex},f_{ex})$ and the convergence estimate

\[
\|\mathbb{T}(a,f)-\mathbb{T}(a_{ex},f_{ex})\|_{X_a\times X_f}\leq
q \|(a,f)-(a_{ex},f_{ex})\|_{X_a\times X_f}
\]

holds for some $q\in(0,1)$ and

\begin{equation}\label{eqn:XaXf_mixed}
\begin{aligned}
&X_a=\{d\in C^{2,\beta}(\Omega)\, : \, d(L)=0\}\,,\\
&X_f=\{j\in C^{1,\beta}(J)\, : \, j(0)=0\}
\end{aligned}
\end{equation}

As opposed to the previous section, here the strategy of proving contractivity for monotone $f$ from [18] would also work here since we can estimate the $C^{1,\beta}$ norm of $\hat{f}^+$. This would also allow to tackle nonhomogeneous initial data by using a maximal parabolic regularity estimate, see <cit.>.

Since no such condition as (<ref>)
(that is sufficiently different slopes of $g_u$, $g_v$) is needed here, the result from Theorem th:contr_fiti_titr remains valid with a single run, i.e., $v=u$.

§.§ Time trace observations of two states

A well-known “folk theorem” suggests that trying to
reconstruct a function depending on a variable $x$
from data measured in an orthogonal direction to $x$ is
inevitably severely ill-conditioned.
This has been borne out in innumerable cases.
The recovery of $a(x)$ in the parabolic equation
$u_t - \nabla.(a\nabla u) = f(u)$ from time-valued data
is an example of the above.
It is known that recovery of spatially-dependent coefficients in an
elliptic operator from time data-valued measurements made on the boundary
$\partial\Omega$ is extremely ill-conditioned with an exponential dependence
of the unknown in terms of the data.
We give a synopsis of the reasons for this that have particular relevance
to our current theme; while recovery of the diffusion coefficient $a(x)$ from
final time data is only a mildly ill-conditioned problem, that from
time-trace data is extremely ill-conditioned.

We will set the problem in one space dimension taking
$\Omega=[0,1]$ and remove all non-essential terms.
Thus we consider the recovery of the uniformly positive  $a(x)$ from
$u_t - (au_x)_x = 0$ with (say) boundary conditions
$u_x(0,t) = 0$, $u(1,t) = 0$ with initial value $u_0(x)$
where our overposed data is the value of $h(t) = u(0,t)$.
Note that this is a single function recovery from a linear problem.
The standard approach here is due to Pierce [24].
The solution to the direct problem has the representation
$u(x,t) = \sum_{n=1}^\infty \langle u_0,\phi_n\rangle\,e^{-\lambda_nt}\phi_n(x)$
where $\{\lambda_n,\phi_n\}$ are the eigenvalues / eigenfunctions of
$\mathbb{L} = -(a u_x)_x$ under the conditions $u_x(0)=u(1)=0$
and we normalise these by
the standard Sturm-Liouville condition $\phi_n'(0)=1$.
For a given $a(x)$ this is transmutable into the norming constants
and more usual Hilbert space norm
[28, 29].
We have to recover $a(x)$ from $h(t) = u(0,t)$, that is from

\begin{equation}\label{eqn:Dirichlet_series}
h(t) = \sum_{n=1}^\infty A_n e^{-\lambda_n t}
\qquad
A_n = \langle u_0,\phi_n\rangle.
\end{equation}

Note that $\lambda_n>0$ and that $h(t)$ is in fact analytic.
It is well-known that the pair $\{\lambda_n,A_n\}$ is uniquely recoverable
from $h(t)$ and this is easily seen by taking Laplace transform
(analytically extending $h(t)$ to all $t>0$ if need be).
Then (<ref>) becomes

\begin{equation}\label{eqn:rat_func_series}
\mathcal{L}(h(t))\to s := H(s) = \sum_{n=1}^\infty \frac{A_n}{s+\lambda_n}
\end{equation}

and shows that knowing the meromorphic function $H(s)$ provides the location
$\lambda_n$ and residues $A_n$ of its poles.
Of course, the extreme ill-conditioning of this step is now apparent;
we know the values of the analytic function $H(s)$ on the positive real line
and have to recover information on the negative real axis – in fact
all the way to $-\infty$ as the eigenvalues have the asymptotic behaviour
$-\lambda_n \sim -c n^2$.
Recovery of just the eigenvalues isn't sufficient to recover $a(x)$
but we have further information, namely the sequence $\{A_n\}$.
The simplest way to use this is to select $u_0$ to be impulsive – say
$u_0(x) = \delta(x)$ for then $A_n = \phi_n(0)$.
Thus with this we have both the spectrum and an endpoint condition
on the eigenfunctions and standard results from inverse Sturm-Liouville theory
show this is sufficient for the recovery of $a(x)$,
[4, 11, 28].

Actually, with this approach the computational cost of an effective
reconstruction of $a(x)$ is relatively low.
Determining the coefficients from a finite sum in (<ref>)
is just Padé approximation of $H(s)$ and there are efficient methods
such as in [3] to effect this recovery.
The reconstruction of $a(x)$ from this finite spectral can be accomplished
extremely quickly, [28].

An entirely similar argument holds if instead our purpose was  to recover
the potential $q(x)$ in $\,\mathbb{L} = -u_{xx} + q(x)u = 0$.
The recovery of two spectral sequences proceeds identically and
this potential form for $\mathbb{L} $  is the canonical version
[29, 6] of the inverse Sturm-Liouville problem
to which all others can be transformed by means of the Liouville transform.
In fact, this part of the process is only mildly ill-conditioned in terms
of mappings between spaces for we have
$\|q_1 - q_2\|_{L_2} \leq C\|\lambda_{n,1}-\lambda_{n,2}\|_{\ell^2}$.

However, the above paints a rather too rosy picture of the difficulties
The basic asymptotic form for the eigenvalues for a $q\in L^2$ is

\begin{equation}\label{eqn:sturm_liouville_eigen}
\lambda_n(q) = \bigl(n-{\textstyle{\frac{1}{2}}}\bigr)^2\pi^2 +
\int_0^1 q(s)\,ds + \epsilon_n, \qquad \epsilon_n \in \ell^2
\end{equation}

The dominant term in this expansion is identical for all $q\in L^2$.
The more regularity assumed of $q(x)$ the more rapid the decay of
the “information sequence” $\{\epsilon_n\}$,
[27, 6, 30].
For example, if $q\in H^{m,2}$ then
$\epsilon_n \sim \sigma_n/n^{2[m/2]+1}$ where again $\sigma_n\ell^2$.
Thus the approximation to the eigenvalues  obtained from inverting
(<ref>) will inevitably contain errors and
from this we have to recover the sequence $\{\epsilon_n\}$ by first
subtracting off the term asymptotic to effectively $n^2\pi^2$.
For $n=10$ this is the order of a thousand and meantime the information
is in a possibly rapidly decreasing sequence of the difference.
In short, the inverse Sturm-Liouville problem is very mildly ill-conditioned from a space domain/range definition. but the strong masking by the
fixed leading term in the eigenvalue expansion which is independent of $q$
makes recovery of anything more than a few Fourier modes impractical from
other than exceedingly accurate spectral data.

From the existence of the Liouville transform one might conclude that
the picture is identical for the determination of the coefficient $a(x)$.
This is indeed so up to the point where we recover the spectral information.
However, the version of (<ref>) for this case is

\begin{equation}\label{eqn:sturm_liouville_eigen_a}
\lambda_n(q) = {\textstyle \frac{(n-\frac{1}{2})^2\pi^2}{L^2}} + \epsilon_n,
\qquad\mbox{where }\
L = \int_0^1 [a(s)]^{-1/2}\,ds.
\end{equation}

Thus unlike the potential case where the leading and dominant term is
independent of $q(x)$, the leading term now contains information about the
coefficient $a(x)$.
For a fixed error in the data, this can make a substantial difference.

All of the above has been predicated on two factors; that we are in one
space dimension and our underlying parabolic equation is linear.
With a nonlinear reaction term $f(u)$ the whole approach above fails
(and this is precisely our case of interest).
In higher space dimensions there is no sequel to the inverse Sturm-Liouville
theory and much more information is required over the one dimensional case.

There is an approach that can transcend these issues.
This is to look directly at the maps from unknown coefficient to the data:
$F_a[a] \to u(0,t)$ or $F_q[q] \to u(0,t)$.
We can get a sense of the practical invertibility by looking at the
Jacobian matrix $J$ arising from a finite dimensional version of
$F'$, where $F := [F_a\; F_q]$, taken in a direction
$\delta a$ or $\delta q$:

\begin{equation}\label{eqn:F'_map}
\begin{aligned}
&\hat u_t - \nabla\cdot(a(x)\nabla \hat u) = f'(u)\hat u -
\nabla.(\delta a\nabla u) \\
&\hat u_t - \triangle \hat u + q(x) \hat{u}= f'(u)\hat u - \delta q\,u \\
\end{aligned}
\end{equation}

An important question is whether this is feasible in respect to even the
theoretical invertibility of the map $F'$:
either as one of its two components or in combination.
We provide a sketch below to show that in the linear case $f(u) = c(x)u$
$F'$ is indeed locally invertible in $\mathbb{R}^d$ for $d>1$,
then later consider the issue of its
practical invertibility by estimating the singular values
of the matrix $J$ representing $F'$.
With this $f$ we can absorb the contribution into the zero'th order term of
$\mathbb{L}$ and use the same notation for the eigenvalues/eigenfunctions
of $\mathbb{L}$ on $\Omega$ with respect to homogeneous impedance boundary
With this, the  solution to (<ref>) in the potential case is given by

\begin{equation}\label{eqn:F'_map_pot}
\hat u(x,t) = \sum_1^\infty \int_0^t e^{-\lambda_n(t-\tau)}
\langle \delta q\, u_0,\phi_n\rangle e^{-\lambda_1\tau}\phi_n(x).
\end{equation}

Suppose now we have $\hat u(\tilde{x},t) = 0$ for $x=\tilde{x}\in\partial\Omega$
(if we only have this prescribed on $[0,T]$ then we simply extend
it by zero for $t>T$) and taking Laplace transforms gives

\begin{equation}\label{eqn:F'_map_zero}
\sum_1^\infty \frac{1}{p+\lambda_n}\frac{\phi_n(\tilde{x})}{p+\lambda_1}
\langle \delta q\; u_0,\phi_n\rangle = 0
\end{equation}

for all $p>0$.
Now we can select the point $\tilde{x}$ so that $\phi_n(\tilde{x}) \not=0$
since $\phi_n(x)$ satisfies homogeneous impedance conditions on
$\partial\Omega$ (or more exactly can select the origin of the eigenfunctions
to accomplish this for a prescribed $\tilde{x}$).
Then analyticity in (<ref>) shows that
$\langle \delta q\, u_0,\phi_n\rangle = 0$ for all $n$ and completeness
implies that $\delta q \,u_0 =0$ and, by our assumption that $u_0$ was
nonzero, implies that $\delta q = 0$.
An entirely analogous argument holds for the case of $a(x)$.
It also trivially  extends to an initial value $u_0$ that is
a finite combinations of eigenfunctions and with a little more work
extends to a general initial value $u_0$ that is not
identically zero in any subset of $\Omega$ of positive measure.
Finally, by considering the situation on the interval $[\epsilon,T]$
we have the new initial value $u(x,\epsilon)$ which is strictly positive
from the maximum principle provided only that $u_0(x)$ is non-negative
and non-trivial.
Thus we have proven,

Suppose the nonlinear term $f(u)$ is zero, the initial condition
$u_0(x)\not=0$ is nonnegative in $\Omega$, then the maps $F':q(x) \to u(\tilde{x},t)$
and $F':a(x) \to  u(\tilde{x},t)$ are injective.

We have carried out the process of computing the singular values of $F'$
taking the case of one spatial dimension and about a constant value for
$a$ or $q$, namely $a=1$ and $q=0$.
The initial condition was the first eigenfunction for
the associated elliptic operator, here $\sin(\frac{\pi}{2}x)$ and
the basis functions used for both $\delta a$ and $\delta q$ were
$\{\sin(n\pi x)\}_{n=1}^{20}$.
The resulting singular values  are shown in the Figure fig:sv_a-q
plotted on a $\log_{10}$ scale.

Logarithms of the singular values for $a(x)$ and $q(x)$
recovery from time trace data using $F:\delta c \to u(1,t)$ with $\{\delta c =\sin(n\pi x)\}$

The expected  exponential decay of these singular values is cleary evident.
From this data $h(t)$ – which arose from a problem with the initial data
being the exact lowest eigenfunction and hence the best possible for
reconstructions, we see that only a few of the lowest singular values
are useable under anything but exceedingly high accuracy data.
These correspond to the lowest frequencies in the basis set.
While the same exponential decay is true for both $a(x)$ and $q(x)$, the rate
is greater for $q$, but more crucially, the first singular values in each
case differ considerably - by roughly a factor of ten.
Note also the effect that smaller time steps makes.
It is much less evident in the case shown  due to the solution $u$
to the base problem having only the lowest spectral mode.
In problems where higher modes are present and one has
to recover $\lambda_n$ from $e^{-\lambda_n t}$; the smaller the time step
$dt$ the better one can recover the eigenvalues and coefficients
(especially the largest index ones sought)
from data error effecting this term.
For example, for $n=10$ and Dirichlet boundary conditions with $a=1$ and $q=0$
this means the factor $e^{-100\pi^2 dt}$ should lie above the noise level.
Of course,  if the diffusion coefficient is substantially smaller
then this changes the above estimate.

Our scaling of this problem where we have assumed the coefficients are of
order unity must be taken into account.
As noted in the introduction, in many physical settings
the diffusion coefficient could be several orders of magnitude smaller
and this would give a corresponding change in the spectrum $\{\lambda_n\}$.
In addition, for diffusion in fluids there is often a significant
dependence of the diffusion coefficient on temperature $u$ in addition
to any spatial dependence.
For example a standard model here is due to Arrhenius
[2, 7]
that supposes that $a = a_0\,e^{-E/R\,u}$ where $R$ is the
universal gas constant and $E$ the activation energy which is typically
such that $E >> R\,u$.
Of course, incorporating this into our model would dramatically increase
the complexity of any analysis.

It is entirely possible to extend these computation to higher
space dimensions and also in fact to include nonlinear reaction terms
although in the latter situation the above theory and hence the proof of
Theorem thm:q-a-injectivity will no longer be valid.

This same effect, namely the superior reconstruction of $a(x)$ against
$q(x)$  from their simultaneous recovery from final time data,
was noticed in [17].
There is another “folk theorem” at play here.
The dependence of the solution of an elliptic equation on changes in its
coefficients is far greater from those appearing in higher order terms
(as in the $a(x)$ here) than in terms appearing in lower order coefficients.
One should expect that this translates into data terms arising from
projections of the solution onto a surface,
and in turn expect that one should be able to
more easily reconstruct coefficients appearing in the highest derivatives of
the operator as distinct from those appearing in the lowest order term.
from knowledge of projections of the solution onto a surface.

The reconstruction of both $a(x)$ and $f(u)$ requires two experiments
providing solutions and the corresponding data obtained from them.
The theoretical results of section <ref> shows this can
be sufficient under appropriate conditions.
In many applications one could attempt an oversampling by taking measurements
from $m>2$ data runs generating $m$ solutions providing overposed data.
Under the unrealistic situation of no data error this would simply lead
to redundancy in terms of uniqueness but might very well lead to
superior reconstructions due to the fact that exact optimality conditions
are difficult to determine and more measurements would likely lead to a
better case being obtained in one of the experiments that in turn would lead
to a larger number of singular values above a given threshold.

In the case of noise in the data one would also expect that oversampling
by utilizing $N$ solutions would lead to an effective lower noise rate
and in consequence superior reconstructions.
In the case of $N=2$ and 1% random noise there is certainly a detectable
variation in the corresponding reconstructions and even a modest increase
in the number of measurements makes a significant reduction.

The extent to which this can be quantified depends on several factors.
One of these is the underlying source of the noise and its probability
density distribution.  If this is Gaussian then $N$ samples would reduce
the effective error by a factor of $\sqrt{N}$, but this holds only
for large $N$ and other distributions would give different answers.
However we are dealing with a far from linear process.
The reconstruction process is itself highly nonlinear in addition to the basic
forwards map being so due to the nonlinear term $f(u)$.
As a result of this, quantification of the overall uncertainty is far
from an easy question.

§ RECONSTRUCTIONS

We will show the results of numerical experiments to recover both
$a$ and $f$  with two different types of data measurements.

The first of these is when we are only able to obtain “census-type”
information and thus measure $g_u := u(x,T)$ and $g_v := v(x,T)$ for some
fixed time $T$ and as the result of two experiments obtained by
altering the initial, boundary or forcing term conditions.

The second is when we are able to measure multiple types of data from a
single experimental run: specifically both the final data and a
boundary measurement of either the solution $u$ or its normal derivative at a
fixed point $\tilde  x\in\partial\Omega$.

To procure data for the reconstructions a direct solver based on a
Crank-Nicolson scheme produced output values and data
values were produced from these by sampling at a relatively small
number $N_x$ and $N_t$ of points in both the spatial and temporal directions:
This sampled data was then both interpolated to a full working size to obtain
data values $g(x)$ and $h(t)$ commensurate to the grid being used by the solver
used in the inverse problem and such that the value of $g(x)$ was
smoothed by an $H^2$ filter while that of $h(t)$ by an $H^1$ filter.

All the reconstructed solutions we show are set in one space dimension.
This is in part to make the graphical results more transparent but also because
some of the algorithms have an $\mathbb{R}^1$ restriction from a
computational aspect and certainly from an analysis one as noted in the
previous section.
However, we will make note when extensions to higher space dimensions are feasible.

§.§ Final time with two runs

We shall pose the setting in one space variable as the exposition is simpler
as is the graphical representation of the reconstructions obtained.

The problem is to identify $a(x)$, $f(u)$ in

\begin{equation*}%\label{eqn:u-v_pair}
\begin{aligned}
u_t-(a u_x)_x &=f(u) + r_u(x,t,u)\\
v_t-(a v_x)_x &=f(v) + r_v(x,t,v)\\
\end{aligned}
\end{equation*}

where $r_u$ and $r_v$ are known driving terms and where $u(x,t)$ and $v(x,t)$
are subject to homogeneous impedance boundary conditions

a\frac{\partial u}{\partial x} +\alpha u =0
\qquad
a\frac{\partial v}{\partial x} +\beta v =0

and initial values

u(x,0)=u_0(x) \qquad
v(x,0)=v_0(x) \qquad

from observations $g_u(x)=u(x,T)$, $g_v=v(x,T)$.

Some restrictions must be imposed on this data.
An ideal situation will include the fact $\nabla g_u$ and $\nabla g_v$
will not vanish on the interior of $\Omega$.
This can be accomplished by varying the impedance coefficients
over the boundary $\partial\Omega$ and by imposing a larger positive flux into
the interior where $\alpha$, $\beta$ are smallest - for example Neumann

Based on the fixed point operator defined in (<ref>), we now
derive some reconstruction schemes based on different ways of resolving
(<ref>) in an efficient way,
and provide computational results obtained with these schemes.

The first method is one of sequential iteration: recovering first
a new approximation to $a(x)$ using one direct solve
then followed by an update of $f(u)$ using a second direct solve,
but now using the new value of $a$.
More precisely,
let $f^0$, $a^0$ to some initial approximation, then for $k=0,1,2,\ldots$

compute $u_t(x,T;a^k,f^k)$ by solving (<ref>)
using the first set of initial/boundary conditions then differentiate
$u$ with respect to time at $t=T$ and for all $x\in\Omega$.

Set $\phi(x):= u_t(x,T;a^k,f^k) - f^k\bigl(g_u(x)\bigr) - r_u(x,T,g_u)$
so that $\bigl((a(x) g_u'(x)\bigr)' = \phi(x)$.
Integrate this equation over $(0,x)$ using the boundary conditions
inherited by $g_u$ and $u(x,T;a^k,f_k)$  to obtain
$\Phi(x) := \int_0^x [u_t(s,T) - f^k(g_u(s)) - r_u(s,T)]\, ds$
which then gives the update

\begin{equation}\label{eqn:iter1}
\end{equation}

Use this new approximation $a^{k+1}$ for $a(x)$ together with the previous
value $f^k$
for $f$ to compute $v(x,T;a^{k+1},f^k)$ from the imposed second set of
initial/boundary conditions and thereby obtain the value of $v_t(x,T)$.

We now update $f$ by

f^{k+1}(g_v)= v_t(x,T)-\bigl(a^{k+1} g'_v(x)\bigr)'

One drawback of this method is the need to differentiate the just-updated
value of $a(x)$ that can lead to numerical instabilities.
Here is a way to avoid this
and works well in a single space dimension.

We again let $f^0$, $a^0$ be some initial approximation, then for

Compute both $u_t(T)$ and $v_t(T)$ by solving
(<ref>), using the two sets of initial/boundary values

\begin{equation}\label{eqn:u-v_pair}
\begin{aligned}
u_t-(a u_x)_x &=f(u) + r_u(x,t,u)\qquad u(x,0)=u_0,\quad
a\frac{\partial u}{\partial x} +\alpha u = k_u(x,t)\\
v_t-(a v_x)_x &=f(v) + r_v(x,t,v)\qquad v(x,0) = v_0\quad
a\frac{\partial u}{\partial x} +\beta v = k_v(x,t)\\
\end{aligned}
\end{equation}

and hence obtain $u_t(x,T;a^k,f^k)$ and $v_t(x,;a^k,f^kT)$.

Let $W(x) = g_u(x) g'_v(x) - g_v(x) g'_u(x)$ and
$W_p(x) = g_u(x) g''_v(x) - g_v(x) g''_u(x)$.
Multiply the first equation in (<ref>) by $g_u$,
the second by $g_v$ then subtract obtaining

\bigl(a(x)g_u'(x)\bigr)'g_v(x) - \bigl(a(x)g_v'(x)\bigr)'g_u(x) = \psi(x)

\begin{equation*}
\begin{aligned}
\psi(x) &= [u_t(x,T;a^k,f^k) - r_u(x,T,g_u) - f(g_u(x))]g_v(x) \\
&\quad -[v_t(x,T;a^k,f^k) - r_v(x,T,g_v) - f(g_v(x))]g_u(x)
\end{aligned}
\end{equation*}

Now integrating over $(0,x)$ we obtain

\begin{equation}\label{eqn:aW}
a(x)W(x) = \int_0^x \psi(s)\,ds
\end{equation}

from which an update $a^{k+1}(x)$ can be obtained by division by $W$.
It is obviously an advantage if $W$ doesn't vanish and this condition can
be obtained by controlling the values of $\alpha$, $\beta$ and
the input flux values $k_u$ and $k_v$.
(Even if $W$ has isolated zeros there are several approaches that can be
taken to still resolve (<ref>) in a stable manner.)

To recover the next approximation for $f$ we multiply the first equation in
(<ref>) by $g'_v$, the second by $g'_u$ then subtract obtaining

\tilde\phi = (u_t - r_u(x,T)g_v' - (v_t - r_v(x,T))g_u' - a(x)W_p(x)

Then $f^{k+1}$ is determined by

\begin{equation}\label{eqn:f-update}
f^{k+1}\bigl(g_u(x)\bigr)g_v'(x) - f^{k+1}\bigl(g_v(x))\bigr)g'_u(x) =
\tilde\phi(x).
\end{equation}

The effect of this is to eliminate the derivative of the just computed $a(x)$.

Solving (<ref>) is not completely straightforward.
It is a delayed argument equation for $f$ if the maximum
value achieved for both $u$ and $v$ occurs at or near the same
point on $\partial\Omega$ but in fact the resolution of the entire scheme
is better if both $u$ and $v$ are “more independent.”

One approach is to use a basis representation for $f$ to  resolve this
Another, which works well if one of $g'_u$ or $g'_v$ is monotone,
is to use the approximation $\tilde f_1(u) = f^k(u)$ and then iterate
successively on (<ref>) by

\begin{equation}\label{eqn:f_successive-update}
\tilde f_{j+1}\bigl(g_u(x)\bigr)g'_v(x) =
\tilde f_j\bigl(g_v(x))\bigr)g'_u(x) + {\tilde\phi}(x).
\end{equation}

This method worked well over a wide variety of examples of functions $f(u)$
and one such pairing of $a(x)$ and $f(u)$ recoveries is shown in
Figure fig:fiti_fiti below.

Recovery of $a(x)$ and $f(u)$ from final data

The initial approximations were $a^0(x)=1$ and $f^0(u)=0$.
An important point to observe is the speed of convergence of the
scheme updating $a(x)$.  The second iteration was already extremely good
and the difference between the third and tenth iteration shown would
not be distinguishable at the figure resolution.
On the other hand, the iterations for $f$ proceeded much more slowly.

One factor that should be noted here is that we used both data sets to
recover $a$ initially, and then similarly for $f$.
Much the same picture would have resulted from reversing this order.
The fact is that the diffusion coefficient $a(x)$ dominates the
equation in terms of its ability to modify solutions $u(x,t:a;f)$ in
comparison to that for $f(u)$.
This situation was true whether pointwise or basis reconstructions schemes
were used for representing the unknowns.

In higher space dimensions it is possible from a reconstruction perspective
to measure less data. For example, in $\mathbb{R}^d$ with $d=2,\,3$ one
data run for $g_u(x) = u(x,T)$ measures this quantity for all $x\in\Omega$, but
the second run measures only $g_v(x) = v(s(x),T)$ for some
one-dimensional curve $C$ parametrized by $s(x)$ that connects points
$x_1,\,x_2$ on $\partial\Omega$.
The strategy is to use the first of these to recover $a(x)$ then since $f(u)$
depends only the single variable $u$ we seek to recover $f(s(x))$ from this
second measurement.
The points $x_1$ and $x_2$ cannot be chosen arbitrarily but must in fact
correspond to the maximum and minimum values of $v(\cdot,T)$ on
$\partial\Omega$ in order to satisfy the range condition necessary
for recovering $f$.
We note here that we have  no theoretical analysis for this case;
neither a uniquness result nor a convergence theorem.

§.§ Final time and time-trace data

As noted previously, an alternative data measurement set is to measure both
boundary values as well as final time information, namely

\begin{equation}\label{eqn:fitititr_data}
u(x,T) = g_u(x) \qquad u(\tilde x,t) = h_u(t)
\end{equation}

where we have assumed that $\tilde x \in\partial\Omega$ is such that that
the maximum range of $u(x,t)$ occurs at that point.
It will also be convenient to take the impedance value to be zero
at $\tilde x$ and arrange that $g_u(x)$ is monotonic in $\Omega$.

We have chosen a particular choice of a time trace measurements.
Other possibilities include a time trace of $u(x^\star,t)$ for
some interior point $x^\star\in\Omega$.
The difficulty here is in ensuring that the range condition is satisfied
and this is likely to be prohibited by the maximum principle.
An alternative is a measurement of the energy
$E(t) = \int_\Omega |\nabla u|^2\,dx$ or some similar functional.
Again, a determining factor would be the range condition.

Again, starting from the fixed point iteration defined by
(<ref>), we derive a numerical reconstruction algorithm.
There are several ways to proceed and we start with pointwise updates
of both $a(x)$ and $f(u)$ in sequence at each full iteration step.

We again let $f^0$, $a^0$ be some initial approximation, then for

Compute $u(x,t;a^k,f^k)$.
Update $a(x)$ using the final time information $g_u(x)$ by setting

\begin{equation}\label{eqn:a_successive-update_2}
\begin{aligned}
\phi(x) &= u_t(x,T;a^k,f^k) - r_u(x,T,gu) - f^k\bigl(g_u(x)\bigr)
\qquad \Phi(x) = \int_0^x \phi(s)\,ds \\
a^{k+1}(x) &= \Phi(x)/g_u'(x) \\
\end{aligned}
\end{equation}

Now update $f(u)$ by projecting onto the boundary at $x=\tilde x$ using the
data $h_u(t)$ and computing $\triangle u(\tilde x,t) = u_{xx}(\tilde x,t)$
so that

\begin{equation}
\psi(t) = h_u(t) - r_u\bigl(\tilde x,t,h_u(t)\bigr)
- a^{k+1}(\tilde x) u_{xx}(\tilde x,t).
\end{equation}

Here the advantage of the imposed boundary condition $u_x(\tilde x,t) = 0$
is evident as there is no need to involve the derivative of $a$.

Now if $h_u(t)$ is monotone then $f = f^{k+1}(u)$ is recovered from

\begin{equation}
f^{k+1}\bigl(h_u(t)\bigr) = \psi(t).
\end{equation}

In Figure fig:fiti_titr_no_noise we show a pair of reconstructions of
$a(x)$ and $f(u)$ using the above scheme.
The spatial interval was taken as $0,L$ with $L=1$ and the final time
of measurement was $T=0.5$.
Thus data consisted of measurements of  $g(x_i) = u(x_i,T)$ and
$h(t_j)=u(L,t_j)$ where the numbers of sampled points were
$N_x=20$ and $N_t=25$.

In this reconstruction the initial values were once again
taken as $a^0 = 1$ and $f^0=0$.
In most applications a better initial approximation would likely be available.

As can be seen from this figure, the convergence of the scheme was rather
slow initially  and in fact the first approximation to the diffusion
coefficient $a(x)$ and the reaction $f(u)$ were both uniformly lower than
the actual values.
One might expect that a uniform too small value of $a$ would result
in an $f$ that was larger then the actual and this did occur on the second
iteration for $f$.  The reason it did not occur at the first iterative step
is at each step the updates for $a$ and $f$ are independent:
the former only using $g(x)$ and the latter only using $h(t)$.

This in fact indicates an almost limiting distance apart of
the initial approximations and the actual values of $a$ and $f$.
If an iterate of $a$ becomes too close to a zero value the scheme becomes
very unstable.
Also, in the presence of noise, one would in general have to take
a much closer initial approximation.

The reason for the final small amount of mis-fit of the reconstructed
and the actual functions is the necessity to use a smoothing scheme on the
data to restore the correct mapping properties.

Recovery of $\,a(x)\,$ and $\,f(u)\,$ from space-and-time data

Recovery of $\,a(x),$ and $\,f(u)\,$ from space-and-time data: 1% noise

Notice here the considerable degradation of the reconstructions under
even modest amounts of noise.
This is not an artifact of the method but inherent in the problem.
Reconstructions using a pair of final values as in the last section
are more forgiving under data noise but still show the limitations
of reconstructions in these inverse problems.

§.§ Non-trivial initial values

In section sec:conv we required that the initial value $u_0(x) = u(x,0) = 0$.
This was to remove a possible singularity arising from the convolution
in equation (<ref>).
This issue has been noted in previous work attempting to recover
a reaction term $f(u)$ from time-valued data on the boundary $\partial\Omega$,
[26, 25, 18]
even when this term is the only unknown in the model.
The question thus arises if this is a fundamental requirement or
merely a delicate technical issue to be surmounted.

Suppose that $u$ satisfies

u_t - (au_x)_x = f(u) + r(x,t),\qquad u(x,0) = u_0(x),\quad
{\mathbb B} u = \phi

where ${\mathbb B}$ is the boundary operator on $\partial\Omega$,
For continuity we should  impose the condition ${\mathbb B} u_0 = \phi$
at $t=0$ but, in particular we suppose $u_0$ is positive and has support
contained in the strict interior of $(0,1)$ to remove all boundary effects.

Now let $v(x,t) = u(x,t) - u_0(x)$ so that $v$ satisfies

v_t - (av_x)_x = \tilde f(v) + s(x,t),\qquad v(x,0) = 0 ,\quad {\mathbb B} v = 0

where $\tilde f(v) = f(u) = f(v+u_0)$ and $s(x,t) = r(x,t) - (a u_0')'$.

The overposed data $u(1,t)=h(t)$ is converted into $v=\tilde h(t)$ but
relating $\tilde h$ to $h$ follows from the fact that $u_0$ vanishes
at $x=0$ and $x=1$ gives that  $\tilde h(t) = h(t)$.

Thus in principle, the problem with nontrivial $u_0$ can be converted into
one with $u_0=0$ by the above transformation:
we get a new, known, right hand side driving term
(assuming $a(x)$ is known) and reconstruct as before to get
$\tilde f$.  Then from this, knowing $u_0$, recover $f$.
Does this then amount to resolution of the issue?
Not quite.  The previous convolution situation is not resolved by this
as the argument of $f$ in $\tilde f$ is now shifted.

However, the transformation does point to a definite issue for nontrivial $u_0$.
We have to ensure that the standard range condition is satisfied:
$h(t)$ must encompass all values that the domain of $f$ requires.
The given data $h(t)$ is the same for both $u$ and $v$ but
$\tilde f$ requires a different range than $f$ (since $u_0>0$)
and thus appears to have missing information.

In Figure fig:fiti_titr_u0 below we show final reconstructions
(taken as the $10^{\hbox{\sevenrm th}}$ iteration which certainly corresponded
to effective numerical convergence) for both $a(x)$ and $f(u)$.
This was under no data-noise conditions and sampling at a large number of
A nontrivial value of $u_0$ was taken; namely $u_0(x) = \beta\,x^2(1-x^2)$.
for the values $\,\beta = 1,\,2,\,5,\,10,\,20$.
This function and its first derivative is zero at both endpoints
so ${\mathbb B}u_0 = 0$. The value $\beta=20$ corresponds to a maximum value of
$1.25$ for $u_0$; $\beta=1$ gives a height of $0.0625$.
The rightmost graphic shows the reconstructed $f(u)$ for these choices.
In the case $\beta =0$ the reconstruction would have been indistinguishable
from the actual $f$ as the data was noise-free and sampled at a large number
of points.
Notice how rapidly the reconstructed $f$ deteriorates  from the actual
with increasing $\beta$ and hence size of $u_0$, as the
range condition violation,
which predominantly affected the smaller values of $u$, became stronger.
For sufficiently large $u_0$ this finally affects the entire reconstruction.

Recovery of $\,a(x)\,$ and $\,f(u)\,$ from space-and-time data when $u_0\not=0$

The coefficient $a(x)$ which is reconstructed solely from the final time
values $g(x)=u(x,T)$ is relatively immune to this effect -
until the reconstructed $f$ becomes sufficiently poor that this has an
effect on subsequent iterations of the combined scheme.
Once again, the smaller values of the data $g(x)$ used are most sensitive
to a nontrivial $u_0$ and since our $g(x)$ function is monotonically
increasing away from $x=0$, the values of $a(x)$ nearer this point are
most affected.
The leftmost graphic shows only values for $\beta=5$, $\beta=10$ and $\beta=20$ since
for smaller $\beta$ the reconstructed $a$ is indistinguishable from the actual
$a(x)$ at the resolution of the figure.
However, this does also illustrate the point that the diffusion term
$\nabla.(a\nabla u) $ plays a dominant role in the behaviour of the equation.

The experiment was repeated but with $u_0$ taken to be a narrow Gaussian
of unit fixed height and center $x_0$.
A similar effect to the above was observed:
a greater disparity between the reconstructions and the actual $f$
occurred as $x_0$ approached the boundary $x=1$ where the time-trace data
$h(t) = u(1,t)$ was measured.
This is in keeping with fact that the initial value is perturbing the
range condition and this is more pronounced the closer the support of the
perturbation is to the measurement boundary.

In short, nontrivial initial values used either directly or through the
lifting device can be used to extract information on both $a$ and $f$,
but this has a definite limitation which becomes more acute as $u_0(x)$
increases in magnitude and so more significantly affects the range condition.

§ ACKNOWLEDGMENTS

The work of the first author was supported by the Austrian Science Fund fwf
under the grant P30054.

The work of the second author was supported in part by the
National Science Foundation through award dms-1620138.

Moreover, the authors would like to thank both reviewers for their careful
reading of the manuscript and their valuable comments and suggestions that
have led to an improved version of the paper.

[1] (MR511740) [10.1016/0001-8708(78)90130-5]
D. G. Aronson and H. F. Weinberger,
Multidimensional nonlinear diffusion arising in population genetics,
Advances in Mathematics, 30 (1978), 33–76.

[2]
R. A. Arrhenius,
über die dissociationswärme und den einflußder temperatur auf den dissociationsgrad der elektrolyte,
Z. Phys. Chem., 4 (1889), 96–116.

[3] (MR1383091) [10.1017/CBO9780511530074]
G. A. Baker and P. Graves-Morris,
Padé Approximants,
Cambridge University Press, Cambridge, second edition, 1996.

[4] (MR15185) [10.1007/BF02421600]
G. Borg,
Eine umkehrung der Sturm-Liouville eigenwertaufgabe,
Acta Mathematica, 76 (1946), 1–96.

[5] [10.1002/9781118788295.ch4]
J. W. Cahn and J. E. Hilliard,
Free energy of a nonuniform system. I. Interfacial free energy,
The Journal of Chemical Physics, 28 (1958), 258–267.

[6] (MR1445771) [10.1137/1.9780898719710]
K. Chadan, D. Colton, L. Päivärinta and W. Rundell,
An Introduction to Inverse Scattering and Inverse Spectral Problems,
SIAM Monographs on Mathematical Modeling and Computation, SIAM, 1997.

[7]
E. L. Cussler,
Diffusion: Mass Transfer in Fluid Systems,
Cambridge University Press., 1997.

[8] (MR855754) [10.1007/BF00251803]
C. M. Elliott and Z. Songmu,
On the cahn-hilliard equation,
Arch. Rational Mech. Anal, 96 (1986), 339–357.

[9] (MR0181836)
A. Friedman,
Partial Differential Equations of Parabolic Type,
Prentice-Hall, 1964.

[10] (MR0445088)
A. Friedman,
Partial Differential Equations,
Holt, Rinehart and Winston, New York, 1969.

[11] (MR0073805) [10.1090/trans2/001/11]
I. M. Gel'fand and B. M. Levitan,
On the determination of a differential equation from its spectral function,
Amer. Math. Soc. Transl., 1 (1951), 253–291.

[12] (MR1423804)
P. Grindrod,
The Theory and Applications of Reaction-Diffusion Equations: Patterns and Waves,
2$^{nd}$ edition, Oxford Applied Mathematics and Computing Science Series, The Clarendon Press, Oxford University Press, New York, 1996.

[13] (MR1085828) [10.1002/cpa.3160440203]
V. Isakov,
Inverse parabolic problems with the final overdetermination,
Comm. Pure Appl. Math., 44 (1991), 185–209.

[14] (MR2193218)
V. Isakov,
Inverse problems for partial differential equations,
2$^{nd}$ edition, Applied Mathematical Sciences, 127, Springer, New York, 2006.

[15] (MR3975371) [10.1088/1361-6420/ab109e]
B. Kaltenbacher and W. Rundell,
On an inverse potential problem for a fractional reaction-diffusion equation,
Inverse Problems, 35 (2019), 065004.

[16] (MR4019539) [10.1088/1361-6420/ab2aab]
B. Kaltenbacher and W. Rundell,
On the identification of a nonlinear term in a reaction-diffusion equation,
Inverse Problems, 35 (2019), 115007.

[17] (MR4002155) [10.1016/j.jmaa.2019.123475]
B. Kaltenbacher and W. Rundell,
Recovery of multiple coefficients in a reaction-diffusion equation,
J. Math. Anal. Appl., 481 (2019), 123475.

[18]
B. Kaltenbacher and W. Rundell,
The inverse problem of reconstructing reaction-diffusion systems,
Inverse Problems, 36 (2020), 065011.

[19] (MR3839508)
C. Kuttler,
Reaction-diffusion equations and their application on bacterial communication,
in Disease Modelling and Public Health, Handbook of Statistics, Elsevier Science, 2017.

[20] (MR1056055) [10.1137/1032046]
H. A. Levine,
The role of critical exponents in blowup theorems,
SIAM Review, 32 (1990), 262–288.

[21] (MR3012216)
A. Lunardi,
Analytic Semigroups and Optimal Regularity in Parabolic Problems,
Modern Birkhäuser Classics, Springer Basel, 1995.

[22]
E. A. Mason,
Gas, State of Matter,
Encyclopedia Britannica, 2020.

[23] (MR1908418)
J. D. Murray,
Mathematical Biology I,
3$^{rd}$ edition, Interdisciplinary Applied Mathematics, 17, Springer-Verlag, New York, 2002.

[24] (MR534419) [10.1137/0317035]
A. Pierce,
Unique identification of eigenvalues and coefficients in a parabolic problem,
SIAM J. Control Optim., 17 (1979), 494–499.

[25] (MR1108124) [10.1002/num.1690030404]
M. Pilant and W. Rundell,
Iteration schemes for unknown coefficient problems in parabolic equations,
Numer. Methods for P.D.E., 3 (1987), 313–325, 1987.

[26] (MR829324) [10.1080/03605308608820430]
M. S. Pilant and W. Rundell,
An inverse problem for a nonlinear parabolic equation,
Comm. Partial Differential Equations, 11 (1986), 445–457.

[27] (MR894477)
J. Pöschel and E. Trubowitz,
Inverse spectral theory,
in Pure and Applied Mathematics, 130, Academic Press, Inc., Boston, MA, 1987.

[28] (MR1166492) [10.1088/0266-5611/8/3/007]
W. Rundell and P. E. Sacks,
The reconstruction of Sturm-Liouville operators,
Inverse Problems, 8(1992), 457–482.

[29] (MR1106979) [10.1090/S0025-5718-1992-1106979-0]
W. Rundell and P. E. Sacks,
Reconstruction techniques for classical inverse Sturm-Liouville problems,
Math. Comp., 58(1992), 161–183.

[30] (MR2311614) [10.1007/s11006-006-0204-6]
A. M. Savchuk and A. A. Shkalikov,
On the eigenvalues of the Sturm-Liouville operator with potentials in Sobolev spaces,
Mat. Zametki, 80 (2006), 864–884.

[31] (MR3671483) [10.1093/imamat/hxx004]
Z. Zhang and Z. Zhou,
Recovering the potential term in a fractional diffusion equation,
IMA J. Appl. Math., 82 (2017), 579–600.

Received March 2020; revised May 2020.

E-mail address<EMAIL_ADDRESS>
E-mail address<EMAIL_ADDRESS>