Recovery of multiple coefficients in a reaction-diffusion equation

TL;DR

This paper addresses the inverse problem of simultaneously recovering spatially-dependent coefficients in a reaction-diffusion equation, demonstrating uniqueness and convergence of an iterative method, including fractional and nonlinear cases.

Contribution

It introduces a novel iterative scheme for recovering multiple coefficients in reaction-diffusion equations, extending to fractional derivatives and nonlinear potentials.

Findings

Proves uniqueness of coefficient recovery under given conditions.

Establishes convergence of the proposed iterative scheme.

Handles fractional derivatives and nonlinear potential terms.

Abstract

This paper considers the inverse problem of recovering both the unknown, spatially-dependent conductivity $a(x)$ and the potential $q(x)$ in a parabolic equation from overposed data consisting of the value of solution profiles taken at a later time $T$. We show both uniqueness results and the convergence of an iteration scheme designed to recover these coefficients. We also allow a more general setting, in particular when the usual time derivative is replaced by one of fractional order and when the potential term is coupled with a known nonlinearity $f$ of the form $q(x)f(u)$.

Tables

Table 1. Table 1: Norm differences for versions of the iteration

Iteration	1	2	3	4	5	6
Parallel scheme
${\Vert a_n-a_{\mathrm{act}}\Vert }_{\infty }$	0.0046	0.0071	0.0052	0.0042	0.0038	-
${\Vert q_n-q_{\mathrm{act}}\Vert }_{\infty }$	0.7723	0.2845	0.3368	0.3561	0.3618	-
${\Vert a_n-a_{\mathrm{act}}\Vert }_2$	0.0039	0.0093	0.0055	0.0038	0.0031	-
${\Vert q_n-q_{\mathrm{act}}\Vert }_2$	0.7372	0.3109	0.1433	0.0980	0.0904	-
Eliminate $q$ scheme
${\Vert a_n-a_{\mathrm{act}}\Vert }_{\infty }$	0.0084	0.0099	0.0076	0.0065	0.0062	0.0061
${\Vert q_n-q_{\mathrm{act}}\Vert }_{\infty }$	0.8080	0.2957	0.3066	0.3258	0.3321	0.3343
${\Vert a_n-a_{\mathrm{act}}\Vert }_2$	0.0053	0.0092	0.0051	0.0036	0.0032	0.0030
${\Vert q_n-q_{\mathrm{act}}\Vert }_2$	0.8510	0.3320	0.1442	0.0961	0.0892	0.0886
Eliminate $a$ scheme
${\Vert a_n-a_{\mathrm{act}}\Vert }_{\infty }$	0.1254	0.0541	0.0518	0.0520	0.0521	-
${\Vert q_n-q_{\mathrm{act}}\Vert }_{\infty }$	1.3512	1.7705	1.8898	1.9318	1.9462	-
${\Vert a_n-a_{\mathrm{act}}\Vert }_2$	0.0347	0.0151	0.0123	0.0119	0.0118	-
${\Vert q_n-q_{\mathrm{act}}\Vert }_2$	0.9042	0.7562	0.7626	0.7788	0.7860	-

Equations

$$u_t(x,t)+\mathbb{L}u(x,t)=r(x,t,u)$$

$$u_t(x,t)-\nabla \cdot (a(x)\nabla u(x,t))=q(x)f(u)+r(x,t,u)$$

$$a\frac{\partial u}{\partial \nu }+\gamma u=h,x\in \partial \Omega ,t>0,$$

$$u(x,0)=u_0(x),x\in \Omega .$$

$$D_t^{\alpha }u(x,t)+\mathbb{L}u(x,t)=r(x,t,u)$$

$$u(x,t)=\sum _{n=1}^{\infty }⟨u_0,{\varphi }_n⟩E_{\alpha ,1}(-{\lambda }_n{t}^{\alpha }){\varphi }_n(x)$$

$$D_t^{\alpha }u-\nabla \cdot (a\nabla u)+qf(u)=r_{u}t\in (0,T),u(0)=u_0$$

$$D_t^{\alpha }v-\nabla \cdot (a\nabla v)+qf(v)=r_{v}t\in (0,T),v(0)=v_0$$

$$a{\partial }_{\nu }u+\gamma u=s(x,t),x\in \partial \Omega u(x,0)=u_0,v(x,0)=v_0$$

$$-\nabla \cdot (a\nabla g_u)+qf(g_u)$$

$$-\nabla \cdot (a\nabla g_v)+qf(g_v)$$

$$\mathbb{M}\left[\genfrac{}{}{0ex}{}{a}{q}\right]=\mathbb{F}(a,q)$$

$$-\nabla \cdot (a^{+}\nabla g_u)+q^{+}f(g_u)=r_u(T)-D_t^{\alpha }u(\cdot ,T;a,q)$$

$$-\nabla \cdot (a^{+}\nabla g_v)+q^{+}f(g_v)=r_v(T)-D_t^{\alpha }v(\cdot ,T;a,q)$$

$$\nabla \cdot (aW)=a\nabla \cdot W+\nabla a\cdot W=\varphi$$

$$W(x)$$

$$\varphi (x)$$

$$a(x)W(x)={\int }_0^x\varphi (s)ds$$

$$q(x)W(x)$$

$$\psi (x)$$

$$-(a(x)u_x{)}_x+q(x)u=\lambda r(x)u0<x<1$$

$$-v_{yy}+Q(y)v=\mu v0<y<L$$

$$\ell (x)={\int }_0^x[r(s)/a(s){]}^{1/2}ds$$

$$Q(y)=[r(x)]^{-1}q(x)-[a(x)/r^{3}(x)]^{1/4}\bigl{[}a(x)\bigl{(}[a(x)r(x)]^{-1/4}\bigr{)}^{\prime}\bigr{]}^{\prime}$$

$$-{\hat{\varphi }}_n^{\prime \prime }+Q_1(x){\hat{\varphi }}_n={\lambda }_n{\hat{\varphi }}_n,-{\hat{\psi }}_n^{\prime \prime }+Q_2(x){\hat{\psi }}_n={\mu }_n{\hat{\psi }}_n,$$

$$|{\lambda }_n-{\mu }_n|\le C\Vert Q_1-Q_2{\Vert }_2$$

$${\lambda }_n=n^2{\pi }^2+{\int }_0^{1}Q(t)dt-{\int }_0^{1}Q(t)\cos 2n\pi tdt+t_n,\{t_n\}\in {\ell }^2$$

$$v(x)=u(x)+{\int }_0^{x}K(x,t)u(t)dt$$

