Augmented data for the Backus problem

TL;DR

This paper investigates an augmented boundary data approach for the Backus problem involving the Laplace equation, using additional boundary information to improve inverse source estimation and numerical recovery of harmonic functions.

Contribution

It introduces an augmented data framework for the Backus problem, deriving a quasi-linear boundary equation and applying finite element methods for numerical solutions.

Findings

Additional boundary data helps estimate sources more accurately.

Derived a quasi-linear boundary equation involving augmented data.

Numerical methods successfully recover harmonic functions from augmented data.

Abstract

Backus considered a boundary value problem for the Laplace equation with the non-linear data in the form of the magnitude $|Du|$ of the gradient of the solution u. We consider this problem with the data expanded by $(\partial/\partial\nu)|Du|$ given on the boundary of the domain. To justify the requirement for additional data, we use them to estimate the number of sources for the related inverse source problem in the plane. We show that, for an arbitrary dimension, a harmonic function satisfies a quasi-linear equation on the boundary of the domain with the coefficients involving the augmented data. We use the finite element method to recover the harmonic function on the boundary by solving numerically the derived equation.

Tables

Table 1. Table 1 . Estimated errors and rates of convergence

	$u_0$		$u_1$		$u_0+u_1$
$H_{max}$	$L^2$	$H^1$	$L^2$	$H^1$	$L^2$	$H^1$
0.2	1.58e-03	2.04e-02	3.00e-04	2.08e-03	1.56e-03	2.17e-02
0.1	2.84e-04	7.90e-03	8.59e-05	1.01e-03	2.86e-04	8.42e-03
0.05	6.40e-05	3.50e-03	2.13e-05	4.05e-04	6.63e-05	3.73e-03
0.025	1.12e-05	1.23e-03	5.31e-06	1.53e-04	1.29e-05	1.32e-03
Rate	2.3562	1.3302	1.9470	1.2602	2.2858	1.3290

Equations

$$p$$

$$p_1$$

$$p_2$$

$$p_1{u}_1-p_2{u}_2-2p{u}_{11}$$

$$p_2{u}_1+p_1{u}_2-2p{u}_{12}$$

$$\Delta u=f,$$

$$f=\sum _{|\alpha |\le N}\sum _{j=1}^{N_{\alpha }}{b}_j^{\alpha }{D}^{\alpha }{\delta }_{x_j^{\alpha }},$$

$$|\alpha |={\alpha }_1+\cdots +{\alpha }_d\text{and}{D}^{\alpha }=\frac{{\partial }^{{\alpha }_1}}{\partial x_1^{{\alpha }_1}}\dots \frac{{\partial }^{{\alpha }_d}}{\partial x_d^{{\alpha }_d}},$$

$$u(x)=u_0(x)+\sum _{|\alpha |\le N}\sum _{j=1}^{N_{\alpha }}{b}_j^{\alpha }{D}^{\alpha }\Gamma (x-x_j^{\alpha }),$$

$$-\frac{1}{2\pi }{\int }_{\partial \Omega }\frac{q}{2p}d\tau \le \sum _{|\alpha |\le N}(|\alpha |+1)N_{\alpha },$$

$$\frac{1}{2\pi i}{\int }_{\gamma }\frac{f^{\prime }}{f}dz=n_{+}-n_{-},$$

$$u_x$$

$$u_y$$

$$u(z)=u(\rho \cos \theta ,\rho \sin \theta )={\int }_0^{\rho }{u}_x\cos \theta +u_y\sin \theta dr.$$

$$u(z)=\frac{{\rho }^n}{n}(c_1\cos n\theta +c_2\sin n\theta )+o({\rho }^n),$$

$$u(z)=u(z_1)-{\int }_{\rho }^r{u}_x\cos \theta +u_y\sin \theta dr,$$

$$\frac{1}{2\pi i}{\int }_{\partial \Omega }\frac{f^{\prime }}{f}dz$$

$$=\frac{1}{2\pi i}{\int }_{\partial \Omega }\frac{1}{2}\left(\frac{\partial }{\partial x}\log p-i\frac{\partial }{\partial y}\log p\right)\cdot (dx+idy)$$

