Schrödingerisation based computationally stable algorithms for ill-posed problems in partial differential equations

In this paper, we are interested in numerically computing ill-posed problems that follow the evolution of a general dynamical system

whose input data given by ${\bm{u}}_{0}\in\mathbb{R}^{n}$. For simplicity, we first consider $H$ being Hermitian with $n$ real eigenvalues, ordered as

We assume that $\lambda_{i}(H)>0$, for some $1\leq i\leq n$, which implies that the system (1.1) contains unstable modes, thus the initial value problem is ill-posed since its solution will grow exponentially in time. This causes significant computational challenges. Normally the numerical errors will grow exponentially unless special care is taken.

Ill-posed or unstable problems appear in many physical applications, for example, fluid dynamics instabilities such as Rayleigh-Taylor and Kevin-Helmholtz instabilities [21, 10], plasma instability [12], and Maxwell equations with negative index of refraction [30, 28]. It also appears in inverse problems [34, 4]. One of the most classical ill-posed problems is the backward heat equation which suffers from the catastrophic Hadamard instability. Usually, some regularization technique is used to make the problems well-posed [8].

In this paper, we study a generic yet simple computational strategy to numerically compute ill-posed problems based on Schrödingerisation, introduced for quantum simulation for dynamical systems whose evolution operator is non-unitary [16, 17]. The idea is to map it to one higher dimension, into Schrödinger-type equations, that obeys unitary dynamics and thereby naturally fitting for quantum simulation. Since the Schrödinger-type equations are Hamiltonian systems that are time-reversible, they can be solved both forwards and backwards in time in a computationally stable way. This makes it suitable for solving unstable problems, as was proposed in our previous work [15]. In this article, we study this method for the classical backward heat equation, and a linear convection equation with imaginary wave speed (or negative index of refraction).

The initial-value problem to the backward heat equation is ill-posed in all three ways: (i) the solution does not necessarily exist; (ii) if the solution exists, it is not necessarily unique; (iii) there is no continuous dependence of the solution on arbitrary input data [27, 18, 13, 23]. This problem is well-posed for final data whose Fourier spectrum has compact support [26]. However, even when the solution exists and is unique, computing the solution is difficult since unstable physical systems usually lead to unstable numerical methods. There have been various treatments of this difficult problem, for example quasi-reversibility methods [35, 32, 22, 3, 25], regularization methods [11, 19, 32, 6], Fourier truncation methods [29], etc. Our approach computes solution consistent with the Fourier truncation method.

The actual implementation of our backward heat equation solver is as follows. First we start with the Fourier transform of the input data denoted by $u_{T}$ and truncate the Fourier mode in finite domain. This is either achieved automatically if the Fourier mode of $u_{T}$ has compact support, or we choose a sufficiently large domain in the Fourier space such that outside it the Fourier coefficient of $u_{T}$ is sufficiently small. The latter can usually be done since the forward heat equation gives a solution that is smooth for $t>0$ so its Fourier coefficient decays rapidly. The Schrödingerisation technique lifts the backward heat equation to a Schrödinger equation in one higher dimension that is time reversible, which can be solved backward in time by any reasonable stable numerical approximation. We then recover $u(t)$ for $0\leq t<T$ by integrating or pointwise evaluation over suitably chosen domain in the extended variable. Since the time evolution is based on solving the Schrödinger equation, the computational method is stable.

We point out that the truncation in the Fourier space regularizes the original ill-posed problem. Although the initial-value problem to the Schrödingerised equation is well-posed even without this Fourier truncation, to recover $u$ one needs finite Fourier mode. We also show that the strategy also applies to variable coefficient heat equation.

We will also apply the same strategy to solve the unstable linear convection equation that has imaginary wave speed.

We will give error estimates on these methods and conduct numerical methods that verify the results of the error analysis. The methods work for both classical and quantum computers. Since the original Schrödingerisation method was introduced for quantum simulation of general PDEs, we will also give the quantum implementation of the computational method.

Moreover, we introduce a smooth initialization for the Schrödingerised equation, so the initial data will be in $C^{k}(\mathbb{R})$ in the extended space for any integer $k\geq 1$ (see section 4.3). This will lead to essentially the spectral accuracy for the approximation in the extended space, if a spectral method is used. Consequently, the extra qubits needed due to the extra dimension, if a qubit based quantum algorithm is used, for both well-posed and ill-posed problems, becomes almost $\log\log{1/\varepsilon}$ where $\varepsilon$ is the desired precision. This reduces the complexity of Schrödingerisation based quantum algorithms introduced in [16, 17] for any non-unitary dynamical system.

