A Carleman-based numerical method for quasilinear elliptic equations with over-determined boundary data and applications

The main aim of this paper is to develop a numerical method based on Carleman estimates to solve quasilinear elliptic PDEs with over-determined boundary data. We consider this new method as the second generation of Carleman-based numerical methods while the first generation is called the convexification, which will be mentioned in detail later. Let $\Omega$ be an open and bounded domain in $\mathbb{R}^{d}$, $d\geq 2$, with smooth boundary $\partial\Omega$. Let $f$ and $g$ be two smooth functions defined on $\partial\Omega$. Let $F:\overline{\Omega}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}$ be a function in the class $C^{2}.$ Let $A=(a_{ij})_{i,j=1}^{d}:\overline{\Omega}\to\mathbb{R}^{d\times d}$ be $C^{2}$, symmetric, and positive definite, that is,

for some fixed $\gamma\in(0,1)$. The following problem is of our interests.

Problem 1.1 is motivated by a class of nonlinear inverse problems in PDEs, in which $f$ and $g$ are the data that can be measured. One important goal of inverse problems is to reconstruct the internal structure of a domain from boundary measurements, which allow us to impose both Dirichlet and Neumann data of the unknown in (1.1). Recently, a unified framework to solve such inverse problems was developed by the research group of the first and second authors, which has two main steps. In the first step, by introducing a change of variables, one derives a PDE of the form (1.1) from the given inverse problem, in which $f$ and $g$ can be computed directly by the given boundary data. In the second step, one numerically solves (1.1) to find $u^{*}$. The knowledge of $u^{*}$ directly yields that of the solution to the corresponding inverse problem under consideration. See [15, 30] and the references therein for some works in this framework. Moreover, this unified framework was successfully tested with experimental data in [13, 14, 22]. Another motivation to study Problem 1.1 is to seek solutions to Hamilton-Jacobi equations under the circumstance that the Neumann data of the unknown can be computed by its Dirichlet data and the given form of the Hamiltonian, see e.g., [25, Assumption 1.1 and Remark 1.1]. Since inverse problems are out of the scope of this paper, we only focus on the applications in solving quasilinear elliptic PDEs and first-order Hamilton-Jacobi equations.

A natural approach to solve (1.1) is based on optimization. That means one sets the computed solution to (1.1) as a minimizer of a mismatch functional, e.g.,

subject to the Cauchy boundary conditions $v|_{\partial\Omega}=f$ and $\partial_{\nu}v|_{\partial\Omega}=g.$ The methods based on optimization are widely used in the scientific community, especially in computational mathematics, physics and engineering. Although effective and popular, the optimization-based approaches have some drawbacks:

1. 1.

   In general, it is not clear that the obtained minimizer approximates the true solution to (1.1).

2. 2.

   The mismatch functional $J$ is not convex, and it might have multiple minima and ravines (see an example in [43] for illustration). To deliver reliable numerical solutions, one must know some good initial guesses of the true solutions.

3. 3.

   The computation is expensive and time consuming.

Drawbacks # 1 and #2 can be treated by the convexification method, which is designed to globalize the optimization methods. The main idea of the convexification method is to employ some suitable Carleman weight functions to convexify the mismatch functionals. The convexity of weighted mismatch functionals is rigorously proved by Carleman estimates. Several versions of the convexification method have been developed in [1, 15, 16, 17, 19, 22] since it was first introduced in [21]. Moreover, we recently discovered that the convexification method can be used to solve a large class of first-order Hamilton-Jacobi equations [25].

In this paper, we introduce a new method to solve (1.1) based on linearization and Carleman estimates. Like the convexification method, our method delivers a reliable solution to (1.1) without requiring a good initial guess. This fact is rigorously proved. Unlike the convexification method which is time consuming, our new method quickly provides the desired solutions. Its converge rate is $O(\theta^{n})$ as $n\to\infty$ for some $\theta\in(0,1).$

We find the numerical solution to (1.1) by repeatedly solving the linearization of (1.1) by a new “Carleman weighted” quasi-reversibility method. The classical quasi-reversibility method was first proposed in [28], and it has been studied intensively since then (see [18] for a survey). By a Carleman weighted quasi-reversibility method, we mean that we let a suitable Carleman weight function involve in the cost functional suggested by the classical quasi-reversibility method. The presence of the Carleman weight function is the key for us to prove our convergence theorem. Our process to solve Problem 1.1 is as follows. We first choose any initial solution that might be far away from the true one. Denote this initial solution by the function $u_{0}$. Linearizing (1.1) about $u_{0}$, we obtain a linear PDE. We then solve this linear PDE by the Carleman weighted quasi-reversibility method to obtain an updated solution $u_{1}$. Using the Carleman weighted quasi-reversibility method rather than the classical one in this step is the key to our success. By iteration, we repeat this step to construct a sequence $\{u_{n}\}_{n\geq 0}$. The convergence of this sequence to the true solution to (1.1) is proved by using Carleman estimates. We then apply this method to numerically solve some quasilinear elliptic equations. It is important to note that our approach works well for systems of quasilinear elliptic PDEs too.