$$=\frac{1}{4\pi }{\int }_{\partial \Omega }\frac{\partial }{\partial x}\log pdy-\frac{\partial }{\partial y}\log pdx$$

$$=\frac{1}{4\pi }{\int }_{\partial \Omega }\frac{\partial }{\partial \nu }\log pd\tau$$

$$=\frac{1}{2\pi }{\int }_{\partial \Omega }\frac{q}{2p}d\tau$$

$$-\frac{1}{2\pi }{\int }_{\partial \Omega }\frac{q}{2p}d\tau =n_{-}-n_{+}\le n_{-}=\sum _{|\alpha |\le N}(|\alpha |+1)N_{\alpha },$$

$${\Delta }_{k}u=\sum _{i=1}^k{D}_{ii}u\text{and}\Delta u={\Delta }_{n}u,$$

$${\Delta }_{n+1}u=\Delta u+D_{n+1,n+1}u=0.$$

$$p=|Du{|}^2+(D_{n+1}u{)}^2\text{and}q=D_{n+1}p.$$

$$\mathrm{div}\frac{Du}{\sqrt{p-|Du{|}^2}}+\frac{1}{2}\frac{\sigma q}{p-|Du{|}^2}=0$$

$$2D_{n+1}u{D}_{i,n+1}u=D_{i}p-2D_{j}u{D}_{ij}u,i=1,\dots ,n,$$

$$2D_{i}u{D}_{i,n+1}u+2D_{n+1}u{D}_{n+1,n+1}u=D_{n+1}p=q,$$

$$D_{i}u(D_{i}p-2D_{j}u{D}_{ij}u)+2(p-(D_{i}u{)}^2)(-\Delta u)=q{D}_{n+1}u.$$

$$(p-|Du{|}^2)\Delta u+D_{i}u{D}_{j}u{D}_{ij}u-\frac{1}{2}{D}_{i}p{D}_{i}u+\frac{q}{2}\sigma \sqrt{p-|Du{|}^2}=0,$$

$$(1+|Du{|}^2)\Delta u-D_{i}u{D}_{j}u{D}_{ij}u=0$$

$$\Delta u-\frac{D_{i}u{D}_{j}u}{1-\frac{\gamma -1}{2}|Du{|}^2}{D}_{ij}u=0,$$

$$\Gamma (x,y,z)=-\frac{1}{4\pi }\frac{1}{\sqrt{x^2+y^2+z^2}}$$

$${\int }_{\Omega }\frac{Du\cdot Dv}{\sqrt{p-|Du{|}^2}}+\frac{1}{2}\frac{\sigma qv}{p-|Du{|}^2}dx=0\text{for }v\in C_0^1(\Omega ).$$

Full text

1. Introduction

Backus [2] considered the problem of determining a harmonic function $u$ given the magnitude $|Du|$ of the gradient. This problem has geophysical applications in gravimetric and geomagnetic (intensity) surveys. For example in magnetometry, in the absence of currents and away from singularities, Maxwell’s equation for the magnetic field reduces to the Laplace equation. The intensity of a magnetic or gravitational field is more readily available than the full field data, which also includes the field direction. Other potential applications can be found in medicine, specifically in magnetoencephalography and magnetocardiography [7]. Traditionally in these fields, the full magnetic field is measured using superconducting quantum interface devices (SQUID). These devices are bulky and their use is associated with a high cost of acquisition and maintenance while scalar magnetometers, which measure only the magnitude of the magnetic field, are in comparison compact, scalable, and capable of achieving a high level of sensitivity.