The rest of the paper is organized as follows. In section 2, we give a brief review of the Schrödingerisation approach for general linear ODEs. In section 3, we study the backward heat equation, for both constant and variable coefficient cases, and show the approximate solution based on the framework of Schrödingerisation. In section 4, the numerical method and error analysis of the backward heat equation are presented. In addition, a $C^{k}$ smooth initial data with respect to the extended variable is constructed to improve the convergence rates in the extended space. In section 5, we apply this technique to the convection equations with purely imaginary wave speed. In section 6, we show the numerical tests to verify our theories. Section 7 shows the quantum algorithms and the corresponding complexity.

Throughout the paper, we restrict the simulation to a finite time interval $t\in[0,T]$. The notation $f\lesssim g$ stands for $f\leq Cg$ where $C$ is independent of the mesh size and time step. Scalar-valued quantities and vector-valued quantities are denoted by normal symbols and boldface symbols, respectively. Moreover, we use a 0-based indexing, i.e. $j=\{0,1,\cdots,N-1\}$, or $j\in[N]$, and ${\bm{e}}_{j}^{(N)}\in\mathbb{R}^{N}$, denotes a vector with the $j$-th component being $1$ and others $0$. We shall denote the identity matrix and null matrix by $I$ and 0, respectively, and the dimensions of these matrices should be clear from the context, otherwise, the notation $I_{N}$ stands for the $N$-dimensional identity matrix.

We start with the basic framework set up in [15] for forward problems which we first briefly review here.

Using the warped phase transformation ${\bm{w}}(t,p)=e^{-p}{\bm{u}}$ for $p>0$ and symmetrically extending the initial data to $p<0$, one has the following system of linear convection equations [16, 17]:

This is clearly a hyperbolic system. When the eigenvalues of $H$ are all negative, the convection term of (2.1) corresponds to waves moving from the right to the left. One does, however, need a boundary condition on the right-hand side. Since ${\bm{w}}$ decays exponentially in $p$, one just needs to select $p=p^{R}$ large enough, so $w$ at this point is essentially zero and a zero incoming boundary condition can be used at $p=p^{R}$. Then ${\bm{u}}$ can be numerically recovered via

or

However, when $\lambda_{i}(H)>0$ for some $i$, some of the waves evolving through (2.1) will instead propagate from left to right. If we were solving the problem in domain $p\in[0,\infty)$, we would need a boundary condition at $p=0$, which is not possible. Instead, we extend the computational domain to also include $[-p^{L},0)$ with a zero incoming boundary at $p=-p^{L}$ for $p^{L}>0$ sufficiently large (so $e^{-p^{L}}\approx 0$). In summary, we use zero boundary condition for (2.1) restricted on the finite domain $[-p^{L},p^{R}]$. It is worth noting that extending to the $p<0$ domain is done for the convenience of employing the Fourier spectral method, which requires periodic boundary conditions, and our zero boundary conditions can be regarded, approximately, as periodic boundary condition. It is important to note that any waves–corresponding to positive eigenvalues of $H$– that start from $p<0$ and travel to the right will induce spurious solutions to the domain $p>0$, therefore, when recovering $u$, one needs to select the correct domain to avoid using these spurious waves. Therefore, the following theorem [15] must be utilized.

Since (2.1) is a hyperbolic system, thus the initial value problem is well-posed and one can solve it numerically in a stable way (as long as suitable numerical stability condition–the CFL condition–is satisfied). Thus although the original system (1.1) is ill-posed, we can still solve the well-posed problem (2.1), and then recover ${\bm{u}}$ using Theorem 2.1, as long as $\lambda_{n}(H)<\infty$!

We also point out here that the discretization of the $p$-derivative in Eq. (2.1) can be achieved not only by using Fourier spectral methods but also through other numerical schemes, such as finite difference or finite element methods. In such cases one can just solve the problem in $p>0$ and, if $\lambda_{i}(H)>0$ for some $i$, one needs to supply boundary condition at $p=0$ for the characteristic fields corresponding to positive eigenvalues of $H$. This boundary data can be set arbitrarily, since the induced spurious wave propagating with speed $\lambda_{i}(H)>0$ will not be used when we recover $u$ using data from $p>p^{\diamond}$, as laid out by Theorem 2.1.

In this section, we consider the following one-dimensional backward heat equation in the infinite domain,