$$\displaystyle K_{tt}-K_{xx}+\bigl{(}Q(x)-P(x)\bigr{)}K$$

$$K(x,\pm x)$$

$$K(x,\pm x)=h+{\frac{1}{2}}\int_{0}^{x}\bigl{(}Q(s)-P(s)\bigr{)}\,ds$$

$${\hat{\varphi }}_n(x)=\frac{\sin (\sqrt{{\lambda }_n}x)}{\sqrt{{\lambda }_n}}+{\int }_0^{x}K(x,t)\frac{\sin (\sqrt{{\lambda }_n}t)}{\sqrt{{\lambda }_n}}dt.$$

$$K(x,t)+{\int }_{R}Q(r)K(r,s)drds={\int }_{\frac{x-t}{2}}^{\frac{x+t}{2}}Q(s)ds$$

Full text

Keywords: Inverse problem, undetermined coefficients, diffusion equation

1 Introduction

We consider the inverse problem of recovering coefficients from the uniformly elliptic operator $-\mathbb{L}$ within a diffusion model. In this case we will include both parabolic as well as anomalous diffusion processes and the situation we describe will be general enough to include known nonlinear reaction terms. Reaction-diffusion equations such as these occur throughout the sciences and we give some specific examples in the next section. Let $\,\mathbb{L}u=-\nabla\cdot(a(x)\nabla u)+q(x)u$ be defined on a domain $\Omega\subset\mathbb{R}^{n}$ with smooth boundary $\partial\Omega$ and where the two coefficients $a$ and $q$ are the quantities to be determined in the inverse problem. In this setting our basic model equation is thus

[TABLE]

where $r(x,t,u)$ is a known forcing function. Extensions of the above are of course possible and we mention the case of a reaction-diffusion model in which (1) becomes

[TABLE]

where the form of the nonlinear driving term $f(u)$ is assumed known. Boundary conditions for (1) will be of the impedance form

[TABLE]

and we impose the initial condition

[TABLE]

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, [9, 19]. Now of course the recovery of the coefficients $a$ and $q$ requires over-posed data and we shall assume this is a spatial measurement at a fixed time $T$ for two different sets of boundary conditions or the value of a single solution at two different later times $t=T_{1},\;T_{2}$. Under different assumptions on the continuous time random walk ctrw model one obtains alternative diffusion processes and we consider the subdiffusion model based on fractional time derivatives. Now the basic equations take the form

[TABLE]

where $D_{t}^{\alpha}$ denotes the Djrbashian-Caputo derivative of order $\alpha$. Setting all of the above in context, we want to also understand how the different diffusive processes effect the ability to recover the coefficients in the inverse problem.

2 Background

Undetermined coefficient problems based on the equation (5) have a long history in the literature. In general, the over-posed data has taken one of two types: the data suggested above, namely the spatial values of $u(x,T)$ for fixed time $T$; or time trace data, typically measured at discrete points $\{x_{i}\}$ on the lateral boundary of the cylinder $\Omega\times(0,t)$ for $t>0$.

In the latter case if the boundary conditions in (3) are of impedance type with $\gamma<\infty$ then this is taken to be Dirichlet values and if $\gamma=\infty$, that is Dirichlet conditions are imposed then the over-posed data is flux values at $\{x_{i}\}$. The latter situation has been the most common, in particular in one spatial dimension, beginning with the work of Cannon and of Pierce, [2, 20] and continuing in the fractional diffusion case by, for example, [4, 23]. The techniques used have mostly revolved around the eigenfunction expansion of the solution (in the homogeneous case)

[TABLE]

where $\{\lambda_{n},\,\phi_{n}(x)\}$ are the eigenvalue/eigenfunctions of $-\mathbb{L}$ on $\Omega$ and $E_{\alpha,\beta}$ is the Mittag-Leffler function. When $\alpha=1$ this recovers the usual exponential function leading to the familiar parabolic solution. This representation is based on the Djrbashian-Caputo derivative from the initial point $a$ $\,{}^{C}_{a}D^{\alpha}_{t}f=I_{a}^{\alpha}\frac{df}{ds}\,$ where $\;I_{a}^{\alpha}f(x)=\frac{1}{\Gamma(\alpha)}\int_{a}^{x}\frac{f(s)}{(x-s)^{1-\alpha}}\,ds\;$ is the Abel fractional integral operator. The subdiffusion case based on this derivative is well-documented in the literature and for background of particular relevance to inverse problems we refer to the papers [17, 15, 14].

Assuming the initial value $u_{0}$ is given, then evaluating (6) at $x_{i}\in\partial\Omega$ from the over-posed data values gives a Dirichlet series which (under specific circumstances) can lead to recovering the spectrum $\{\lambda_{n}\}$ and certain norming constants of the eigenfunction. This offers little in higher space dimensions, but in $\mathbb{R}^{1}$ this conversion to an inverse Sturm-Liouville problem can lead to a uniqueness proof and, in theory, a reconstruction algorithm. However there are serious difficulties. The inversion of the Dirichlet series to obtain $\{\lambda_{n}\}$ given the asymptotic form $\lambda_{n}\approx Cn^{2}$, is an extremely ill-conditioned problem. Thus the ability to effectively recover many eigenvalues is limited, even with extremely small values of $t$ being measured. In order to even accomplish this one must ensure that the initial data is chosen so that $\langle u_{0},\phi_{n}\rangle\not=0$ for any $n$. Since we don’t know the eigenfunctions this is difficult to guarantee other than through an argument that this is expected to occur with probability zero or using very special $u_{0}$ such as a delta distribution as in [4]. In addition, from any collection of such spectral data one can only determine a single coefficient of $-\mathbb{L}$. Indeed, the Liouville transform shows that the entire operator $-\mathbb{L}$ can be mapped into one with only a composite potential term $Q(x)$ appearing in such a manner that the original spectrum is preserved, [3].

The specification of spatial information avoids many of these drawbacks and again has been well studied. See [21, 12] for the parabolic case and the recovery of the single coefficient $q(x)$ and also [24] for the subdiffusion case. In [15] this was generalized to include a nonlinear term, namely $q(x)f(u)$ where $f(u)$ was known and the spatial factor $q(x)$ had to be determined from initial and final data. In the other direction [15] showed that a reaction term $f(u)$ could be recovered from such data. As far as the authors are aware the current work is the first attempt at showing a uniqueness theorem and a reconstruction algorithm for the case of two independent coefficients even in the parabolic case. However in the case of the elliptic operators this has been accomplished albeit in a limited setting. See, for example, [13].