The focus of Backus’ paper is on the data in the form of $|Du|$ given on the boundary of the domain. Backus considered the problem for the cases of a bounded domain, an unbounded domain with the bounded exterior, and the half-plane. For the exterior of a sphere in $\mathbb{R}^{3}$, Backus proved [3] that the problem admits multiple solutions and additional data are required for well-posedness. In two dimensions, the problem for bounded domains has a unique solution [2] provided we are given the locations of all zeros of $Du$ in $\Omega$ and the direction of $Du$ at one point in $\Omega$. This additional requirement results in a well-posed problem, however the data of this type may be viewed as exogenous. Another example of expanded data is the requirement of a definite sign for the normal derivative of the solution on the boundary of an unbounded domain. In this case, Backus established uniqueness under the constraint that the solution vanishes at infinity [2], which itself can be interpreted as an additional boundary condition. Cherkaeva [4] considered a linearized problem for the exterior of a ball in $\mathbb{R}^{3}$. She obtained uniqueness by requiring additional data in the form of the potential or its normal derivative on the magnetic equator of the sphere. The magnetic equator in this context is defined as the curve on the sphere on which the vertical component of the main field vanishes. A related result was obtained by Jorge and Magnanini [8]. They proved that, if $u$ and $v$ are two harmonic functions defined outside the unit ball and regular at infinity such that $|Du|=|Dv|$ on $S^{2}$ and $u=v$ on the set $\{x\in S^{2}:\ (\partial/\partial\nu)(u+v)=0\}$, then $u=v$. Also for the exterior problem, Backus proved that, if $|Du|=|Dv|$ throughout the domain, then $u=v$. The last result sets a precedent for assuming the knowledge of $|Du|$ and hence its derivatives on a larger set in a special formulation of the problem.

We propose an intermediate assumption, namely, we require the data in the form of $|Du|$ and its normal derivative of $(\partial/\partial\nu)|Du|$ to be given on the boundary of the domain. The data in the equivalent form $(p,Dp):=(|Du|,D|Du|)$ was considered by this author in [6] for planar domains. In the case of two independent variables, the non-linear equations relating the data $(p,Dp)$ and the solution $(Du,D^{2}u)$ take the form

[TABLE]

where the subscripts denote partial derivatives, i.e., $u_{j}=\partial u/\partial x_{j}$, $u_{jk}=\partial^{2}u/\partial x_{j}\partial x_{k}$, etc. These equations are equivalent to a system of linear equations for the components of the solution $(Du,D^{2}u)$ with the components of the data $(p,Dp)$ appearing as coefficients:

[TABLE]

This system can be reduced to a system of two ordinary differential equations and further rewritten as a single ordinary differential equation in one complex variable, which simplifies the analysis and lends itself readily to reconstruction methods including analytic solution formulas and numerical computation.

The goal of the present paper is to extend the recovery of $(u,Du)$ from $(p,Dp)$ to higher dimensions. The organization is as follows. In Section 2, we consider a relevant inverse source problem and revisit an estimate from [6] on the number of sources in terms of the expanded data in the two-dimensional case. In the case of a general dimension $d\geq 2$, we derive in Section 3 a non-linear elliptic equation satisfied by a harmonic function $u$ on a hyperplane with coefficients involving $|Du|$ and $(\partial/\partial\nu)|Du|$. We report in Section 4 on numerical studies of recovering $u$ (and $Du$) on $\partial\Omega$ from the knowledge of $|Du|$ and $(\partial/\partial\nu)|Du|$ (or $D|Du|$). The method is based on solving the non-linear equation derived in Section 3.

In addition to the already mentioned references [2, 3, 8, 6], the Backus problem was studied by Díaz, Díaz, and Otero [5]. These authors considered an oblique derivative problem related to the exterior Backus problem by the Kelvin transform. In particular, they proved the maximality of the solution of the exterior problem with a definite sign of the normal derivative. Lieberman [11] addressed the regularity of solutions for a class of nonlinear boundary value problems that includes the Backus problem on bounded domains. Payne and Schaeffer [14] obtained bounds for solutions of parabolic equations and their gradient with the magnitude of the gradient prescribed on the boundary of bounded domains. Kaiser and Neudert [10] and Kaiser [9] considered the exterior problem with the direction of the gradient given on the boundary and the magnitude unknown. The data in the form of the magnitude of the gradient was also studied in connection with the problem of electrical impedance tomography by Nachman, Tamasan, and Timonov [12, 13]. These authors considered the problem of finding the conductivity $\sigma$ from the measurement of the magnitude of the current $|J|$ where the current is given by $J=\sigma Du$ and $u$ is a solution of the elliptic equation $D\cdot(\sigma Du)=0$. In addition, either the Cauchy data [12] or the Dirichlet data [13] are assumed to be given on the boundary.