where $\Omega_{x}=(-\infty,\infty)$. We want to determine $u(t,\cdot)$ for $T>t\geq 0$ from the data $u_{T}$. The change of variable $v(t,x)=u(T-t,x)$ leads to the following formulation of (3.1),

It is easy to see that the eigenvalues of $H$ in (1.1) are positive after spatial discretization. Let $\hat{g}(\eta)$ denote the Fourier transform of $g(x)\in L^{2}(\Omega_{x})$ and define it by

Let $\|g\|_{H^{s}(\Omega_{x})}$ denote the Sobolev norm $H^{s}(\Omega_{x})$ defined as

When $s=0$, $\|\cdot\|_{H^{0}(\Omega_{x})}$ denotes the $L^{2}(\Omega_{x})$-norm, and $s$ can be noninteger [9]. For the above solution to be well-defined, we make the following assumptions of $u_{T}$ for analysis:

• (H1)

  Assume ${u}_{T}\in H^{s}_{0}$, namely $H^{s}$ with compact support. This yields a well-posed problem in a subspace of $L^{2}$ in the sense of Hadamard [26].

• (H2)

  For more general $\hat{u}_{T}$, we will truncate the domain in $\eta$ and begin with the truncated–thus compactly supported–final data.

We remark here that truncation is inherently achieved in our methodology (see (S1) and (S2)), thus providing a form of regularization to the original ill-posed problem.

If the solution of (3.1) in this Sobolev space exists, then it must be unique [9]. We assume $u(t,x)$ is the unique solution of (3.1). Applying the Fourier transform, one gets the exact solution $u(t,x)$ of problem (3.1):

Denoting $u_{0}(x)=u(0,x)$, then it is easy to see from (3.3) that

The Schrödingerisation for (3.1) gives

where $\Omega_{p}=(-\infty,\infty)$. Again using the Fourier transform technique to (3.5) with respect to the variables $p$ and $x$, one gets the Fourier transform $\hat{w}(t,\eta,\xi)$ of the exact solution $w(t,x,p)$ of problem (3.5) satisfying

The solution of $w$ is then found to be

Note the initial-value problem (3.5) is well-posed, even for $u_{T}\in H^{s}$ not being compactly supported in the Fourier space, as can be easily seen from (3.6), which is an oscillatory ODE with purely imaginary spectra. Following the proof in [9, Theorem 5, section 7.3], one has the regularity of $w(0,x,p)\in H^{s}(\Omega_{x})\times H^{1}(\Omega_{p})$ under the assumption in (H1) or (H2).

Using the Fourier transform only on $x$ gives

where $\hat{w}^{x}$ is the Fourier transform of $w$ on the $x$-variable, i.e.

Here $\hat{w}^{x}$ is a linear wave moving to the right, if one starts from $t=T$ and going backward in time (one can also use the change of variable $t\mapsto T-t$ to make the problem (3.5) forward in time for notation comfort), and the solution is computed by the inverse Fourier transformation

This corresponds to the case of (2.1) in which all eigenvalues of $H$ are positive. Therefore, to recover $u(0,x)$ from Eq.(3.10), one needs to choose $p>\eta^{2}T$ for time $T>0$ to obtain

This is the basis for the introduction of our stable computational method: we solve the $w$-equation (3.5) with a suitable–stable– computational method, and then recover $u$ by either integration or pointwise evaluation using Theorem 2.1.

However, in the continuous space $\eta$ maybe unbounded, thus the condition $p>\eta^{2}T$ may not be satisfied. We consider two scenarios here:

1. (S1)

   The input data $\hat{u}_{T}$ has compact support, namely $\hat{u}_{T}=0$ for $|\eta|>\eta_{0}$ for some $\eta_{0}<\infty$. While this is not generally true, when $u_{T}(x)$ is smooth in $x$, its Fourier mode decays rapidly, so one can truncate $\eta$ for $|\eta|>\eta_{0}$ with desired accuracy.

2. (S2)

   One discretizes equation (3.5) numerically. This usually requires first to truncate the $x$-domain so it is finite, followed by some numerical discretization of the $x-$derivative, for example the finite difference or spectral method. Then $|\eta|\leq\eta_{0}\leq\mathscr{O}(1/(\Delta x)^{2})$ where $\Delta x$ is the spatial mesh size in $x$.

In this section, we consider the discretization of (3.5) and the corresponding error estimates. For simplicity, we consider $u$ to be a periodic function defined in $\Omega_{x}=[-\pi,\pi]$.