The outline of this paper is as follows. In the next section we present the recovery algorithm for the pair of coefficients and then proceed by giving conditions that lead to both uniqueness and convergence. The final sections shows some quite detailed numerical experiments that show both the feasibility and the constraints of the method.

We will use a variety of spaces and norms during the analysis and make some comments on their notation here.

By $C^{k,\beta}(\Omega)$ we mean the Schauder spaces of those functions whose $k^{\mbox{{th}}}$ derivative is Hölder continuous of order $\beta$ on the set $\Omega$. The Hölder norm of a function being \;\sup_{x\in\Omega}|f(x)|+\sup_{x,y\in\Omega}\bigl{(}(|f(x)|-|f(y)|)/(|x-y|)\bigr{)}.

Sobolev spaces will be denoted by the usual notation $H^{k,p}(\Omega)$, and in case of $p=2$ simply $H^{k}(\Omega)$. However, we will have no use here for $p$ values other than $p=2$ or $p=\infty$ and the latter almost always coupled with $k=0$ giving $L^{\infty}$. Then $L^{2}=H^{0,2}$ and $L^{\infty}=H^{0,\infty}$ norms will be denoted by $\|\ldots\|_{2}$ and $\|\ldots\|_{\infty}$, and the $L^{2}$ inner product by $\langle\cdot,\cdot\rangle$. Additionally we will use the Hilbert spaces $\dot{H}^{s}(\Omega)$ defined by $\dot{H}^{s}(\Omega)=\left\{v\in L^{2}(\Omega):\sum_{j=1}^{\infty}\lambda_{j}^{s}|(v,\varphi_{j})|^{2}<\infty\right\}$ with the norm $\|v\|_{\dot{H}^{s}(\Omega)}^{2}=\sum_{j=1}^{\infty}\lambda_{j}^{s}|(v,\varphi_{j})|^{2}$, which is equivalent to the $H^{s}$ Sobolev norm for $s\in[0,2]$ for $a\in H^{1,\infty}(\Omega)$, $q\in L^{2}(\Omega)$, and $\mathbb{L}$ elliptic.

3 Algorithms and their analysis

In this section we develop the main algorithmic scheme to show uniqueness for recovery of both $a(x)$ and $q(x)$ from the overposed data under appropriate conditions. This scheme is constructive being iterative in nature with a convergence analysis possible.

3.1 An iteration scheme to recover $a$ and $q$

To recapitulate, we seek to determine both $a(x)$, $q(x)$ in

[TABLE]

with prescribed impedance boundary and given initial conditions

[TABLE]

and known forcing functions $r_{u}=r_{u}(x,t,u)$, $r_{v}=r_{v}(x,t,v)$ and reaction term $f(u)$ from observations $g_{u}(x):=u(x,T)$, $g_{v}:=v(x,T)$.

For a given $a$, $q$ we can evaluate (7) on the surface $t=T$ to obtain

[TABLE]

or

[TABLE]

where $\mathbb{M}$ is a linear operator depending only the data $g_{u}$ and $g_{v}$ and its derivatives and $\mathbb{F}$ is a nonlinear operator on $(a,q)$.

The strategy will be to provide conditions under which $\mathbb{M}$ is invertible and the combined nonlinear operator $\mathbb{T}(a,q):=\mathbb{M}^{-1}\mathbb{F}(a,q)$ is contractive.

To establish an iterative reconstruction scheme, we therefore let $\mathbb{T}(a,q)=(a^{+},q^{+})$, where

[TABLE]

and $u(x,t):=u(x,t;a,q)$, $v(x,t):=v(x,t;a,q)$ solve equation (7).

In the case $f(u)=u$ there is an obvious approach to the above; multiply the first equation in (9) by $g_{v}$, the second by $g_{u}$ and subtract thereby eliminating $q$ from the left hand side.

This gives

[TABLE]

where

[TABLE]

The value of $W$ on $\partial\Omega$ is known from the boundary conditions imposed on the system so that (12) gives an update for $a(x)$ in terms of $W$ and $\phi$. This shows that the scheme (11) can be inverted for $a$ and then by substitution, also for $q$ provided that $W$ does not vanish in any subset of $\Omega$ with nonzero measure.

If now, for example, we are in one space dimension and $u$ and $v$ share the same boundary conditions at the left endpoint, then

[TABLE]

where we used the fact that $W(0)=0$. The above will be the basis of one implementation of our reconstruction process in Section 4: namely the eliminate $q$ version.

This also works the other way around so an alternative is to multiply the first equation by $\nabla g_{v}$ and the second by $\nabla g_{u}$ and subtract giving

[TABLE]

where $\psi(x)$ has already been computed from the previous iteration.

There is a seeming symmetry between this uncoupling of $a$ and $q$ but this first impression could be misleading. In (14) we obtain an updated $a$ directly from previous iteration values of both $a$ and $q$ and this involves only the function $W$. The inversion of $a$ will go smoothly if $W$ does not vanish in $\Omega$ and even zeros of measure zero can be handled as we will see in Section 4. In (15) we obtain an updated $q$ that also depends not only on $W$ but also $\nabla W$. As we shall see, this makes the uncoupling of $q$ less stable than the other way around.

However, the important point is that the above shows that the linear operator $\mathbb{M}$ can be inverted by eliminating either of $q$ or $a$. Given this, we could also use (10) directly by inverting the linear operator $\mathbb{M}$ and solving simultaneously for $a$ and $q$ after representing these functions in a basis set. An implementation of this approach will also be shown in Section 4 and in general turns out to be the most effective approach, more clearly avoiding some of the difficulties noted above by using a least squares setting.

3.2 Contractivity in one space dimension: eigenfunction expansion

It is well known that the full Sturm-Liouville form

[TABLE]

can be placed into canonical Schrödinger form

[TABLE]

where the form of the boundary conditions is preserved and the equivalent $Q$ is given by the classical Liouville transformation, achieved by setting $Q(x)=\frac{f^{\prime\prime}}{f}(x)+L^{2}q(x)$ where $f(x)=[a]^{1/4}$, $\,L=\int^{1}_{0}[a(s)]^{-1/2}\,ds$ and $\mu=L^{2}\lambda$, see [3]. However, we do not need such regularity assumptions and the following version due to Everitt can be found in [8].

Lemma 3.1**.**

Let $a$, $r$ be such that $a$, $a^{\prime}$, $r$, $r^{\prime}$ be absolutely continuous on $(0,1)$ with $a(x)$ and $r(x)$ strictly positive. Then the mapping

[TABLE]