We point out that the collection of data in the form of the gradient of $|Du|$ for a potential $u$ is already implemented in magnetic surveys in archaeological geophysics and in gravimetric surveys with the help of devices called gradiometers. These devices consist of two or more sensors forming a basis for a finite difference approximation of the gradient of the field magnitude.

2. Estimate for the number of poles

The recovery of $(u,Du)$ from $(p,Dp)$ on the boundary may be viewed as a first step in the identification of points sources located inside the domain given the magnitude of the field data on the boundary. The second step of finding the potential assuming a specific form of the solution can be classified as an inverse source problem for the Poisson equation with point sources.

To state this problem more precisely let $\Omega$ be a domain in $\mathbb{R}^{d}$ where $d\geq 2$ and consider the equation

[TABLE]

with the source term $f$ of the form

[TABLE]

where $\alpha=(\alpha_{1},\dots,\alpha_{d})\in(\mathbb{N}\cup\{0\})^{d}$ is a multiindex with

[TABLE]

$N$, $N_{\alpha}\in\mathbb{N}$; $b^{\alpha}_{j}\in\mathbb{R}$, $x^{\alpha}_{j}\in\Omega$, for $j=1,\dots,N^{\alpha}$; and $\delta_{x}$ is the Dirac delta function supported at $x$. The problem is as follows: given the potential and its normal derivative on the boundary of the domain, find the locations $x_{j}^{\alpha}$ and moments $b_{j}^{\alpha}$ of the point sources. The problem has been extensively studied and many results are available.

The number of point sources in the inverse source problem is often unknown. In the two-dimensional case studied in [6], the necessary and sufficient conditions for the existence of solution to the ordinary differential equation arising from (2) are related to the estimation of the number of dipoles and monopoles counting their multiplicities. Below, we extend the estimate to the case of poles of arbitrary degrees for dimension two and make it more precise.

Theorem 1**.**

Suppose $\Omega$ is a smooth domain in $\mathbb{R}^{2}$ and $u$ is a solution of (3)-(4) in $\Omega$, that is,

[TABLE]

where $u_{0}$ is a harmonic function on $\bar{\Omega}$, $\Gamma$ is the fundamental solution of the Laplace equation, and $x_{j}^{\alpha}\in\Omega$. Let $p=|Du|^{2}$, $q=\partial p/\partial\nu$ and assume that $p>0$ on $\partial\Omega$. Then,

[TABLE]

Moreover, the equality holds if $p$ does not vanish in $\Omega$.

Proof.

To obtain estimate (6) we apply the argument principle of complex analysis  [1, p. 152]. Namely, if $f$ is a meromorphic function and $\gamma$ is a simple closed curve that does not pass through zeros or poles of $f$, then

[TABLE]

where $n_{+}$ and $n_{-}$ are the numbers of zeros and poles of $f$, respectively, located inside the region with the boundary $\gamma$.

We apply the argument principle to the function $f=u_{y}+iu_{x}$. First, we discuss how the poles and zeros of $u$ correspond to those of $f$. Since $u$ is harmonic away from its singularities, $f$ is meromorphic and its poles coincide with the poles of $u$. Moreover, if $u$ has a zero or pole at $z_{0}=(x_{0},y_{0})$, then, in polar coordinates $(\rho,\theta)$ at $z_{0}$ and with $z=x+iy=\rho e^{i\theta}$, we claim that

• •

for $n\neq 0$, $u=\mathcal{O}(\rho^{n})$ and $u=o(\rho^{n-1})$ as $\rho\to 0$ if and only if $f(z)\sim z^{n-1}$ as $z\to z_{0}$;

• •

$u=\mathcal{O}(\log\rho)$ and $u=o(\rho^{-1})$ as $\rho\to 0$ if and only if $f\sim z^{-1}$ as $z\to z_{0}$.