Next, we use our method to solve some first-order Hamilton-Jacobi equations. More precisely, to find viscosity solutions to the first-order equation

we use the vanishing viscosity procedure and consider, for $0<\epsilon_{0}\ll 1$,

with given Cauchy boundary data. Our new method is robust in the sense that it works for general nonlinearity $F({\bf x},u,\nabla u)$ that might not be convex in $\nabla u.$ We refer the readers to [9, 8, 47] and the references therein for the theory of viscosity solutions. A weakness of our new approach in computing the viscosity solutions to Hamilton-Jacobi equations is that we need to require both Dirichlet and Neumann data of the unknown $u$. See [25, Remark 1.1] for some circumstances that this requirement is fulfilled. There have been many important methods to solve Hamilton-Jacobi equations in the literature. For finite difference monotone and consistent schemes of first-order equations and applications, see [2, 10, 40, 44, 45] for details and recent developments. If $F=F({\bf x},u,\nabla u)$ is convex in $\nabla u$ and satisfies appropriate conditions, it is possible to construct some semi-Lagrangian approximations by the discretization of the Dynamical Programming Principle associated to the problem, see [11, 12] and the references therein.

The paper is organized as follows. In Section 2, we recall a Carleman estimate and two examples about inverse problems in which Problem 1.1 appears. In Section 3, we introduce the new iterative method based on linearization and Carleman estimates. In Section 4, we prove the convergence of our method. Some numerical results for quasilinear elliptic equations and first-order Hamilton-Jacobi equations are presented in Section 5. Concluding remarks are given in Section 6.

In this section, we recall a Carleman estimate, which plays a key role for the proof of the convergence theorem in this paper. We then present an inverse scattering problem in which Problem 1.1 appears.

Our approach to solve (1.1) is based on linearization and iteration. Assume that the solution $u^{*}$ to (1.1) is in the space $H^{p}(\Omega)$ for some $p>\lceil d/2\rceil+2$ where $\lceil d/2\rceil$ is the smallest integer that is greater than $d/2.$ Define the set of admissible solutions

Then, the assumption in Problem 1.1 implies that $W\neq\emptyset$, and $u^{*}\in W$. We now construct a sequence $\{u_{n}\}_{n\geq 0}$ that converges to the solution $u^{*}$ to (1.1). Take a function $u_{0}\in W.$ Assume by induction that we have the knowledge of $u_{n}$ for some $n\geq 0$. We find $u_{n+1}$ as follows. Assume that $u^{*}=u_{n}+h$ for some $h\in H_{0}$ where

Plugging $u^{*}=u_{n}+h$ into (1.1), we have

for all ${\bf x}\in\Omega.$ Heuristically, we assume at this moment that $h$ is small, that is, $\|h\|_{C^{1}(\Omega)}\ll 1$.

By Taylor’s expansion, we approximate (3.2) as

for all ${\bf x}\in\Omega$ where

for all $v\in H^{p}(\Omega)$. Here, $F_{s}$ and $\nabla_{{\bf p}}F$ are the partial derivative of $F$ with respect to its second variable and its gradient vector with respect to the third variable, respectively.

The next step is to compute a function $h\in H_{0}$ satisfying (3.3). Since there is no guarantee for the existence of such a function $h$, we only compute a “best fit” function $h$ by the Carleman-based quasi-reversibility method described below. For each $v\in H^{p}(\Omega)$, $\lambda>1$, $\beta>0$ and $\eta>0$, define the functional $J_{v}^{\lambda,\beta,\eta}:H_{0}\to\mathbb{R}$ as

Here, the function $\mu_{\beta}$ is defined in (2.1) and $\eta\|v+\varphi\|^{2}_{H^{p}(\Omega)}$ is a regularization term.

For each $n\geq 0$, thanks for Proposition 3.1, we can minimize $J_{u_{n}}^{\lambda,\beta,\eta}(\varphi)$ on $H_{0}$. The unique minimizer is the desired function $h.$ We then set

The construction of the sequence $\{u_{n}\}_{n\geq 0}$ above is summarized in Algorithm 1. We will prove that the sequence $\{u_{n}\}_{n\geq 0}$ converges to $u^{*}$ in Section 4 as $n\to\infty$ and $\eta\to 0$. The presence of the Carleman weight $e^{2\lambda\mu_{\beta}({\bf x})}$ is a key point for us to prove this convergence result.

In this section, we prove that the sequence $\{u_{n}\}_{n\geq 0}$ generated by Algorithm 1 converges to the solution $u^{*}$ to (1.1). The following result is the main theorem in this paper.