taking $(0,1)$ into $(0,L)$ has an inverse mapping $\ell^{-1}$ where $L=\int_{0}^{1}[r(s)/a(s)]^{1/2}\,ds$. Then the Sturm-Liouville equation (16) can be transformed into (17) by the Liouville transformation $\,y\to\ell(x)$, $\;v(y)\to v(\ell(x))=[a(x)r(x)]^{1/4}u(x)$ where

[TABLE]

Remark 3.1**.**

We will of course take $r(x)=1$ as it is not involved in our current problem and this makes the full transformation simpler. We note that the coefficients $a$ and $r$ play a complementary role and the more general case can be used if indeed our inverse problem was to determine the specific heat $r(x)$ and the potential $q(x)$ in the parabolic equation $r(x)u_{t}-u_{xx}+q(x)u$.

Remark 3.2**.**

One can extend the spaces involved even further by accepting potentials $Q(x)$ in the distribution space $H^{-1}(0,1)$ and still retain the essential features required of the eigenvalues/eigenfunctions, see [11]. This in turn allows weaker assumptions to be placed on the coefficient $a$: for example, to lie in $H^{1}(0,1)$. See, [22].

Now let $Q_{1}$ and $Q_{2}$ be in $L^{2}[0,1]$ and let $\{\hat{\phi}_{n}(x),\lambda_{n}\}_{1}^{\infty}$ and $\{\hat{\psi}_{n}(x),\mu_{n}\}_{1}^{\infty}$ denote the corresponding eigenfunction/eigenvalue pairs:

[TABLE]

where we choose the normalizations $u^{\prime}(0)=1$ if the left boundary condition is Dirichlet and $u(0)=1$ in the case of an impedance boundary condition with a finite impedance value.

Lemma 3.2**.**

For some constant $C=C(M_{Q})$

[TABLE]

Proof.

We have the asymptotic expansion [3]

[TABLE]

in the case of Dirichlet conditions where each subsequent $k^{th}$ term in the expansion is the inner product of $Q$ with $\cos 2k\pi x$ and the decay rate of each term is as indicated. The estimate (19) now easily follows

In fact this above asymptotic rate holds for $Q$ only in $L^{2}(0,1)$ but the term $t_{n}$ is replaced by $\eta_{n}\in\ell^{2}$. In the case of non-Dirichlet conditions the cosine terms are replaced by sines. ∎

Suppose that $P(x)$ and $Q(x)$ are two potentials in $L^{2}(0,1)$ with the same spectrum $\{\lambda_{n}\}$: that is, $-u^{\prime\prime}+Qu=\lambda u$ and $-v^{\prime\prime}+Pv=\lambda v$ where $u$ and $v$ share the same conditions at $x=0$. We will take these to be Dirichlet and impose the normalization that $u^{\prime}(0)=v^{\prime}(0)=1$. Then it follows from the well-known Gel’fand-Levitan formulation that

[TABLE]

where

[TABLE]

In the case of impedance boundary conditions at $x=0$, $u^{\prime}(0)-hu(0)=0$ we would normalize by $u(0s)=1$ and (3) above would be replaced by

[TABLE]

The key point here is that $K(x,t)$ does not depend on $\lambda$, only on $P$ and $Q$. It satisfies a hyperbolic partial differential equation with characteristics given by the lines $x=\pm t$. The boundary conditions in the above form a Goursat problem for the hyperbolic equation.

Suppose now that $P(x)=0$ and we have (18). For definiteness we assume Dirichlet conditions at $x=0$.

Then

[TABLE]

and the solution to (21) can be written as

[TABLE]

where $R$ is the rectangular region with corners $(x,t)$, $(\frac{x-t}{2},\frac{x-t}{2})$, $(\frac{x+t}{2},\frac{x+t}{2})$, and $(0,0)$, From the Volterra equation (22) it follows that $K$ depends on $Q$ through a Lipschitz bound $\|K(\cdot;Q_{1})-K(\cdot;Q_{2})\|_{\infty}\leq C\|Q_{1}-Q_{2}\|_{2}$. Indeed, if $Q>0$ then we have immediately from the above that $\|K(\cdot;Q_{1})-K(\cdot;Q_{2})\|_{\infty}\leq\|Q_{1}-Q_{2}\|_{2}$.

Putting the above together and noting the correspondence $\phi_{n}=\sqrt{2\lambda_{n}}\hat{\phi}_{n}$, we have the lemma

Lemma 3.3**.**

Given the above, the eigenfunctions corresponding to $Q_{1}$ and $Q_{2}$ must satisfy

[TABLE]

The point of all of this is that the given the regularity assumptions on functions $a_{1}(x)$, $a_{2}(x)$, $q_{1}(x)$, $q_{2}(x)$, the transformed functions $Q_{1}$ and $Q_{2}$ must satisfy

[TABLE]

for any $\beta>0$, and we can therefore apply the above lemmas to see that without loss of generality we may assume that both the original eigenfunctions and eigenvalues must depend on a Lipschitz manner on $a_{1}-a_{2}$ and $q_{1}-q_{2}$.

Now returning to the solution representation (6) of $u(x,t)$

[TABLE]

we see that for $t>0$ and the parabolic case

[TABLE]

This converges uniformly for $t\geq t_{0}>0$ and shows $u_{t}(x,t)$ for fixed $t$ is Lipschitz continuous in both the eigenfunctions $\{\phi_{n}\}$ and eigenvalues $\{\lambda_{n}\}$ and hence in light of the above also in the functions $q$ and $a$ in the stated norms. In addition, the exponential decay of the term $e^{-\lambda_{n}T}$ shows that for $T$ sufficiently large, this Lipschitz constant will be less than unity.

In the case of the fractional time operator we need a modification. There is no longer exponential decay of the solution as the function $E_{\alpha,1}(-\lambda T^{\alpha})$ has only linear decay in time and is a well-known difference between the classical and fractional cases especially in regards to inverse problems, [14, 17].

Now taking the time derivative of the solution means that the sequence $\{b_{n}\}$ defined in (25) must be assumed to be at least in $L^{2}(\Omega)$ and in consequence we must add the additional regularity assumption that the initial data $u_{0}\in H^{2}(\Omega)$. This follows directly from the fact that multiplying by $\lambda$ in the Fourier coefficients is equivalent to taking two derivatives in space.

We now formally show these statements below.

The difference $\mathbb{F}(a_{1},q_{1})-\mathbb{F}(a_{2},q_{2})$ can be decomposed as follows

[TABLE]

This allows an estimate, for example, of the supremum norm of the difference by

[TABLE]