To prove the claim in one direction, let $v$ be such harmonic conjugate of $u$ that, for $F=iu-v$, we have $F(z)\sim z^{n}$ for $n\neq 0$ or $F(z)\sim\log z$ as $z\to z_{0}$. Then $f=F^{\prime}$ and, therefore, $f(z)\sim z^{n-1}$. For the other direction, suppose $f(z)=z^{n-1}(c_{2}+ic_{1})+o(z^{n-1})$ for some constants $c_{1}$ and $c_{2}$, not both zero. Then

[TABLE]

We consider the cases of zeros and poles separately. In the case of a zero at $z_{0}$, i.e., for $n\geq 1$, integrating along the line segment connecting $z_{0}$ and $z$, we have

[TABLE]

The substitution of (7) into the right-hand side of this representation yields

[TABLE]

i.e., $u=\mathcal{O}(\rho^{n})$ and $u=o(\rho^{n-1})$ as $\rho\to\infty$. In the case of a pole at $z_{0}$, i.e., for $n\leq 0$, suppose $f$ has no other poles in the ball of radius $r>0$ centered at $z_{0}$ so that $\sup_{|z-z_{0}|=r}u(z)<\infty$. Integrating along the line segment connecting $z$ and $z_{1}=z_{0}+r(z-z_{0})/|z-z_{0}|$, we have

[TABLE]

which, after the substitution of (7) into the right-hand side, becomes (8) when $n\leq-1$ and $u(z)=c_{1}\log\rho+o(\log\rho)$ when $n=0$ and the claim is proved. It follows in particular from this claim that $f$ has a pole of order $n+1$, for some $n\geq 0$ if and only if $u\sim D^{\alpha}\Gamma$ with $|\alpha|=n$.

Performing the computation for the argument principle, we have

[TABLE]

Hence

[TABLE]

where the last identity follows from the claim above.

∎

3. A quasi-linear equation

In this section, we show that harmonic functions satisfy a quasi-linear elliptic equation on a flat portion of the boundary of a domain. With this result, we reduce the linear Laplace equation to a non-linear equation but in a lower dimensional space.

To distinguish the Laplace operators in $\mathbb{R}^{n}$ and $\mathbb{R}^{n+1}$ we introduce the following notation:

[TABLE]

so that the Laplace equation in $\mathbb{R}^{n+1}$ becomes

[TABLE]

Theorem 2**.**

Suppose $U$ is an open set in $\mathbb{R}^{n+1}$ such that a portion of its boundary $\Omega\subset\{x_{n+1}=0\}\cap\partial U$ is an open set in $\mathbb{R}^{n}$. Let $u$ be a harmonic function in $U\cup\Omega$. Denote by $Du$ the tangential component of the full gradient of $u$ on $\Omega$, i.e., $Du=(D_{1}u,\dots,D_{n}u)$. Let

[TABLE]

Then $u$ satisfies

[TABLE]

in $\Omega$, where ${\rm div\,}$is the divergence operator in $\mathbb{R}^{n}$ and $\sigma={\rm sgn\,}D_{n+1}u$.

Proof.

Taking partial derivatives of the first equation in (10), we obtain

[TABLE]

[TABLE]

where, here and throughout, the summation from 1 to $n$ is assumed for repeated indices $i$ and $j$. Multiplying the last equation by $D_{n+1}u$ and replacing in it the terms involving the partial derivative with respect to $x_{n+1}$ using (12) and (9) we arrive at

[TABLE]

Rearranging the terms further, we obtain

[TABLE]

the equation involving the partial derivatives with respect to $x_{1},\dots,x_{n}$ only. Finally, equation (11) is the divergence form of equation (13). ∎

Equation (11) is elliptic but not uniformly elliptic: the eigenvalues are $\lambda=p-|Du|^{2}$ and $\Lambda=p$. In particular, the ellipticity is lost where the full gradient of $u$ is tangent to the boundary, i.e., on the magnetic equator as referred to in [4].

The equation has an outward resemblance to the minimal surface equation

[TABLE]

in the divergence part of the operator. Also equation (11) is similar to the equation of gas dynamics

[TABLE]

in that the former loses ellipticity when $|Du|^{2}=p$ and the latter changes the type from elliptic when $|Du|<[2/(\gamma+1)]^{1/2}$ to hyperbolic when $[2/(\gamma+1)]^{1/2}<|Du|<[2/(\gamma-1)]^{1/2}$.

4. Numerical studies

In our numerical experiments, we consider equation (11) on the unit square in $\mathbb{R}^{2}$. The data is generated from functions of the form (5) that are harmonic in $\mathbb{R}^{3}$ with the exception of singularities, if any. Specifically, we choose the function $u$ from among the following choices:

• •

$u_{0}(x,y,z)=(x-x_{0})^{2}+(y-y_{0})^{2}-2(z-z_{0})^{2}$, where $(x_{0},y_{0},z_{0})=(-2,-3,-2.5)$;

• •

$u_{1}(x,y,z)=\Gamma(x-x_{1},y-y_{1},z-z_{1})$, where

[TABLE]

is the fundamental solution of the Laplace equation in $\mathbb{R}^{3}$ and $(x_{1},y_{1},z_{1})=(0.2,0.1,0.5)$;

• •

$u_{0}+u_{1}$.

The function $u_{0}$ plays the role of a potential for the background field; the function $u_{1}$ is a monopole and it models a small perturbation; and $u_{0}+u_{1}$ generates the combined field.

The values of the data: $p$, $q$, $\sigma$, and of the exact solution (for boundary values and error estimates) are collected on the unit square embedded in the coordinate $xy$-plane, i.e., $\Omega=[0,1]\times[0,1]\times\{z=0\}$. We supplement equation (11) with the Dirichlet boundary conditions. Although this choice of data requires the knowledge of $u$ on $\partial\Omega$ and it is not realistic in practice, we use it to take advantage of the available boundary value solver. In this sense, the function $\sigma$ is another exogenous assumption on the data. In all of our examples, $\partial u/\partial\nu$ is sign definite, i.e., $\sigma$ is constant.

We plot the function $u$ and corresponding $|Du|^{2}$ and $(\partial/\partial\nu)|Du|^{2}$ in Figure 1. We observe that the range of values of both the potential $u$ and the magnitude of the gradient $|Du|^{2}$ in the perturbation is two orders of magnitude smaller than those in the background field. The plots of these functions in the background field and in the combined field are virtually indistinguishable by visual inspection. The effect of the perturbation is more pronounced however for $(\partial/\partial\nu)|Du|^{2}$. This function is constant in the background field and has a local extremum near the origin in the perturbation. The combined field inherits the magnitude from the background field and the shape from the perturbation.

We apply the finite element method as it is implemented in the PDE Toolbox of Matlab. Being in divergence form, equation  (11) admits a weak formulation:

[TABLE]

This equation is discretized using the linear elements. We compute the solution for the values of the target maximum mesh edge length of the triangular elements $H_{max}$ ranging from 0.025 to 0.2 . The stock solver in the PDE Toolbox uses the Gauss-Newton iteration which requires an initial guess. On the coarsest mesh, we solve the two-dimensional Laplace equation with the same Dirichlet boundary values as in the non-linear problem and supply the resulting harmonic function as an initial guess to the non-linear solver. For successive mesh refinements, the initial guess is an interpolation of the solution from the previous mesh. The estimates of the $L^{2}$- and $H^{1}$-norms of the errors are based on the comparison with the exact solution. They are summarized along with the estimated rates of convergence in Table 1. The log-log plots of these errors appear in Figure 2.

Conclusions

The Backus problem is not well-posed and requires extra conditions to guarantee uniqueness. We consider this problem with additional data that is in the same paradigm as the non-linear boundary values in the original formulation of the problem. We show how the expanded data may be used to estimate the number of sources in the corresponding inverse source problem on the plane. In addition, we derive a non-linear equation in terms of the augmented data that is satisfied by the harmonic function on a manifold of lower dimension. We demonstrate that the resulting equation is amenable to numerical solution by standard methods and we recover some solutions using commonly available software.