for $\sigma>1/2$, where we have used the fact that

[TABLE]

with

[TABLE]

and, using $\langle\phi_{j},\phi_{n}\rangle=0$ for $j>n$,

[TABLE]

Higher norms, which are needed to recover $a$ in $C^{1,\beta}(\Omega)$ can be estimated by using continuity of the embedding $\dot{H}^{\sigma}(\Omega)\to C^{1,\beta}(\Omega)$. In spite of the arising additional powers of $\lambda_{n}$, the exponentially decaying factors will always dominate, even if the $C^{k}(\Omega)$ norm with arbitrary $k\in\mathbb{N}$ is taken. This is not the case for the Mittag-Leffler function, where due to the estimate $\frac{1}{1+\Gamma(1-\alpha)x}\leq E_{\alpha,1}(-x)$ for all $x>0$, the attainable smoothness is limited to $C^{1,\beta}$ with $\beta<1/2$.

3.3 Contractivity from the pde directly

We will reproduce the contractivity estimates from the previous section but without direct reference to the eigenfunction expansion (25) since in higher space dimensions the relationship between the coefficients and the eigenfunctions/eigenvalues is much less clear. Instead we will use the differential equations directly for the differences of both $a_{1}$, $a_{2}$ and $q_{1}$, $q_{2}$. This leads to a non-homogeneous parabolic/subdiffusion system where the right hand side depends on the coefficients themselves. Caution is needed at several places as even in the parabolic case there are difficulties. For example, in a strong solution interpolation there are regularity concerns. For a right hand side function $F\in C^{k,\beta}(\Omega)\times C^{r,\beta}(0,T)$ where $F$ contains the values of $a$ and $q$, the corresponding solutions $u_{t}$ will lie in $C^{k,\beta}(\Omega)\times C^{r,\beta/2}(0,T)$ and this regularity drop causes difficulties with the mapping properties of $\mathbb{T}$. On the other hand, weak solutions in Sobolev spaces with, say, $q\in L^{2}(\Omega)$ and $u\in H^{1}(\Omega)$ have the issue that $q\,u$ is undefined in $\mathbb{R}^{d}$ with $d>3$ and other embedding estimates needed often make further restrictions unless $d=1$. This is exactly the situation even in the case of an unknown potential as shown in [15] and clearly is more complex with a conductivity $a$ involved.

For two different coefficient pairs $(a,q)$, $(\tilde{a},\tilde{q})$ (with corresponding solutions $u$, $v$, $\tilde{u}:=u(x,t;\tilde{a},\tilde{q})$, $\tilde{v}:=v(x,t;\tilde{a},\tilde{q})$), the difference $(da^{+},dq^{+})=\mathbb{T}(a,q)-\mathbb{T}(\tilde{a},\tilde{q})$ satisfies

[TABLE]

where $\hat{u}$, $\hat{v}$ solve

[TABLE]

with $da=a-\tilde{a}$, $dq=q-\tilde{q}$.

The system (26) of PDEs for $da^{+}$ and $dq^{+}$ and its stable solution has been considered in different contexts, see, e.g., [1], see also [7] and [10], for the single coefficient case.

Multiplying the first equation in (26) by $g_{v}$, the second one with $g_{u}$, and subtracting, we get, with

[TABLE]

and using some cancellation leading to $\nabla(da^{+}W)=\nabla da^{+}\cdot W+da^{+}(g_{v}\triangle g_{u}-g_{u}\triangle g_{v})$ that

[TABLE]

While considering the simultaneous identification of $a$ and $q$, we restrict ourselves to one space dimension $\Omega=(0,L)$, since this more easily allows to resolve (29) and, in the computation of $dq^{+}$, to eliminate $\nabla da$; 333space derivatives will therefore simply denoted by a prime, whereas, to avoid additional notation, we will stay with the nabla notation when it is a partial space derivative of a space and time dependent function Moreover, we will need this restriction on the space dimension for being able to use the embedding $H^{1}(\Omega)\to L^{\infty}(\Omega)$. Later on we will consider the identification of the potential $q$ alone, which works in higher space dimensions as well, see Section 3.4. Additionally, to prove decay of $\tilde{u}$ as needed for establishing contractivity, we will focus on the homogeneous case $r_{u}=r_{v}=0$ and assume that $W$ is bounded away from zero.

Lemma 3.4**.**

Let $\Omega=(0,L)\subseteq\mathbb{R}$, $r_{u}=r_{v}=0$, $\alpha\in(0,1)$ and let $g_{u},g_{v}\in H^{2}(\Omega)$ satisfy $|\frac{1}{W(x)}|\leq C_{0}$ for $W$ defined in (28). Then any solution $(da^{+},dq^{+})$ of (26) with $da^{+}(0)=0$ satisfies the estimate

[TABLE]

Proof.

From (29)

[TABLE]

Multiplying the first equation in (26) by $g_{v}^{\prime}$, the second one with $g_{u}^{\prime}$, and subtracting, we can analogously eliminate ${da^{+}}^{\prime}$ to get, with $\tilde{W}=g_{v}g_{u}^{\prime\prime}-g_{u}g_{v}^{\prime\prime}=W^{\prime}$

[TABLE]

This yields the estimate (30). ∎

Remark 3.3**.**

From (30) it seems that the problem of recovering $q$ is more ill-posed since the estimates for $dq^{+}$ – even in a weaker norm – require higher derivatives of the data $g_{u}$, $g_{v}$. This is confirmed by the computational results.

Estimate (30) could obviously as well be obtained by assuming $da^{+}(L)=0$ instead of $da^{+}(L)=0$. To achieve that $da^{+}$ vanishes at one of the boundary points, we assume that $a(0)=:a_{0}$ or $a(L)=:a_{L}$ is known and prescribe this value as an additional condition in (11).

To estimate $\|D^{\alpha}_{t}\hat{u}(T)\|_{2}$ (and likewise $\|D^{\alpha}_{t}\hat{v}(T)\|_{2}$) appearing in the right hand side of (30), we use eigensystems $(\lambda_{j},\varphi_{j})$ $(\tilde{\lambda}_{j},\tilde{\varphi}_{j})$ of the operators defined by $\mathbb{L}w=-\nabla\cdot(a\nabla w)+q\,w$ and $\widetilde{\mathbb{L}}w=-\nabla\cdot(\tilde{a}\nabla w)+\tilde{q}w$, to obtain the representations

[TABLE]

where we have assumed $r_{u}=0$ in order to obtain a proper decay of $D^{\alpha}_{t}\tilde{u}(t)$.

We first of all provide an estimate of the right hand side of the equation (27) for $\hat{u}$, i.e., of the inhomogeneity in (31).

For this purpose we will make use of the Poincaré-Friedrichs type inequality

[TABLE]

and the assumptions

[TABLE]

for some constant $\underline{q}$ (not necessarily positive) with

[TABLE]

as well as

[TABLE]

Lemma 3.5**.**

The function $\tilde{u}$ defined by (32), with (34), (35), (36) and $\widetilde{\mathbb{L}}u_{0}\in L^{2}(\Omega)$ satisfies the estimate

[TABLE]

for some $C_{2}$ depending only on $\Omega$ and $\frac{1}{\underline{a}}$, where $\;{\displaystyle e(\tilde{\lambda}_{1},t):=\sup_{\mu\geq\tilde{\lambda}_{1}}\max\{1,\mu\}E_{\alpha,1}(-\mu t^{\alpha})\,.}$

Proof.

We use the fact that $\tilde{u}$ solves the first of the two equations (7) with $a=\tilde{a}$ and $r_{u}=0$, as well as homogeneous impedance boundary conditions and start by estimating

[TABLE]

for $C_{1}=\max\{1,C_{H^{1}\to L^{\infty}}^{\Omega}\|(-\triangle+\mbox{id})^{-1}\|_{L^{2}(\Omega)\to H^{2}(\Omega)}\}$. For the $\nabla D^{\alpha}_{t}\tilde{u}(t)$ term in (38), from the fact that $\tilde{z}:=D^{\alpha}_{t}\tilde{u}$ solves

[TABLE]

with impedance boundary conditions and using integration by parts, we get

[TABLE]

With (33), (34), (35) we obtain

[TABLE]

Moreover, for the $\triangle D^{\alpha}_{t}\tilde{u}(t)$ term in (38), we use the fact that for any $w\in H^{2}(\Omega)$

[TABLE]

hence by our assumption (36) on the smallness of $\tilde{a}^{\prime}$,

[TABLE]

which we apply with $w=D^{\alpha}_{t}\tilde{u}(t)$. In here, by (39),

[TABLE]

where

[TABLE]

so that all spatial derivatives of $D^{\alpha}_{t}\tilde{u}$ that are needed for further estimating (38) can be expressed via $D^{\alpha}_{t}\tilde{u}$ and $(D^{\alpha}_{t})^{2}\tilde{u}$.

For these, we get from

[TABLE]

that

[TABLE]

Using this and (40), (41), (42), (43) in (38), we get (37). ∎

Lemma 3.6**.**

The function $\hat{u}$ defined by (31), with $\tilde{u}$ according to Lemma 3.5 satisfies the estimate

[TABLE]

with a constant $C_{3}$ depending only on $\Omega$, $\frac{1}{\underline{a}}$ and $\|\widetilde{\mathbb{L}}u_{0}\|_{2}$ and

[TABLE]

Proof.

Using the fact that $\hat{z}:=D^{\alpha}_{t}\hat{u}$ solves

[TABLE]

with initial conditions $\hat{z}(0)=D^{\alpha}_{t}\hat{u}(0)=\nabla(da\nabla u_{0})-dq\,u_{0}$, we get

[TABLE]

Take the $L^{2}$ norm of the above expression for $D^{\alpha}_{t}\hat{u}$ to get, using an inequality of the form

[TABLE]

that

[TABLE]

Thus from (37) we get (44) with (45). ∎

It is therefore crucial for contractivity to prove that $\Phi(T)$ as defined in (45) tends to zero for increasing $T$.

Lemma 3.7**.**

For $\Phi(T)$ according to (45) we have $\Phi(T)\to 0$ as $T\to\infty$.

Proof.

Using the identities

[TABLE]

and the bound

[TABLE]

that hold for every $\alpha\in(0,1)$ and all $s\in\mathbb{R}$ and all $x\in\mathbb{R}^{+}$, as well as the complete monotonicity of the function $x\mapsto E_{\alpha,1}(-x)$ on $\mathbb{R}^{+}$, see, e.g., [5, 14], we can estimate as follows

[TABLE]

where we have used the substitution $s=Tr$. Hence we get

[TABLE]

∎

This together with (30), (44) yields contractivity of $\mathbb{T}$ for sufficiently large $T$ on the closed convex set

[TABLE]

with the norm

[TABLE]

In the definition of $B$, we assume that $\underline{a}>0$, $\underline{q}\in\mathbb{R}$, $\rho>0$ satisfy (35) and

[TABLE]

Boundedness away from zero of $\lambda_{1}=\lambda_{1}(a,q)$, uniformly for $(a,q)\in B$ follows from the identity

[TABLE]

hence

[TABLE]

which together with the constraints (35) and (50) on $\underline{a}$, $\underline{q}$ and $\varrho$ yields

[TABLE]

Thus (44) remains valid with $\min\{\lambda_{1},\tilde{\lambda}_{1}\}$ replaced by $\underline{\lambda}$ and we get contractivity with a uniform constant on $B$, provided $T$ is sufficiently large.

Theorem 3.1**.**

*Let $\Omega=(0,L)\subseteq\mathbb{R}$, $r_{u}=r_{v}=0$, $\alpha\in(0,1)$ and let $g_{u},g_{v}\in H^{2}(\Omega)$ satisfy $|\frac{1}{W(x)}|\leq C_{0}$ for $W$ defined in (28).

Then for $T>0$ sufficiently large, the operator $\mathbb{T}$ defined by $\mathbb{T}(a,q)=(a^{+},q^{+})$ solving (11) with $a^{+}(0)=a_{0}$ or $a^{+}(L)=a_{L}$, is a contraction on the set (48) with respect to the norm defined by (49).*

Remark 3.4**.**

Under the conditions of Theorem 3.1 we have convergence of the fixed point iteration defined by $(a_{k+1},q_{k+1})=\mathbb{T}(a_{k},q_{k})$ for any starting value in $B$. Moreover, we have uniqueness of a solution to the inverse problem on $B$.

Remark 3.5**.**

In the realistic setting of noisy data $g_{u}^{\delta}\approx g_{u}$, $g_{v}^{\delta}\approx g_{v}$, $g_{u}^{\delta},g_{v}^{\delta}\in L^{2}(\Omega)$ with $\|g_{u}-u(\cdot,T;a_{act},q_{act})\|_{2}\leq\delta$, $\|g_{v}-v(\cdot,T;a_{act},q_{act})\|_{2}\leq\delta$, we first of all filter the data to get approximations $\tilde{g}_{u}^{\delta},\tilde{g}_{v}^{\delta}\in H^{2}(\Omega)$ with

[TABLE]

where $\tilde{\delta}=O(\psi(\delta))$ with an index function $\psi$444i.e., a continuous non-decreasing function $\psi:(0,\infty)\to\mathbb{R}^{+}$ with $\lim_{x\to 0}\psi(x)=0$ depending on the smoothness of both $v(\cdot,T;a_{act},q_{act})$, $u(\cdot,T;a_{act},q_{act})$. Then we perform the fixed point iteration up to a stopping index $k_{*}(\tilde{\delta})\sim\log(1/\tilde{\delta})$, see, e.g., [16, Section 3.5] to obtain a convergence rate $|||(a_{k_{*}(\tilde{\delta})},q_{k_{*}(\tilde{\delta})})-(a_{act},q_{act})|||=O(\tilde{\delta})$.

3.4 Contractivity in the potential only case

Contractivity, hence convergence of the fixed point iteration and uniqueness can be proven also in higher space dimensions $\Omega\subseteq\mathbb{R}^{3}$ in case the diffusion coefficient $a$ is known, see also [15]. This is due to the fact that the regularity requirements on $D^{\alpha}_{t}u(T)$ can be weakened to just an $L^{\infty}$ estimate as we only recover an $L^{2}$ coefficient $q$, while we needed to estimating $D^{\alpha}_{t}\hat{u}(T)$ in $H^{2}(\Omega)$ for obtaining also $a$ in $H^{1}(\Omega)\cap W^{1,\infty}(\Omega)$. Since the proof in fact very much follows the line of the previous section, we will here only show the key steps.

To this end, we consider the problem of identifying $q(x)$ in

[TABLE]

with homogeneous impedance boundary conditions

[TABLE]

from observations $g(x)=u(x,T)$. Note that $-\triangle$ can be replaced by an arbitrary second order elliptic differential operator with sufficiently smooth coefficients.

Define a fixed point operator $\mathbb{T}$ by

[TABLE]

For two different potentials $q$, $\tilde{q}$ (with corresponding solutions $u$, $\tilde{u}:=u(x,t;\tilde{q})$, the difference $dq^{+}=\mathbb{T}(q)-\mathbb{T}(\tilde{q})=\frac{-D^{\alpha}_{t}u(T)+D^{\alpha}_{t}\tilde{u}(T)}{g}:=-\frac{D^{\alpha}_{t}\hat{u}(T)}{g}$ where $\hat{u}$, solves

[TABLE]

with $dq=q-\tilde{q}$.

We use eigensystems $(\lambda_{j},\varphi_{j})$ $(\tilde{\lambda}_{j},\tilde{\varphi}_{j})$ of the operators defined by $\mathbb{L}w=-\triangle w+q\,w$ and $\widetilde{\mathbb{L}}w=-\triangle w+\tilde{q}w$, to obtain the representations

[TABLE]

Here

[TABLE]

for $C_{0}=C_{\dot{H}^{\sigma}\to L^{\infty}}^{\Omega}\|\widetilde{\mathbb{L}}u_{0}\|_{2}$. Hence

[TABLE]

which for $\sigma\leq 2$, by Lemma 3.7, tends to zero as $T\to\infty$. Due to the fact that $H^{\sigma}(\Omega)$ continuously embeds into $L^{\infty}(\Omega)$ for a bounded domain $\Omega\subseteq\mathbb{R}^{3}$, this is works in up to three space dimensions and gives contractivity for $T$ large enough.

Theorem 3.2**.**

Let $\Omega\subseteq\mathbb{R}^{3}$ be a bounded $C^{1,1}$ domain, $r=0$, $\alpha\in(0,1)$ and let $g\in H^{2}(\Omega)$ satisfy $|\frac{1}{g(x)}|\leq C_{0}$. Then for $T>0$ sufficiently large, the operator $\mathbb{T}$ defined by (52) is a contraction with respect to the $L^{2}$ norm.

4 Reconstructions

We will show the results of numerical experiments with the three versions of the basic iterative scheme: compute $a$, $q$ in parallel; eliminate $q$ recover $a$; eliminate $a$ recover $q$,

In the reconstructions to be shown we used the initial values $\alpha=1$, $T=0.5$ and a noise level (uniformly distributed) of 1% as a basis for discussion. At the end of the section we will indicate the effective dependence of the reconstruction process on these quantities. The reconstructions we show will be in $\mathbb{R}$ as the graphical illustration is then more transparent and there is little to be gained technically or visually from higher dimensions.

We will also take the following actual functions to be reconstructed as

[TABLE]

As data we took two differing initial values $u_{0}(x)$ and $v_{0}(x)$ and as boundary conditions we used (nonhomogeneous) Dirichlet at the left endpoint and Neumann at the right; typically different for each of $u$ and $v$.

One such data set is shown in Figure 1

4.1 Performance of the three schemes

In the parallel scheme with actual values, since we make no constraints on the form of the unknown functions other than sufficient regularity, we do not choose a basis with in-built restrictions as would be obtained from an eigenfunction expansion. Instead we used a radial basis of shifted Gaussian functions $b_{j}(x):=e^{-(x-x_{j})^{2}/\sigma}$ centered at nodal points $\{x_{j}\}$ and with width specified by the parameter $\sigma$. The linear operator $\mathbb{M}$, in equation (10) then takes the matrix form

[TABLE]

where $A_{1}$ denotes the representation of $a(x)$ using the values $g_{u}$ and $Q_{2}$ the representation of $q(x)$ using the values $g_{v}$.

The sequential schemes are based on eliminating one of $a(x)$ or $q(x)$ and having $\mathbb{M}$ represented through pointwise values of the functions $W$ and $W^{\prime}$.

The singular values of the component matrices $A_{1}$ and $Q_{1}$ are shown in the leftmost figure in Figure 2 and the functions $W$ and $W^{\prime}$ in the rightmost figure .

Even before seeing the resulting reconstructions it is clear the far superior conditioning of the $A$ matrix over that of $Q$ – a factor of over 10 in the larger singular values and of 100 in the lower ones – is going to significantly favour the reconstruction of the $a(x)$ coefficient. This is also borne out from the rightmost figure here: while the values of the function $W$ are modest, there is a much larger range in the values of $W^{\prime}$. Note that from equation (14) the reconstruction of $a(x)$ requires only $W$, but that of (15) which updates $q(x)$ requires both $W$ and $W^{\prime}$. Indeed this turns out to be the situation in both cases.

Using the parallel scheme reconstructions of $a$ and $q$ from two initial/boundary values under $1\%$ random uniform noise are shown in Figure 3 The initial approximations were $a(x)=1$ and $q(x)=0$. The first iteration resulted in an already near perfect reconstruction of the $a(x)$ coefficient but that of $q(x)$ lagged significantly behind and in the end the error in the data measurements were predominately in this coefficient.

In the eliminate $q$ and update $a$ strategy some care must be taken here as the function $W(x)$ has two zeros; an interior one around $\tilde{x}=0.17$ and at the endpoint $x=1$. Thus a straightforward division by $W$ to recover $a$ from (14) isn’t possible. In theory the right-hand side $\Phi(x)=\int_{0}^{x}\phi(s)\,ds$ term should also vanish at these points including the interior one $\tilde{x}$, but with data noise this isn’t going to be the case. In practice we found the following device worked very well. Remove a small segment of several grid points from an interval $I_{\tilde{x}}=(\tilde{x}-\delta,\tilde{x}+\delta)$, set $\Phi(\tilde{x})=0$ and then fill in the interval by interpolation using say a smoothing spline with the level of smoothing chosen depending on the assumed noise in $\phi$. Then a straight division recovers $a(x)$. A new direct solve then is used to recover the next iteration of $q(x)$. As Figure 4 shows the results are comparable to the previous reconstruction.

It is worth observing from Figures 1 and 4, that those regions where the reconstructions of $q$ are poorest coincide with regions of smaller values of $g_{u}(x)=u(x,T)$ and $g_{v}(x)=v(x,T)$, namely near the left hand point of the interval. This is in keeping with the fact that both $W$ and $W^{\prime}$ are smaller in magnitude at these points.

We do not show a reconstruction for the version based on the update equation (15). While a somewhat satisfactory reconstruction of the coefficient $a(x)$ was obtained, the scheme failed to converge for the coefficient $q(x)$. This fact alone shows how dominant a role the diffusion coefficient plays in the process at the expense of the much weaker potential term. It is only under significantly less data noise that an effective reconstruction of the latter was possible. The relative rates of convergence of all three versions is shown in the displayed table.

The difference in the iteration counts in the table is due to the use of a discrepancy principle as a stopping rule, which terminates the iteration as soon as the residual drops below the noise level. The discrepancy principle is well established as a regularization parameter choice for ill-posed problems see, e.g., [6], [18]. Here note that regularization is done by smoothing the data here (cf. Remark 3.5) so the stopping index does not act as a regularization parameter; merely the discrepancy principle appears to find a good final iterate.

4.2 Changing some of the parameters

The above reconstructions were set with the final time taken to be $T=0.5$; the question is how the schemes would progress for different $T$ In particular, we are interested in this as one should expect the contraction constant (if indeed there is a contraction) to be smaller with increasing $T$.

The answer is much as expected; the contraction constant varies with $T$, and more generally, so does the strength of the nonlinear contribution terms $u_{t}(x,T;a,q)$ and $v_{t}(x,T;a,q)$. We would thus expect more iterations to be required as $T$ was reduced and indeed this is the case. As an illustration of this effect, Figure 5 shows the first iterate of $a_{1}(x)$, for the values of $T=0.1$ and $T=0.05$. Recall from Figure 3 that even the first iteration $a_{1}$ was sufficiently close to the actual when $T=0.5$ so that it was barely indistinguishable from the actual. We do not show the reconstructions of $q(x)$ here as in neither case of $T=0.1$ nor of $T=0.05$ did these converge. Indeed for these smaller value of $T$ the iterations oscillated widely without any sense of convergence, it being quite clear that we were in a region of non-convergence.

Since the coefficient $a(x)$ is not the determining factor in the reconstruction process, as a function of $\alpha$ the scheme reverts to our ability to determine $q$ exactly as in [15]. Thus, as expected, the contraction constant increases with decreasing $\alpha$; quite sharply at first. Hence the convergence rate decreases with decreasing $\alpha$.

The variation of the schemes with noise level is now predictable. With an error much in excess of $1\%$ the schemes degrade rapidly. Even for errors less than this value, the schemes based on pointwise evaluation and requiring $W^{\prime}$ did quite poorly, especially for smaller values of the time measurement $T$ and the non-parabolic case.

4.3 Measurements at two later times

Finally, we indicate a little about the possibility of a single data run but taking spatial measurements at two later times $t=T_{1}$ and $t=T_{2}$. Some reflection will show the inherent difficulties here. Unless there is some dynamic change in the solution profile for $T_{1}<t<T_{2}$ the likelihood is that $g_{1}(x)=u(x,T_{1})$ and $g_{2}(x)=u(x,T_{2})$ will not differ sufficiently to allow effective division by the “Wronskian” function $W=g_{1}g_{2}^{\prime}-g_{1}^{\prime}g_{2}$. Such a change could occur by an input of a large nonhomogeneous forcing function or a change in the boundary conditions between these time values.

Consider the case when there are no such changes; we use the first solution $u(x,t)$ above but take $g_{1}=u(x,T/2)$ and $g_{2}=u(x,T)$, with $T$ again the value $T=0.5$. The story here is told by looking at a plot of the resulting functions $W$ and $W^{\prime}$ as shown in Figure 6 which highlights the weak linear independence of the two data values.

This shows $W$ is very small in magnitude with several zeros; $W^{\prime}$ is in consequence highly oscillatory and relatively large in magnitude. These are exactly the features one would try to avoid in making up a data set.

4.4 Inclusion of a known nonlinearity

In [15] it was shown that it was possible to reconstruct the potential $q(x)$ occurring in the term $q(x)f(u)$ where $f(u)$ is assumed known. Thus equation (2) seeks to generalize this to also include the determination of the conductivity $a(x)$.

Some of the methods described in section (3.1) will no longer quite work as, for example, the possibility to eliminate $q$ isn’t directly feasible as in the case of $f(u)=u$. However the scheme based on basis functions and determining $a$ and $q$ in parallel goes through exactly as before since the operator $\mathbb{M}$ remains linear due to the fact that we are evaluating the known $f$ at the data values $g_{u}$ and $g_{v}$.

Reconstructions for this case do depend quite strongly on $f$. If $f(u)$ decays more rapidly than $u$ itself then it is quite likely the recovery of $q(x)$ will be even more challenging. On the other hand if $f(u)$ is substantially larger than $u$ (but assuming the resulting solutions of equation (2) remain bounded) then the term involving $q(x)$ plays a more substantial role and the recovery of $q$ is much improved.

Acknowledgment

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

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